nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。

任意の基底関数における正規方程式

多項式による正規方程式の導出を、こちらにサイトでお世話になった。
線形回帰の Normal Equation(正規方程式)について
しかし、別に多項式じゃなくても何でも良いと思ったので、導出してみた。

何かしらの入力 xに対する出力 y N個の組 (x_1, y_1), \cdots, (x_N, y_N) で得られているとする。
この結果から、 (x,y)の関係を基底関数 \phi^{(0)}(x), \cdots, \phi^{(N_P)}(x)の組み合わせ(線形結合)で表現したい。
言い換えれば、各基底関数に掛ける係数 c^{(0)}, \cdots, c^{(N_p)}をどのように取れば最も誤差が小さくなるか、を求めたい。
ここで、関数 \phi xに依存するが、係数 c xに依存しないことに注意。全ての xで上手いこと行くような cが欲しいのである。

ある点 (x_i, y_i)における誤差 \Delta_iは、差分の二乗を取って大小の比較が出来るように定義する。(要素数で割って分散で表しても良い)


\displaystyle
\Delta_i = \left( c^{(0)} \phi^{(0)}(x_i) + \cdots + c^{(N_p)} \phi^{(N_p)}(x_i) - y_i \right)^2 \\
\displaystyle
\qquad = \left( \left( \sum^{N_p}_{j=0} c^{(j)} \phi^{(j)}(x_i) \right) - y_i  \right)^2 \\
\displaystyle
\qquad = \left( \left( X \vec{c} \right)_i  - y_i \right)^2 =  \left( X \vec{c} \right)^2_i - 2 \left( X \vec{c} \right)_i y_i  + y^2_i \\
\displaystyle
\qquad \equiv \delta^2_i

ここで、行列 Xとベクトル \vec{c}, \vec{y}, \vec{\delta}を以下の様に定義した。


\displaystyle
\begin{align}
X &= \left( \begin{array}{cccc} \phi^{(0)}(x_1) & \phi^{(1)}(x_1) & \cdots & \phi^{(N_p)}(x_1) \\ \phi^{(0)}(x_2) & \phi^{(1)}(x_2) & \cdots & \phi^{(N_p)}(x_2) \\ \vdots & \vdots & \vdots & \vdots \\ \phi^{(0)}( x_N ) & \phi^{(1)}( x_N ) & \cdots & \phi^{(N_p)}( x_N ) \end{array} \right) \\
\vec{c} &= \left( \begin{array}{c} c^{(0)} \\ c^{(1)} \\ \vdots \\ c^{(N_p)} \end{array} \right)
\end{align}

\displaystyle
\begin{align}
\vec{ y } = \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_N \end{array} \right)
\quad
\vec{ \delta } = \left( \begin{array}{c} \delta_1 \\ \delta_2 \\ \vdots \\ \delta_N \end{array} \right)
\end{align}

 \vec{c}はベクトルのなので、 X \vec{c}もベクトルであることに注意。
また、行列やベクトルのサイズが N_p+1なのか Nなのかにも注意。

行列を使って表現することを意識して、全誤差 \Deltaを次の様に表す。

 
\displaystyle
\Delta = \sum^N_{i=1} \Delta_i = \sum^N_{i=1} \delta_i \delta_i = {}^t \! \vec{\delta} \, \vec{\delta}
 = {}^t \!\! \left(  X \vec{c} -\vec{ y } \right) \left( X \vec{c} - \vec{ y } \right)

ここで、 \Delta c^{(k)}微分したものを \Delta^{(k)}と定義すると、


\displaystyle
\frac{d}{d c^{(k)}} \left( X \vec{c}\right)_i = \frac{d}{d c^{(k)}} \sum^{N_p}_{j} c^{(j)} \phi^{(j)}(x_i) \\
\displaystyle
\quad = \phi^{(k)}( x_i )
\\
\displaystyle
\frac{d}{d c^{(k)}} \left( X \vec{c}\right)^2_i = \frac{d}{d c^{(k)}} \sum^{N_p}_{j} \sum^{N_p}_{j'} c^{(j)} c^{(j')} \phi^{(j)}(x_i) \phi^{(j')}(x_i) \\
\displaystyle
\quad = \phi^{(k)}( x_i ) \sum^{N_p}_{j'} c^{(j')} \phi^{(j')}(x_i) + \phi^{(k)}( x_i ) \sum^{N_p}_{j} c^{(j)} \phi^{(j)}(x_i) \\
\displaystyle
\quad = 2 \phi^{(k)}( x_i ) \sum^{N_p}_{j} c^{(j)} \phi^{(j)}(x_i) = 2 \phi^{(k)}( x_i ) \left( X \vec{c}\right)_i
\\
\displaystyle
\therefore \Delta^{(k)}_i = 2 \phi^{(k)}( x_i ) \left( \left( X \vec{c}\right)_i - y_i \right) = 2 \phi^{(k)}( x_i ) \delta_i
\\
\displaystyle
\therefore \Delta^{(k)} =  2 \sum^N_i \phi^{(k)}( x_i ) \delta_i = 2 \left( {}^t \! X \vec{\delta} \right)_k

ただし、 \vec{ \Delta }の要素数は係数 cと同じであることに注意。

今考えているのは、全誤差 \Deltaが最小になる様な \vec{c}を求めることなので、 \Delta極値を取る条件、すなわち全ての \Delta^{(k)}がゼロになる条件から、係数 cが満たすべき方程式が得られる。


\displaystyle
\vec{ \Delta } =  2 {}^t \! X \vec{\delta} = 2 {}^t \! X \left( X \vec{c} - \vec{ y } \right) = \vec{0} \\
\displaystyle
\therefore \vec{ c } = \left( {}^t \! X X \right)^{-1} \, {}^t \! X \, \vec{ y }

これが、正規方程式と呼ばれる関係式であり、結果だけ見ると、


\displaystyle
\vec{y} = X \vec{c}

を満たす \vec{c}を求める問題にあたかもすり替わったかの様に見える。
実際には、極値問題をちゃんと解いているわけで、「えいやっ!と置いてしまえ!」ということではなく、安心して使える。

導出的には最小値ではなく極値を求めているため、式だけでは停留点や最大値が求まってしまう可能性を除けないが、実際には誤差は無限に大きくなれるため、少なくとも最大値は除かれている。

この導出では、基底関数の形を特に指定していないため、(個人的には)見通しが良くなっている(と思う)。
例えば、二つの基底関数 \phi^{(0)}( x ) = 1,  \phi^{(1)}( x ) = x を使えば線形フィッティングだし、様々な周期の三角関数(もしくは平面波)を使えばフーリエ変換になる。

また、面白いのが、基底関数行列 Xは上で示したように N \times (N_p + 1 )行列なので、一般に正方行列ではなく、正則でもない。
そのため、 \vec{y} = X \vec{c} \vec{c}について解くためには、まずは {}^t \! Xを掛けて \vec{c}の係数行列を正方にするという着眼点が、個人的には新鮮だった。

散乱振幅からサイトT行列を復元する。

実験結果か何かで、散乱振幅 f( \theta )の角度分布が得られているとする。
ここからサイト T行列を復元する方法を考える。
サイト T行列が求まれば、そこからArgand diagramを作って、(準)共鳴準位があるとかないとか、追加で情報が得られる(もし角度分布のエネルギー依存性まで測っていればの話だが、、、)。
サイト T行列のArgand diagramは前回紹介した方法が使える。
koideforest.hatenadiary.com

散乱振幅は、以下の形で与えられる。

\displaystyle
f( \theta ) = \sum_l ( 2 l + 1 ) t_l P_l( \cos \theta )

 P_l( \cos \theta ) はLegendre多項式である。
Legendre多項式には、直交性が成り立つ。それは以前の記事で確認した。
koideforest.hatenadiary.com

したがって、


\displaystyle
\int^1_{-1} d( \cos \theta ) f( \theta ) P_l( \cos \theta )
  = \sum_{l'} ( 2 l' + 1 ) t_{l'} \int^1_{-1} d( \cos \theta ) P_{l'}( \cos \theta ) P_{l}( \cos \theta )  
\\
\displaystyle
  = \sum_{l'} ( 2 l' + 1 ) t_{l'} \frac{ 2 }{ 2 l + 1 } \delta_{ l l' }
  = 2 t_l
\\
\displaystyle
\therefore
t_l = \frac{ 1 }{ 2 } \int^1_{-1} d( \cos \theta ) f( \theta ) P_l( \cos \theta )
  = \frac{ 1 }{ 2 } \int^1_{-1} dx \, f( \cos^{-1}( x ) ) \, P_l( x )

例えば、角運動量 l = 2のサイト T行列 t_2のみからなる散乱振幅の角度分布は次のようになる。
f:id:koideforest:20180831011856p:plain

当たり前であるが、思いっ切り P_2( \cos \theta ) の形である。
この時、 t_2 = 0.1234 + 0.5678 i と適当に設定したが、導出した式からちゃんと t_2が復元されることが以下のソースで確認出来る。

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import legendre
from scipy.integrate import quad

def Plx( l, x ):
    return legendre( l )( x )

cp_f = 1.
t_2 = 0.1234 + 0.5678j
f = lambda theta: cp_f * ( 2 * 2 + 1 ) * t_2 * Plx( 2, np.cos( theta ) )
thetas = np.linspace( 0., 2. * np.pi, 100 )

fig = plt.figure( figsize = ( 4.72, 4.72 ) )
ax = fig.add_subplot( 1, 1, 1, polar = True )
ax.plot( thetas, abs( f( thetas ) ) )
plt.show()
#plt.savefig( 'scattering_amplitude_2.png' )
#plt.close()

r0 = quad( lambda x: np.real( f( np.arccos( x ) ) * Plx( 0, x ) ), -1, 1 )[ 0 ] / ( 2 * cp_f )
i0 = quad( lambda x: np.imag( f( np.arccos( x ) ) * Plx( 0, x ) ), -1, 1 )[ 0 ] / ( 2 * cp_f )
print( r0, i0 )
# 0.0, 0.0

r2 = quad( lambda x: np.real( f( np.arccos( x ) ) * Plx( 2, x ) ), -1, 1 )[ 0 ] / ( 2 * cp_f )
i2 = quad( lambda x: np.imag( f( np.arccos( x ) ) * Plx( 2, x ) ), -1, 1 )[ 0 ] / ( 2 * cp_f )
print( r2, i2 )
# 0.1234, 0.5678

部分散乱振幅(サイトT行列)のアーガンド図

仕事で図を作る必要があったので、練習がてらまとめた。

位相シフトなどの概念にについては、wikipedia等を参照して頂きたい。
散乱振幅 - Wikipedia
簡単に言うと、「ポテンシャルが無い時の波は、ポテンシャルによってどれくらい変調されたか?」を表す量である。

ある部分波(角運動量 lの位相シフト \delta_lに対するサイト T行列 t_lは、次の式で与えられる。

\displaystyle
t_l = \frac{ e^{ 2 i \delta_l } - 1 }{ 2 i k }

位相シフトは、ある程度エネルギーが高くなると単調に減少する。
ここでは簡略化し、一次で減少するとする。
f:id:koideforest:20180829065145p:plain
この時の t_lの実部(赤点)と虚部(青点)、及びそのArgand diagramは次の図のようになる。
f:id:koideforest:20180829065152p:plainf:id:koideforest:20180829065202p:plain
左巻きでグルグルする。

仮に、一次で単調に増加するようにすると、右巻きでグルグルする。
向きは後で大事になる。

次に、(準)共鳴準位がある場合を考える。
共鳴準位のエネルギー位置で位相は急激に変化し、最終的に \piずれる。
これを簡単に \tan^{-1}で表すとする。ただし、原点を通る様に下駄を履かせている。
f:id:koideforest:20180829074635p:plainf:id:koideforest:20180829074644p:plainf:id:koideforest:20180829074651p:plain
Argand diagramにおいて、共鳴準位のある20 eV近傍で、急激に(分かりにくいが)右巻きに回転しているのがわかる。

最後に、両方を組み合わせた位相シフトの変化を考える。
f:id:koideforest:20180829074709p:plainf:id:koideforest:20180829074721p:plainf:id:koideforest:20180829074729p:plain
急に右巻きの成分が現れ、共鳴準位のから遠ざかるとまた左巻きに戻っているのがわかるかと思う。
この様な挙動から、共鳴準位の有無を観察することが出来る。

以下ソースコード

import numpy as np
from matplotlib import pyplot as plt

def wavenumbers( energy_ev ):
    RYD2EV = 13.605
    EV2RYD = 1. / RYD2EV
    return np.sqrt( energy * EV2RYD )

energy = np.linspace( 1., 100, 199 )

#delta = - energy / 20.
#delta = np.arctan( energy - 20 ) + np.pi / 2.
delta = - energy / 20. + np.arctan( energy - 20 ) + np.pi / 2.

plt.plot( energy, delta )
plt.xlabel( "Energy [eV]" )
plt.ylabel( "$\delta_l$" )
plt.show()
#plt.savefig( "delta.png" )
#plt.close()

t_l = ( np.exp( 2.j * delta ) - 1. ) / ( 2.j * wavenumbers( energy ) )
plt.plot( energy, t_l.real, "ro" )
plt.plot( energy, t_l.imag, "bo" )
plt.xlabel( "Energy [eV]" )
plt.ylabel( "$t_l$" )
plt.show()
#plt.savefig( "tl.png" )
#plt.close()

fig = plt.figure( figsize = ( 4.72, 4.72 ) )
ax = fig.add_subplot( 1, 1, 1 )
ax.plot( t_l.real, t_l.imag, "go" )
ax.set_xlim( - 1.1 * m_t, 1.1 * m_t )
ax.set_ylim( - 0.1 * m_t, 2.1 * m_t )
ax.set_xlabel( "Re( $t_l$ )" )
ax.set_ylabel( "Im( $t_l$ )" )
for ie, ene in enumerate( energy ):
    if ene % 10 == 0:
        ax.text( t_l.real[ie], t_l.imag[ie], "{:4.0f} eV".format( energy[ie] ) ) 
plt.show()
#plt.savefig( "arg.png" )
#plt.close()

ベイズの定理:薬物検査の難しさ

ベイズの定理のWikipediaで、応用例として薬物検査で陽性が出た時に本当に使用者でる確率が取り上げられている。

ベイズの定理
ベイズの定理 - Wikipedia

\displaystyle
P( B | A ) = \frac{ P( A | B ) P( B ) }{ P( A ) }

全確率の定理
Law of total probability - Wikipedia

\displaystyle
P( A ) = \sum_n P( A | B_n ) P( B_n )

感度:薬物使用者のうち99%が陽性になる。

\displaystyle
P( + | U ) = 99 / 100, \quad P( - | U ) = 1 / 100

特異度:薬物非使用者のうち99%が陰性になる。

\displaystyle
P( + | N ) = 1 / 100, \quad P( - | N ) = 99 / 100

全国民から一人無作為に選んだ結果が0.01%が薬物使用者である確率。

\displaystyle
P( U ) = 0.01 / 100, \quad P( N ) = 99.99 / 100

これから、「陽性が出た時にその人が薬物使用者である確率 P( U | + )」を求める。

\displaystyle
P( U | + )
  = \frac{ P( + | U ) P( U ) }{ P( + ) }
  = \frac{ P( + | U ) P( U ) }{ P( + | U ) P( U ) + P( + | N ) P( N ) }
\\
\displaystyle
  = \frac{ 0.99 \times 0.0001 }{ 0.99 \times 0.0001 + 0.01 \times 0.9999 }
  =  0.0098... \sim 0.01

つまり、陽性が出てもその人が薬物使用者である確率は1パーセントでしかない。
何の確証も無しに「こいつ薬やってるっぽい!」って言うのに比べて、当たる確率は100倍にはなっているが、全然検査としてダメじゃんって感じがするだろう。

薬をやっている人を、0.01%から0.1%に増やすと、

\displaystyle
P( U | + )
  = \frac{ 0.99 \times 0.001 }{ 0.99 \times 0.001 + 0.01 \times 0.999 }
  =  0.090...
同様に当たる確率も10倍程度される。そもそもの確率が高ければ、陽性のチェックの精度も高くなる。

このチェックの精度を上げるには、特異度が重要だとWikipediaに書いてある。
特異度を99.9%にすると( P( + | N ) = 0.001

\displaystyle
P( U | + )
  = \frac{ 0.99 \times 0.0001 }{ 0.99 \times 0.0001 + 0.001 \times 0.9999 }
  =  0.090...
精度が10倍近く改善された。
特異度が99.99%になると、精度は50%近くまで上がり、99.999%で90%を超える。

一方で感度はあまり重要ではなく、99.99%にしても精度は1%を超えずほぼ変化しない。

答えがわかっているものに対して近い反応が示せたとして、それを使って答えがわからないものをチェックするのは容易ではないのかがよく分かった。

誕生日のパラドックス:「私」か「全員」か

たまたま見つけたやつ。
【オリジナル】「誕生日が同じだった運命のふたりの漫画。」漫画/山本アリフレッド [pixiv]
30人のクラスの中で、同じ誕生日のペアが出来る確率は7割だから、誕生日が同じなのは奇跡ではない、というもの。
これについて、「『私と』ではなく『クラスの中で』と言っているから、指摘は正しい」というコメントが有って、確かにと思った次第である。

365人の集団が全てバラバラの誕生日を持っていたとする。
この時、

「無作為に一人選んだ時に、その人がある特定の誕生日を持つ確率」 = 1 / 365

であるが、

「この365人の中で、ある特定の誕生日を持つ人がいる確率」= 365 / 365 = 1

であるため、「私」なのか「集団」なのかで話が全然違ってくる。

つまり

「30人のクラスの中で、自分と同じ誕生日を持つ人がいる確率」= ( 30 - 1 ) /365

であるから、自分と同じ誕生日の人がいる確率は数パーセントである。
更にこの集団を、「自分が好きな人」というフィルターをかければ、その確率はグッと低くなる。十分奇跡と言えるだろう。

ちなみに、「30人のクラスの中で、複数の人が同じ誕生日を持つ確率」は、「全員が違う誕生日を持つ確率」の余事象を使って求めることが出来る。
同じ誕生日の二人組がいる確率について | 高校数学の美しい物語

 \displaystyle
P = 1 - \bar{ P } = 1 - \frac{ 365 - 1 }{ 365 } \frac{ 365 - 2 }{ 365 } \cdots \frac{ 365 - ( 30 - 1 ) }{ 365 } \sim 0.706

bar_p = 1.
for i in range( 30 ):
  bar_p = ( 365 - i ) / 365
p = 1. - bar_p
# p = 0.706316...

個人的には、パラドックスというよりかは説明不足のような気もする。
問題文は大事です。

結合法則と交換法則が成り立たない簡単な例

加法や乗法は、(それらの群の上では)結合の法則と交換の法則を満たす。


\displaystyle
( a + b ) + c = a + ( b + c )
\\
\displaystyle
a + b = b + a

減法や除法も加法群および乗法群に含まれてしまうため、何かこれらを満たさない別の簡単な例はないかなぁと思っていたら、パッと冪演算を思い付いた。
冪演算を \hat{} で定義すると、


\displaystyle
( a \hat{} b ) \hat{} c = ( a^b )^c = a^{ bc } \neq a \hat{} ( b \hat{} c ) = a^{ b^c }
\\
\displaystyle
a^b \neq b^a

少なくともパッと思い付く実数の範囲では、冪演算は半群すら作らない。
こう考えても(リー群とか色々思うと)指数はなかなか奥が深いもんだなぁと感じる。

第二量子化と場の演算子と波動関数

以前、第二量子化についての記事を書いたが、誤って消してしまったので、再度投稿。

任意の波動関数 \psi(x) を、基底関数 u( x )で展開し、その展開係数を cとする。
これに対応する場の演算子 \hat{ \psi }(x) は、消滅演算子 \hat{ a }を用いて、次のように書ける。


\displaystyle
\psi( x ) = \sum_{ \alpha } u_{ \alpha }( x ) \, c_{ \alpha }
\\
\displaystyle
\hat{ \psi }( x ) = \sum_{ \alpha } u_{ \alpha }( x ) \, \hat{ a }_{ \alpha }

したがって、第二量子化は展開係数を量子化したものと言える。

基底関数 u(x)は完備性により次の関係を満たす。


\displaystyle
\sum_{ \alpha } u^*_{ \alpha }( x' ) u_{ \alpha }( x )
  = \sum_{ \alpha } < \alpha | x' > < x | \alpha >
  = \sum_{ \alpha } < x | \alpha >< \alpha | x' > 
\\
\displaystyle
  = < x | x' > 
  = \delta( x - x' )

これを用いて、場の演算子 \hat{\psi}^{\dagger}(x')を真空状態 |> に作用させた波動関数を考えると、

\displaystyle
< x | \hat{ \psi }^{\dagger}( x' ) |> 
  = < x | \sum_{ \alpha } u^*_{ \alpha }( x' ) \, \hat{ a }^{\dagger}_{ \alpha } |>
\\
\displaystyle
  = \sum_{ \alpha }  u^*_{ \alpha }( x' ) < x | \alpha >
  = \sum_{ \alpha } u^*_{ \alpha }( x' ) u_{ \alpha }( x )
  = \delta( x - x' )

したがって、場の演算子 \hat{\psi}^{\dagger}(x')は、本当に位置 x'にだけ粒子の存在確率を与えることがわかる。
波動性は一切無い。

「粒子を足す」という言葉にも色々あって、「エネルギー \varepsilonの粒子を足す」とかだと波動性が発生していることになり、注意が必要だと思う。