一日一膳(当社比)

数理系の記事など。

1次元Ising模型をシミュレーションする.

$$\newcommand{\Expo}[1]{ \mathrm{exp}\Bigl(#1 \Bigr)}$$

$$\newcommand{\optim}[1]{ \underset{#1}{\mathrm{arg \ min}} }$$

目次

1次元Ising模型を解析的に求める

1次元Ising模型

Ising模型は,統計力学でよく題材にされる数理モデルです.
Nを3以上の整数として,X_0, X_1, \ldots, X_{N-1}を実数値確率変数とします.以下,パラメータ\sigma^2\ >\ 0, Jを導入します.{\sf X}=(X_0, \ldots, X_{N-1})^{\sf T}(ただし,{}^{\sf T}は転置をあらわす)のエネルギーを,
E({\sf X})=\displaystyle{\sum_{j=0}^{N-1}X_j{}^2 -2J \sum_{j=0}^{N-1}X_jX_{j+1} }
\quad \quad \  ={\sf X}^{\sf T}Q{\sf X}
で定義します,ここでQは次で定まるN \times N行列です,
Q_{ij} = \begin{cases}  1, \quad (i=j)\\  -J, \quad  (i-j \equiv \pm1 \ {\rm mod}\ N) \\ 0, \quad {\rm otherwise.}\quad \end{cases}.
また,パラメータJ結合定数と呼ぶことにします.
1次元Isingモデルを次の{\mathbb R}^Nに値をとる確率モデルで,{\sf X}=(X_0, \ldots, X_{N-1})^{\sf T}の同時分布p({\sf X}|\sigma^2,J)が次で与えられるものと定義します.
 p({\sf X}|\sigma^2,J)= \frac{1}{Z(\sigma^2)}\Expo{-\frac{1}{2 \sigma^2}E({\sf X})}

ただしZ(\sigma^2)は規格化定数で次で定義される

\displaystyle{Z(\sigma^2)=\int_{{\mathbb R}^N}d^N {\sf X}\ \Expo{-\frac{1}{2 \sigma^2} E({\sf X})}.}

J=0の場合は,X_0, \ldots, X_{N-1}は(正規)ホワイトノイズ,つまり標準正規分布からなる確率分布列であって互いに無相関であるものになります.

分配関数を求める

 この小節では,分配関数を実際に求めてみます.分配関数は,定義より

Z(\sigma {}^2)

\displaystyle{=\int_{{\mathbb R}^N}d^N {\sf X}\ \Expo{-\frac{1}{2 \sigma^2} {\sf X}^{\rm T}Q{\sf X}}}

=(2 \pi )^{N/2} (\sigma {}^2)^{N/2} (\det{Q})^{-1/2}

\det{Q}が必要となりますが,対角化によりN\geq 3のとき

\displaystyle{\det{Q}= \prod^{N-1}_{j=0}(1-2J\cos{(j \alpha)}),\quad (\alpha := \frac{2\pi}{N})}で得られます.よって,分配関数が求まりました.

1次元Ising模型の分配関数
Z(\sigma^2)
\displaystyle{=(2\pi\sigma^2)^{N/2}\prod^{N-1}_{j=0}(1-2J\cos{(j \alpha)})^{-1/2}}, 
ただし\alpha = \frac{2\pi}{N}

 分配関数が求まると,そこから様々な統計量を計算することができます.

分散\langle X^2_i \rangle

分散\langle X_i \rangle , i=0, \ldots, N-1は次のようにして求まります.

\langle X_i^2 \rangle

\displaystyle{ = -\frac{1}{N} \frac{\partial }{\partial (1/2\sigma^2)}\ln{Z(\sigma^2)}}

 \displaystyle{ = \frac{-2}{N} \frac{\partial }{\partial (1/\sigma^2)}{\frac{N\ln{\sigma^2}}{2}}}

\displaystyle{=\sigma^2}.

 

共分散\langle X_iX_{i+1} \rangle

共分散\langle X_i X_{i+1}\rangle , i=0, \ldots, N-1は次のようにして求まります.

\langle X_i X_{i+1}\rangle

\displaystyle{=-\frac{1}{N} \frac{\partial}{\partial (J/\sigma^2)} \ln{Z(\sigma^2)}}

\displaystyle{=-\frac{\sigma^2}{2N}  \frac{\partial}{\partial J}\sum_{j=0}^{N-1}\ln{(1-2J\cos{j\alpha})}}

\displaystyle{=\frac{\sigma^2}{N}\sum_{j=0}^{N-1}\frac{\cos{j\alpha}}{1-2J\cos{j\alpha}}}

\displaystyle{\overset{N \to \infty}{\longrightarrow} \sigma^2 \int_0^1 \frac{\cos{(2\pi q)}}{1-2J\cos{(2\pi q)}}dq}.

1次元Isingモデルの統計量

分散

\displaystyle{\langle X_i^2\rangle = \sigma^2}

共分散

\displaystyle{\langle X_iX_{i+1}\rangle \sim  \sigma^2 \int_0^1 \frac{\cos{(2\pi q)}}{1-2J\cos{(2\pi q)}}dq}

Rによるシミュレーション

この小節では,前節で計算した分散,共分散をR, rstanを利用してシミュレーションしてみたいと思います.

シミュレーションの概要は次の通りです.

1. Ising_model_siml.stanで,1次元Ising模型を定義.

2. Ising_model_worksheet.Rで,Ising_model_siml.stanを呼び出して,適当なN, \sigma^2, Jで乱数発生.

3. 発生させた乱数から,monte carlo原理で分散,共分散をシミュレート.

stanコード

 Ising_model_siml.stan 

こんな感じです.

functions {
real Ising_lpdf(vector X,int N,real sigma, real J) {
//対数尤度
real val = 0;
for(i in 1:N){
if(i<N){
val+= X[i]^2-2*J*X[i]*X[i+1];
} else {
val+= X[i]^2-2*J*X[i]*X[i-N+1];
}
}
return - N*log(sigma) -(1/(2*sigma^2))*val;
}
}

data {
int N;
real sigma;
real J;
}

parameters {
vector[N] X;
}

model {
target+= Ising_lpdf(X|N, sigma, J);
}

Rコード

Ising_model_worksheet.R

 シミュレーションの実行は,次のコードです.ここでは,N= 500, sigma = 1, J = 0.30としてみました.

#Ising_model_worksheet.R
getwd();
library(rstan)
#Simulate N= 500, sigma= 1, J= 0.30
data <- list(N= 500, sigma=1, J= 0)
fit <- stan(file = 'Ising_model_siml.stan', data = data)

 

rstanのサンプリングが終わったら,Rで相関係数を出してみましょう.

# correlation <X_1X_2> ----------------------------------------------------
#get sample data
sample <- extract(fit)
vals <- sample$X[,1] * sample$X[, 2]
#monte carlo
corr <- mean(vals)
corr

コンソールに,corrの値(=\langle X_1X_2 \rangleの近似値)が出力されます.

R.console
R> corr
[1] 0.4249355

ちなみに理論値は,

 \displaystyle{\int_{0}^1\frac{\cos{(2 \pi q)}}{1-2 \cdot 0.30 \cdot \cos{(2\pi q)}} dq= \frac{5}{12}= 0.4166\cdots} 

 ですから,誤差は2%程でしょうか.

 

 

 

 

 

 

pymcを使ってみた

  \require{AMScd,amsmath,amssymb,amsfonts,amscd}
 \def\CC{{\mathbb C}}
 \def\ZZ{{\mathbb Z}}
 \def\PP{{\mathbb P}}
\def\QQ{{\mathbb Q}}
$$\newcommand{\expo}[1]{ \mathrm{exp}(#1)}$$
$$ \newcommand{\coup}[1]{\langle \[ X \] ,#1\rangle} $$

$$\newcommand{\diff}[1]{\frac{d}{dy}#1} $$
$$\newcommand{\pdiff}[1]{\frac{\partial}{\partial y} #1} $$

keywords :ベイズ統計,MCMC(マルコフ連鎖モンテカルロ法),Python

目次

記号:以下,確率変数Xに対してその確率密度関数p(x)で表します.

ベイズ推定の要点

Bayesの公式を理解する

この小節では,Bayesの公式をおさらいすることが目的です.Bayesの公式は一言でいうならば,確率変数の(1)因果関係を利用して,その(2)原因(3)観測結果から推定し直すための公式です.とくに,つぎの3つがキーワードになります:

  1. 因果関係(= 尤度関数)
  2. 原因(= ハイパーパラメーター,事前分布)
  3. 観測結果(= 観測値)

Bayesの公式から得られる,(4)原因の推定結果は,事後分布といいます.言い換えれば,Bayesの公式は(1)-(3)を利用して事後分布を得るために必要な公式といえます.

 

(1)確率変数の因果関係とは

   まず,確率変数の因果関係について詳しく見ていきます.観測のできる量(=観測量)確率変数Xがあるとします.目的は,Xを扱いやすいモデルであてはめてその特徴や傾向を理解することです.扱いやすいモデルの例を挙げておきましょう.

  • (コイン投げ)p\ (0 \leq p \leq 0)をパラメーターとする.Xはパラメーターに依存して確率的以下のようにきまる:

p(X ~|~ p) = \begin{cases} p  \quad ,\text{if}\ X= 1,\\ 1-p \quad,\text{otherwise.}\end{cases}

  • (一次元ガウス分布)\mu, \sigma >0をパラメーターとする.実数値確率変数Xはパラメーターに依存して確率的にきまる.その確率密度pは以下で与えられる:

p(x~|~\mu,\sigma) = \frac{1}{\sqrt{2\pi  \sigma^2}} \exp{\bigl(-\frac{(x-\mu)^2}{2\sigma^2}\bigr)}

以上の2つの例ように,確率変数Xがなんらかのパラメーターに依存している状況において,パラメーター自体をXを引き起こす原因と考えれば良いのです.以下この小節では確率変数Xは,パラメーター\thetaに依存しているとして話を進めます.とくに注意すべきは確率密度f(X~|~\theta)は因果関係を表していることです.つぎのように図で考えるとわかりやすいと思います.

\hspace{100pt}   \text{原因} \ \theta \overset{\text{確率}p(x~|~\theta)}{\longrightarrow} \text{結果}\ x

 

(2)原因

 ベイズ推定の最も大きな特徴は,その原因の捉え方にあります.原因をたくさん考えたいケースを考えます.例えば,(1)のコイン投げの例では原因p(=表の出やすさ)は,0.5(=均等なコイン)だったり,0.8(=よく表が出るコイン)だったり,いくつもの原因があり得ます.もし,あり得る原因が有限個のケースならば各原因ごとに統計検定を実行,比較すればよいですが,そうでないなら違う考え方が必要でしょう.そこでベイズ統計のだす処方箋は,

(HYP)「原因はいろいろありうる.それを確率変数で扱おう」

というものです.この,原因(=パラメーター)をあらわす確率変数をハイパーパラメーターとよびます.もちろん,これだけでは問題を変えただけで原因をどう調べるかという問いに答える必要があります.

 

(3)観測

では,どのようにしたらハイパーパラメーター(=原因)を調べられるでしょうか.一つの方法は,統計検定の考え方を援用する方法です.その考え方とは,

(SAM)「確率変数を詳しく調べたければ,観測を繰り返し行えばいい」

 という考え方です.いま私たちが知りたい確率変数は,ハイパーパラメーターです.そこで(SAM)を考慮すると,次のように考えるのが自然です,

(BAY)「ハイパーパラメーターを詳しく知りたければ,観測を繰り返せばいい」

 これが,Bayes推定の基礎になる考え方です.観測値を考慮(=条件付き確率)したハイパーパラメーターは \theta ~|~ D\quad (D\text{は,観測値})となるのでこの確率分布が知りたい対象です.これは,事後分布(posterior distribution)と呼ばれています.事"後"分布があるなら事"前"分布もあるはずですが,それは単に観測前の\thetaの分布のことです.symbolicに書けば,

\hspace{100pt}   \text{事前分布} \ \theta \overset{\text{観測値}D}{\longrightarrow} \text{事後分布}\ \theta~|~D

 ということです.

 

以上でBayesの公式を述べる準備ができました. 

 Bayesの公式

\bf{Bayes' formula}

[\tex:p(\theta ~|~D) = \frac{p(D~|~\theta)p(\theta)}{p(D)}]

 左辺が,今知りたい事後分布で,それは因果関係p(D~|~\theta),事前分布p(\theta)p(D)(周辺分布という)からわかるというのがBayes' formulaのご利益です.尚,量p(D~|~\theta)尤度関数と呼ばれています.

 

ここで問題が生じます,それは

  1. 尤度関数p(D~|~\theta)は,どう計算すればよいか,
  2. 規格化定数p(D)は,どう計算すればよいか,
  3. 事前分布p(\theta)はどのように選ぶべきか,

という問題です.

1.について.もちろん,考慮する因果関係によるのですが,できるだけ簡単な(=尤度が計算しやすい)モデル(ベータ二項モデル,線形回帰モデル,ロジスティック回帰モデル,etc)を選ぶのが鉄則です.しかし,どうしても簡単なモデルで説明できないならば,個別対応となるでしょう.

 

2.こちらの問題について.p(D)積分\int_{\theta}p(D~|~\theta)p(\theta)d\thetaで求められますが,この積分が求まるケースはむしろ稀です.ですが,パソコンによる数式処理環境が整った現代なら,この積分も簡単なコードで求まります.そのシミュレーション方法の一つが後程でてくるマルコフ連鎖モンテカルロ法(Markov Chain Monte Carlo)です.

 

3.この問題は私の力量不足により答えられません.(勉強したらまた記事に....)

pymcによるシミュレーション

 pymcとは,全節で挙げたBayes推定時の問題のうち,2つめ,つまり規格化定数p(D)の計算を実行するソフトパッケージです.pymcは,python上にimportしてつかうのでpythonを利用したデータ分析時に威力を発揮します.今回はpymcの動かし方をソースコードとともに紹介します.尚,今回のソースコードは次の書籍を参考にしています.

 

Pythonで体験するベイズ推論 PyMCによるMCMC入門

Pythonで体験するベイズ推論 PyMCによるMCMC入門

 

 

例(ロジスティック回帰分析)

つぎの問題を通して,pymcの使い方を紹介してみます.

問題:Aさんは毎日,9:00の始業時間に合わせ,8時台に家をでる.sample.csvは,Aさんの家の出発時間(8:00からの経過時間(分))departure_time,その日の天気weather(雨:1,晴:0),そして遅刻したかどうかresult(遅刻:1 ,セーフ:0)を連ねた10件のデータからなる.Aさんの遅刻する確率を出発時間からもとめよ.

今回は簡単のため,天気を考慮せずに考えてみます.

使いたいライブラリをimportします.


 #libraries
import numpy as np import pandas as pd from scipy import stats #library for plot from matplotlib import pylab as plt

csvファイルを読み込むために,pandasライブラリのread_csvを利用します.データ数はたかが知れているので,表示させて確認もしておきます.


    sample = pd.read_csv("sample.csv")#csv読み込み
sample#データの中身

つぎに散布図を書きます.


#散布図
plt.scatter(sample['departure_time'],
    sample['result'])
plt.ylabel("1:delay 0:in time")
plt.xlabel("departure_time(min)")
plt.title("departure time versus delay")
plt.show()

 

f:id:kimigayoseishou:20180522192600j:plain

出発時間と遅刻の関係

今回はデータが簡単なため,散布図を書いただけでも「8:20-30分あたりが遅刻するかどうかの分かれ目っぽい」と推測できます.

logistic回帰?

#事前分布 
beta0 = pm.Uniform("beta0", -100, 100, value = 1 / 200)#(-100,100)上の一様事前分布 
beta1 = pm.Uniform("beta1", -100, 100, value = 1 / 200)

 


#logistic回帰
#逆logit関数
@pm.deterministic
def invlogit(d = departure, beta0 = beta0, beta1 = beta1): 
    return 1.0 / (1.0 + np.exp(beta0 + beta1 * d)) 
#モデル式.ここから,尤度が計算される.
obs = pm.Bernoulli("bernoulli_obs", invlogit, value = result, observed = True)

#MCMCによるBayes推定
#Model classを生成:Model(モデル式, 事前分布)
model = pm.Model([obs, beta0, beta1])

#サンプリング
mcmc = pm.MCMC(model)
mcmc.sample(40000, 10000)    

    #結果を参照
from pymc import Matplot 
Matplot.plot(beta0)#事後分布の平均は14らへん
Matplot.plot(beta1)#事後分布の平均は-0.5らへん
#png形式でディレクトリに保存される

f:id:kimigayoseishou:20180522194227j:plain

beta0の推定結果

f:id:kimigayoseishou:20180522194235j:plain

beta1の推定結果

では得られた推定値と,元のデータを比較してみましょう.8:20分台後半からが,勝負の分かれ目であることがわかります.

f:id:kimigayoseishou:20180522211254j:plain

推定値と元のデータの比較

Hirzebruch's Riemann-Roch theorem とBorisov's formula(2)

<script type="text/javascript" src="https://c328740.ssl.cf1.rackcdn.com/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
 \require{AMScd,amsmath,amssymb,amsfonts,amscd}
 \def\CC{{\mathbb C}}
 \def\ZZ{{\mathbb Z}}
 \def\PP{{\mathbb P}}
\def\QQ{{\mathbb Q}}
$$\newcommand{\expo}[1]{ \mathrm{exp}(#1)}$$
$$ \newcommand{\coup}[1]{\langle \[ X \] ,#1\rangle} $$
$$\newcommand{\diff}[1]{\frac{d}{dy}#1} $$
$$\newcommand{\pdiff}[1]{\frac{\partial}{\partial y} #1} $$

今回記事は,前回記事で紹介した式

sonnamonyaro.hatenablog.com

{\sf Theorem}[ {\rm Borisov\ proposition2.2} ]
Xn次元コンパクトK\overset{..}{a}hler多様体とする.\Omega_XX上の正則微分形式の層をあらわす.整数p,qに対してHodge数をh^{p,q}(X):={\rm dim}_\CC H^q(X,\Omega^p_X)とするとき,等式
$$ \sum_{p,q \in \ZZ}(-1)^{p+q}\Bigl(q-\frac{n}{2} \Bigr)^2h^{p,q}(X) $$
$$=\frac{1}{12}n \chi^{\rm top}(X)+\frac{1}{6}\int_Xc_1(X)\cup c_{n-1}(X),\quad (1)$$
が成り立つ,ここでc_j(X)は接束TXの第j\text{-}chern類をあらわす.

の証明の後半となります.

 

 


{\sf Lemma}
Xn次元compact K\overset{..}{a}hler多様体とする.このとき,つぎの等式が成り立つ:
\sum_{p,q\in \ZZ}(-1)^{p+q}\Bigl(q-\frac{n}{2}\Bigr)^2h^{p,q}(X)=\sum_{p,q\in \ZZ}q^2(-1)^{p+q}h^{p,q}(X)-\frac{n^2}{4}\chi^{{\rm top}}(X).
{\sf proof}
前回記事で示した等式
$$\sum_{p,q\in \ZZ}(-1)^{p+q}q\cdot h^{p,q}(X)=\frac{n}{2}\cdot \chi^{\rm top}(X),\quad \text{(1)}$$
を使えばよいです.\square

 

つぎに求めるべきは\sum_{p,q\in \ZZ}q^2(-1)^{p+q}h^{p,q}(X)ですが,これは\chi \text{-}y種数を用いれば-\diff{\Bigl(y\diff{\chi(X;y)} \Bigr)}\Bigr|_{y=-1}
とかけます.これは,書き直すと\frac{d^2}{d^2y}\chi(X;y)\Bigr|_{y=-1}+\Bigl\{\frac{n}{2}-\frac{n^2}{4}\Bigr\}\chi^{{\rm top}}(X)です.そこで,\frac{d^2}{d^2y}\chi(X;y)\Bigr|_{y=-1}をHirzebruchのRiemann-Rochに持ち込みます.
{\sf Lemma}
Xn次元compact K\overset{..}{a}hler多様体とする.このとき,つぎの等式が成り立つ:
\frac{d^2}{d^2y}\chi(X;y)\Bigr|_{y=-1}=\frac{1}{6}\coup{c_1\cup c_{n-1}}+\frac{1}{2}\binom{n}{2}\coup{c_n}\quad \text{(2)}.
{\sf proof}
$$\frac{d^2}{d^2y}{\rm td}_y(X)$$
$$=\sum^n_{i=1}Q(\gamma_1)\cdots \underbrace{\frac{d^2}{d^2y}Q(\gamma_i)}_{\text{ここだけ}} \cdots Q(\gamma_n)\Bigl|_{y=-1}$$
$$+\sum_{i\neq j}Q(\gamma_i)\cdots \underbrace{\pdiff{Q(\gamma_i)}}_{\text{ここだけ}} \cdots \underbrace{\pdiff{Q(\gamma_j)}}_{\text{ここだけ}}\cdots Q(\gamma_n)\Bigr|_{y=-1}$$
$$=\sum^n_{i=1}( 1+\gamma_1)\cdots \underbrace{\frac{1}{6}\gamma_i{}^2}_{\text{ここだけ}} \cdots ( 1+\gamma_n)$$
$$+\sum_{i\neq j}( 1+\gamma_i)\cdots \underbrace{\frac{-1}{2}\gamma_i}_{\text{ここだけ}} \cdots \underbrace{\frac{-1}{2}\gamma_j}_{\text{ここだけ}}\cdots ( 1+\gamma_n)$$
である.よって,\frac{d^2}{d^2y}{\rm td}_y(X)H^{2n}(X;\QQ)の項は,
$$\frac{1}{6}( c_1(X)\cup c_{n-1}(X)-nc_n(X))+\frac{1}{2}\binom{n}{2}c_n(X)$$
HirzebruchのRiemann-Rochより,
$$\frac{d^2}{d^2y}\chi(X;y)\Bigr|_{y=-1}$$
$$=\Bigl\langle \[ X \] ,\frac{d^2}{d^2y}{\rm td}_y(X) \Bigr\rangle $$
なので,主張はしたがう.\square

Ligober=Wood,Borisovの式の証明


(1)(2)より,
$$\sum_{p,q\in \ZZ}(-1)^{p+q}\Bigl(q-\frac{n}{2}\Bigr)^2h^{p,q}(X)$$
$$=\sum_{p,q\in \ZZ}q^2(-1)^{p+q}h^{p,q}(X)-\frac{n^2}{4}\chi^{{\rm top}}(X),\quad \text{(1)より}$$
$$=\frac{1}{6}\coup{c_1\cup c_{n-1}}+\frac{1}{2}\binom{n}{2}\coup{c_n}-\frac{n^2}{4}\chi^{{\rm top}}(X)\quad \text{(2)より}$$
$$=\frac{1}{6}\int_Xc_1(X)\cup c_{n-1}(X)$$
$$+\Bigl\{ \frac{1}{2}\binom{n}{2}-\frac{n}{6}+\frac{n}{2}-\frac{n^2}{4}\Bigr\}\chi^{{\rm top}}(X)$$
$$=\frac{1}{12}n \chi^{\rm top}(X)+\frac{1}{6}\int_Xc_1(X)\cup c_{n-1}(X).\quad \square$$

Hirzebruch's Riemann-Roch theorem とBorisov's formula(1)

<script type="text/javascript" src="https://c328740.ssl.cf1.rackcdn.com/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
 \require{AMScd,amsmath,amssymb,amsfonts,amscd}
 \def\CC{{\mathbb C}}
 \def\ZZ{{\mathbb Z}}
 \def\PP{{\mathbb P}}
\def\QQ{{\mathbb Q}}
$$\newcommand{\expo}[1]{ \mathrm{exp}(#1)}$$
$$ \newcommand{\coup}[1]{\langle \[ X \] ,#1\rangle} $$

$$\newcommand{\diff}[1]{\frac{d}{dy}#1} $$
$$\newcommand{\pdiff}[1]{\frac{\partial}{\partial y} #1} $$
今回から2回に分けて次の公式の導出を紹介したいと思います:
{\sf Theorem}[ {\rm Borisov\ proposition2.2} ]
Xn次元コンパクトK\overset{..}{a}hler多様体とする.\Omega_XX上の正則微分形式の層をあらわす.整数p,qに対してHodge数をh^{p,q}(X):={\rm dim}_\CC H^q(X,\Omega^p_X)とするとき,等式
$$ \sum_{p,q \in \ZZ}(-1)^{p+q}\Bigl(q-\frac{n}{2} \Bigr)^2h^{p,q}(X) $$
$$=\frac{1}{12}n \chi^{\rm top}(X)+\frac{1}{6}\int_Xc_1(X)\cup c_{n-1}(X),\quad (1)$$
が成り立つ,ここでc_j(X)は接束TXの第j\text{-}chern類をあらわす.

\chi_y種数

{\sf definition}[{\rm  Hir}15.5.]
Xn次元compact複素多様体とする.このとき不定元yに関する多項式\chi(X;y)を次で定める:
$$\chi(X;y)=\sum_{p\in \ZZ}\chi^{\rm hol}(\Omega_X^p)y^p.$$
このとき\chi(X;y) \chi_y\text{-}{\bf \text{種数}}という.
正則Euler数の定義から
$$\chi(X;y)$$
$$=\sum_{p\in \ZZ}\chi^{\rm hol}(\Omega_X^p)y^p$$
$$=\sum_{p\in \ZZ}\sum_{q \in \ZZ}(-1)^q{\rm dim}_\CC H^q(X,\Omega_X^p)y^p$$
$$=\sum_{p,q\in \ZZ}(-1)^qh^{p,q}y^p,\quad \text{(2)}$$とかくことができる.<\p>

{\sf Proposition}[{\rm Hir}15.8.]
Xがcompact K\overset{..}{a}hler多様体であるとする.\chi(X;y)の引数として-1を選ぶと,その値はXの位相的Euler数になる:
$$\chi(X;-1)=\chi^{\rm top}(X),\quad (3).$$

Xがcompact K\overset{..}{a}hler多様体ならば,Hodge分解により等式
$${\rm dim}_{\CC}H_{\rm dR}^j(X;\CC)=\sum_{p+q=j}h^{p,q}(X)$$
,(H_{\rm dR}^j(X;\CC)Xの複素係数de Rham cohomology群)がなりたつ.
これ用いて次のように変形すればよい:
$$\chi(X;-1)$$
$$=\sum_{p,q\in \ZZ}(-1)^{p+q}h^{p,q}(X)$$
$$=\sum_{j\in \ZZ} (-1)^j\sum_{p+q=j} h^{p,q}(X)$$
$$=\sum_{j\in \ZZ}(-1)^j{\rm dim}_{\CC}H_{\rm dR}^j(X;\CC)$$
$$=\chi^{\rm top}(X) .\quad \square$$

{\sf Example}
複素射影空間\CC \PP^nについて,その\chi_y\text{-}種数はつぎのようにもとまる.
\CC \PP^nのHodge数はつぎのようになる
$$h^{p,q}(\CC\PP^n)=\begin{cases}1 &(0 \leq p=q\leq n)\\0 &(\text{otherwise})\end{cases}. $$
よって(2)より
$$\chi(\CC\PP^n;y)=\sum_{q=0}^n(-1)^qy^q=\frac{1-(-y)^{n+1}}{1+y}.$$

一般Todd種数


Xを向き付けられた実2n次元微分可能多様体に対してH^i(X;\QQ)i
係数コホモロジー群とする.cup積
$$ \cup:(a,b)\mapsto a\cup b\in H^{i+j}(X;\QQ),\quad a\in H^i(X;\QQ),b\in H^j(X;\QQ) $$
を考えるとH^{\text{even}}(X;\QQ)=\underset{i\in 2\ZZ}{\bigoplus }H^{i}(X;\QQ)\ZZ次数つき可換\QQ\text{-}代数になる.Xの向きから定まる基本類[X]\in H_{2n}(X;\QQ)との a\in H^{2n}(X;\QQ)とのカップリングを\coup{a}\in \QQとかく.
{\sf definition}
Xn次元compact複素多様体とし,複素構造から誘導された向きをあたえる.Xの正則接束TXの第ichern類(i=0,\ldots,n)c_i\in H^{2i}(X;\QQ)であらわす.chern root \gamma_i,\ldots,\gamma_n\in H^2(X;\QQ)を関係式
$$(1+\gamma_1)\cdots(1+\gamma_n)=1+c_1+\cdots+c_n$$
で定めておく.
\QQ [y]を係数環にもつ冪級数R:=(\QQ [ y] )[ [ w] ] の元
$$Q(y;w):=\frac{w(1+y)}{1-\expo{-w(1+y)}}-yw \in R$$
を用いてX{\bf \text{一般化} Todd\text{種数}}多項式
$$\text{Td}(X;y):=\coup{\text{td}_y(X)}\in \QQ[y],$$

$$\text{ただし},\quad \text{td}_y(X):=Q(\gamma_1)\cdots Q(\gamma_n)$$
で定義する.

Hirzebruch genera of R-R theorem

 

定理(1)の証明に使う次の定理を紹介する
{\sf theorem}[{\rm  Hir}21.3.]
Xを射影的複素多様体とする.このとき,
$${\rm Td}(X;y)=\chi(X;y).$$
この定理の証明はここではしないです.しないです.

y=-1のとき,(3)より\chi^{\rm top}(X)={\rm Td}(X,-1)=\coup{c_{\text{dim}_\CC X}}となる.これは,Xに対するGauss-Bonne-Chernの公式である.

定理(1)の証明(前編)


さて定理(1)の証明に入る.
{\sf Lemma}
Xn次元compact K\overset{..}{a}hler多様体とする.このとき,つぎの等式が成り立つ:
$$\sum_{p,q\in \ZZ}(-1)^{p+q}q\cdot h^{p,q}(X)=\frac{n}{2}\cdot \chi^{\rm top}(X).$$
{\sf proof}
Q(y;w)y=-1で展開すると,
$$Q(y;w)$$
$$=\Bigl\{ 1+\frac{(1+y)}{\expo{(1+y)w}-1}\Bigr\}\cdot w $$
$$=\Bigl\{1+ \frac{1}{w}\Bigl( 1-\frac{(1+y)w}{2}+\frac{(1+y)^2w^2}{12}+O((1+y)^3)\Bigr)\Bigr\}\cdot w $$
$$=1+w-\frac{(1+y)w}{2}+\frac{(1+y)^2w^2}{12}+O((1+y)^3)$$
である.よって,
$$\diff{\text{td}_y}(X)\Bigr|_{y=-1}$$
$$=\sum_{j=1}^nQ(-1;\gamma_1)\cdots \underbrace{\pdiff{Q(y;\gamma_j)}}_{\text{ここだけ}}\Bigr|_{y=-1}\cdots Q(-1;\gamma_n)$$
$$=-\frac{1}{2}\sum_{j=1}^n(1+\gamma_1)\cdots \underbrace{\gamma_j}_{\text{ここだけ}} \cdots (1+\gamma_n)$$
となる.せやから
$$\diff{\chi(X;y)}\Bigr|_{y=-1}$$

$$=\Bigl\langle \[ X \] {}_,\diff{ {\rm td}_y }(X) \Bigr| _{y=-1} \Bigr\rangle $$

$$=-\frac{1}{2}\Bigl\langle \[ X \] {}_,\sum_{j=1}^n(1+\gamma_1)\cdots \underbrace{\gamma_j}_{\text{ここだけ}}\cdots (1+\gamma_n)\Bigr\rangle$$

$$=-\frac{1}{2}n\cdot \coup{\gamma_1\cdots\gamma_n}$$
$$=-\frac{n}{2}\cdot \chi^{\text{top}}(X)$$
がわかる.
$$\diff{\chi(X;-1)}\Bigr|_{y=-1}=-\sum_{p,q\in \ZZ}(-1)^{p+q}q\cdot h^{p,q}(X)$$
なのでえた.\square
今日のところはここまでとします(・\omega・)/

参考文献


[{\rm Hir}]:\text{F.Hirzebruch},{\it Topological\  Methods \ in \ Algebraic\  Geometry },\text{Third Edition}

 \text{,Graduate Texts in Mathematics ,No.52,New York ,Springer Verlag }(1978).

[{\rm Bor}]:\text{L. Borisov}{\it  On \ Betti\  numbers \ and \ Chern\  classes }

{\it of \ varieties\  with \ trivial \ odd \ cohomology \ groups.\ }\text{arXiv: alg-geom/9703023.}

 

 

ブラウン運動の最尤推定

 \require{AMScd,amsmath,amssymb,amsfonts,amscd}

 \require{color}

keywords :時系列解析 , 最尤法 , Kullback-Leibuler情報量 , ブラウン運動

今回のテーマは、前回に引き続き時系列解析の話題です。

前回の記事ではワイトノイズ、ブラウン運動、AR(1)過程について 

その平均、自己共分散がどのように変化するかについて書きました。 

sonnamonyaro.hatenablog.com

 

今回の記事では、これらのモデルを実際に応用する局面で最も基本的な問題

 

問題:与えられた時系列データを最も適切に再現するように時系列モデルを決定せよ 

 

について、最尤法という手法により上の問題を扱いたいとおもいます

 

目次

 

最尤法

最尤法は統計学で広く使われる考え方で、パラメーターを持った確率変数

に従うサンプルから、適当なパラメーターを決定する方法 と言えます。

例えば、ここに表が出る確率がp\ (0<p<1)のコインがあったとして、

pを求めろと上司に言われたとします。そうしたら、あなたはおそ

らく(上司に逆らわずに)10回、100回、・・・とコインをひたすら投げ、

表の出た割合をの値として報告するでしょう。

実はそれが、最尤法なのです。

それでは、最尤法の説明にはいります

尤度関数

\theta \in \mathbb{R}をパラメータに持つ確率変数

X_\thetaを考えます。またはX_\theta連続な確率密度f(x,\theta)

を持つと仮定します。先の例でいうと、pがパラメータ\theta

コインを1回投げた時、表が出たら1裏が出たら0とするカウントがX_\theta 

になります。

{\bf definition}

尤度関数(Likelihood function)とは、上のように用意されたパラメーター

(\theta)付き確率変数と、観測されたN個のデータ\mathsf{D}=(x_1,x_2,\cdots,x_N)に対して定まる、\thetaの関数でつぎで定義される

{\bf (LF)}\ \mathcal{L}(\theta|\mathsf{ D})=\mathsf{D}\text{が観測される確率}

もし、個のデータが一回ずつ採取され、それらが独立と仮定する場合には

(先ほどのコインの例)積の法則により

{\bf (LF')}\ \mathcal{L}(\theta|\mathsf{ D})=\prod_{i=1}^{N}f(x_i,\theta)

 この式が、よく皆さんが目にする式ではないかと思います。なお、section2で

は(LF')の式ではなく(LF)の式を使うので注意してください

最尤法は次のような考え方です

 尤度関数を最大化するパラメーター\theta^\astは"適切"な選び方である

ここで導入された、尤度関数を最大化するパラメーターを最尤推定量 

({\bf Maximum \ Likelihood \ Estimator/\  M.L.E.})といいます。覚えておきましょう。

また、確率密度が微分可能ならば、最尤推定量での尤度関数の微分は消えることから

次の尤度方程式が成り立つことがわかります

 \frac{d\ \text{ln}\ \mathcal{L}(\theta^\ast|\mathsf{ D})}{d\ \theta}=0
 

 

 

ここで、適切の意味は次のような事実をさします

{\bf Fact}

最尤推定量は、サンプルを死ぬほど採りまくった場合に、正しい値を教えてくれる

ことが保証されている(強一致性,漸近正規性,...)

これらの事実は統計学のまともな本に譲るとして、うえの事実を納得するための

考え方を紹介するにとどめます 。

Kullback-Leibuler情報量

Kullback-Leibuler 情報量は、2つの確率分布が等しいか否かをはかる量です。

今回の記事では、Kullback-Leibuler情報量と最尤法のかかわりをみてみます。

{\bf definition}
\mathsf{A},\mathsf{B}を確率分布、それらの確率密度f(x),g(x)\ (x \in \mathbb{R})
は連続とする。このときKullback-Leibuler 情報量は次の式で定義される
I(\mathsf{A}||\mathsf{B})=\int_{\mathbb{R}}dx\ f(x){\rm ln} \Bigl(\frac{f(x)}{g(x)}\Bigr)

 

次のことが知られています。

重要性質
\mathsf{A},\mathsf{B}が確率分布であるならば、
I(\mathsf{A}||\mathsf{B})\geq 0\ ,\text{等号成立は}\ \mathsf{A}=\mathsf{B}\ \text{のときに限る}

 

うえの性質より、2つの確率分布が等しいかどうかをみるには、それらのKullback-Leibuler情報量が0かどうかに注目すればよいのです。

 

さてこのことを次のパラメーター推定問題

{\bf Problem}

確率分布{\sf A}(確率密度f(x))と、パラメーター付き確率分布
{\sf A}(\theta)(確率分布f(x,\theta))が与えられたとき、
\mathsf{A}_\thetaの分布が\mathsf{A}に最も分布が近くなる\thetaをあたえよ

に応用してみます。
0\leq I(\mathsf{A}||\mathsf{A}_\theta)
=\int_{\mathbb{R}}dx\ f(x){\rm ln}\ f(x)-\int_{\mathbb{R}}dx\ f(x){\rm ln}\ f(x,\theta)
なので分布\mathsf{A}_\theta\mathsf{A}に近づけるためには、上の式の第二項
\int_{\mathbb{R}}dx\ f(x){\rm ln}\ f(x,\theta)
を最大にすればよいことがわかります。

 

すこしずつ最尤法がにちかずいてきました。もうすこし頑張りましょう。
結局上の積分直接計算するには、f(x)を知らなければなりません。しかし、\mathsf{A}から採ったサンプル\mathsf{D}=\{x_1,x_2,\cdots,x_N\}]があれば、

そこからモンテカルロ法を利用して間接的積分を計算できます:
\int_{\mathbb{R}}dx\ f(x){\rm ln}\ f(x,\theta)\sim \frac{1}{N}\sum_{i=1}^N{\rm ln}\ f(x_i,\theta)

これって対数尤度ですね。

つまり、最尤法の原理「尤度関数を最大化せよ」というのは、真の分布と
モデル分布のKullback-Leibuler情報量が小さいことが期待されるという点で、適当なパラメーターの選び方であると言えるわけです。

ブラウン運動最尤推定

実際に最尤法のコツをつかむために一次元ブラウン運動最尤推定について
解説します。

今回考えるブラウン運動は、2つのパラメーター\alpha,\sigma^2
を持つ時系列モデルX_t\ (t=0,1,2,\cdots)で次の式で与えられるものです

 ({\bf BM1})X_0=0

 ({\bf BM2})X_{t}=\alpha +X_{t-1}+w_t
 ({\bf BM3})w_t \underset{\text{i.i.d.}}{\sim} N(0,\sigma^2)
 
さて、このデータに従う時系列データが得られたとします。

簡単のため、時系列データは時刻1から$T \geq 1$まで連続的に取ったものとし、

それを\mathsf{D}=\{x_0=0, x_1,\cdots,x_T\}
とします。 

尤度関数の計算

今回は、パラメーターは\alpha,\sigma^2なので尤度関数は2変数に
なります。それを、 \mathcal{L}( \alpha,\sigma^2|\mathsf{D})
とおきます。

\mathcal{L}( \alpha,\sigma^2|\mathsf{D})
=\mathbb{P}[ X_1=x_1\  \text{and}\  X_2=x_2 \ \text{and}\  \cdots \ \text{and}\  X_T=x_T]
=\mathbb{P}[ X_1=x_1 \ \text{and}\  X_2-X_1=x_2-x_1 \ \text{and}\  \cdots \ \text{and} \ X_T-X_{T-1}=x_T-x_{T-1}]
ここで式({\bf BM2})よりX_i-X_{i-1}=w_i+\alphaで、これらが({\bf BM3})

より独立であったことから、尤度関数が
\mathcal{L}( \alpha,\sigma^2|\mathsf{D})
\mathbb{P}[ w_1+ \alpha =x_1\  \text{and}\  w_2 +\alpha=x_2-x_1 \ \text{and}\  \cdots \ \text{and}\  w_T +\alpha=x_T-x_{T-1}]
=\prod_{i=1}^T \Bigl(\frac{1}{\sqrt{2 \pi \sigma^2}}\text{exp}\bigl( \frac{(x_i-x_{i-1}-\alpha)^2}{2\sigma^2}\bigr) \Bigr)
ともとまりましたね。目がちかちかするので、対数を取りましょう
\text{ln}\ \mathcal{L}( \alpha,\sigma^2|\mathsf{D})
=-\frac{T}{2}\text{ln}(2\pi \sigma^2)-\frac{1}{2\sigma^2}\sum_{i=1}^T(\alpha-(x_i-x_{i-1}))^2

尤度方程式
ここまでこれば、あとは尤度方程式を解けばオッケーです
0=\frac{\partial \ \text{ln}\ \mathcal{L} }{\partial \alpha}=-\frac{1}{ \sigma^2}\sum_{i=1}^T(\alpha-(x_i-x_{i-1}))
0=\frac{\partial \ \text{ln}\ \mathcal{L} }{\partial \sigma^2}=-\frac{T}{2\sigma^2}+\frac{1}{2(\sigma^2)^2}\sum_{i=1}^T(\alpha-(x_i-x_{i-1}))^2

これを連立して解けば次のように最尤推定量がもとまります
\alpha^\ast=\frac{x_T}{T}
\sigma^{\ast 2}=\frac{1}{T}\sum_{i=1}^T(x_i-x_{i-1})^2- \frac{x_T^2}{T^2}

とめでたく求まりました。

式の形を見れば気づくかもしれませんが、この最尤推定は暗算でやる方法

があるので、考えてみるとよいでしょう。

 

時系列解析入門

 \require{AMScd,amsmath,amssymb,amsfonts,amscd}

 \require{color}

key words:時系列解析 ,  ブラウン運動( Brownian Motion) ,AR過程

今回のテーマは時系列解析です。

時系列解析の目的は、確率的に変化する情報(確率変数)を時刻ごとに

追って分析することです。数学、物理学、経済学などの各分野で実用の

機会も多い題材です。今回はその中でも最も単純な例である、

について、基本的なことをまとめます。予備知識は仮定しませんので

ご心配なく。

目次

時系列解析

統計学のおさらい

このセクションでは、時系列解析で広く使われている統計学のツールを導入します。

まずは時系列解析の対象である、時刻ごとに変化し、しかも確率的である情報 というものをはっきりさせる必要があります。そこで次の定義を採用します。

definition 

 t= 0,\pm1,\pm2,\cdots を時刻として, 各 t に対して、割り当てられた 

実数値確率変数  X_t (t=0,\pm1,\pm2,\cdots) があるとき、これを時系列データ 

という。

 

X_0 が例えば現在の情報であるなら、 X_{ \textcolor{red}{1}}は1秒の情報、 X_{ \textcolor{red}{-1}}

は1秒の情報と理解すればよいです。時系列解析によって、現在の情報、及び過去

の情報から、未来の情報を予測(=確率的に知る)できるわけです。

さて、確率変数Xが与えられたとき、その平均、分散を考えることができます。

それらをそれぞれ

 \begin{equation} {\mathbb E}[X] ,  {\mathbb V} [X] \end{equation}

と書くことにします。これらは1つの確率変数に対して定まる量ですが、時系列解析

では無限個の確率変数が登場しているため、2つの確率変数 (例えば、現在X_0と1秒後X_1に対して定まる量も重要になります。特に、共分散は、時系列解析で特に重要です。共分散は、{\mathbb R^2}上の確率変数  (X,Y) に対して次で

定義されます

  \text{Cov}[X,Y ]

={\mathbb E}[(X-{\mathbb E}[X])(Y-{\mathbb E}[Y])].

また共分散をー1から1の間に入るよう正規化した相関係数は次で与えられます

 \rho(X,Y)=\frac{ \text{Cov} [X,Y]}{ \sqrt{{\mathbb V} [X] {\mathbb V} [Y] }}.

 この数字が1に近いとき、[ ex:X,Y]は協調関係にあり(つまり、X,Yの大小が同時に起きる)、また-1に近いときトレード・オフの関係(つまり、X,Yの大小が相反的に起こる)になります。

今回は主に時系列データ[tex; X_t]および時刻t,t'に対して(X_t,X_{t'})の共分散をしばしば考えます。これを自己共分散autocovarianceといいます。

 

このセクションの最後に正規分布の再生性をおさらいします

Theorem

確率変数の組(X_1,X_2)があり、各X_i (i=1,2)が平均 \mu_i、分散 \sigma^2_i正規分布N(\mu_1,\sigma_1^2)に従っているとする。このとき、確率変数X_1\pm X_2正規分布N(mu_1 \pm \mu_2,\sigma^2_1+\sigma^2_2)

 に従う。

 

定常性

 

さてここでは、時系列データに対して定常性(stationarity)という概念を定めます。

定常というのは、より詳しくは時刻に関して一定であるという意味だと考えてください。それでは定義。

definition

時系列データX_t (t=0,\pm 1,\pm 2,\cdots)定常であるとは次の3条件

  • (SA){\mathbb E}[X_t]tに依らない
  • (SV){\mathbb V}[X_t]tに依らない
  • (SAC) \text{Cov}[X_t,X_{t'}]は差の絶対値|t-t'|のみに依存する

が満たされるときをいう。☐

以降この記事の中では、条件(SA)/(SV)/(SAC)をそれぞれ定常平均を持つ /定常分散を持つ/定常自己相関を持つ、と呼びますが一般的な呼称でないこと注意しておきます。
definition

時系列データX_t (t=0,\pm 1,\pm 2,\cdots)定常増分(stationary increments)をもつとは、時刻t>t'に対してX_t-X_{t'}X_{t-t'}と同分布であるときをいう。

では次セクションでこれらの性質を例と照らし合わせながら見ていきます。

 

時系列解析の実例

さていよいよ時系列解析で実際に使われる例について話します。今回紹介するのは、ホワイトノイズ、ブラウン運動、AR(1)過程の3つですが、これらそれぞれ、高校で登場した2項漸化式

  • x_{n+1}=c (c: \text{定数}) 定数
  • x_{n+1}-x_n=c (c: \text{定数}) 線形
  • x_{n+1}-\alpha x_n=c (\alpha,c: \text{定数}) 指数型

の時系列解析における類似物と考えられます。 

 

ホワイトノイズ 

definition2.1.

 ホワイトノイズ白色雑音)は、次で与えられる時系列データX_t (t=0,\pm1,\pm2,\cdots)のことである

  • (WN1)X_t =w_t ,w_t:正規分布 N(0,\sigma^2_{ \text{noise}}) に従う確率変数。
  • (WN2)確率変数の組(X_t,X_{t'}) (t, \neq t')は独立。

ここで、\sigma^2_ \text{noise}は>0は定数である。□

{\mathbb E}[X_t]=0 ,{\mathbb E}[X_t]=\sigma^2_{ \text{noise}} となるので、

ホワイトノイズは、定常平均と定常分散を持つことになります・・・①。 

また、時刻t,s (t>s)についてX_t-X_s の分布は(WN1,2)と正規分布

再生性より常に平均0分散2\sigma^2_{ \text{noise}}正規分布であるため、定常増分を持たないことがわかります。

次に、ホワイトノイズの自己共分散\phi(t,t’)はというと、次のようになります

 \phi(t,t')=\begin{cases} \sigma^2_{ \text{noise}} ( \text{if } t=t') \\0 ( \text{if } t \neq t')\end{cases}

 これは、t=t'のときはWN1から、t \neq t'のときはWN2からわかります・・・②。

①②よりホワイトノイズは定常過程であることがわかります

ブラウン運動

definition2.2.

ブラウン運動とは次で定まる時系列データX_t (t=0,1,2,cdots)のことである。

 

  • (BM1) X_0=0
  • (BM2) w_t:=X_{t}-X_{t-1} (t=1,2,\cdots) は常に正規分布 N(0,\sigma^2_ \text{noise})に従う
  • (BM3) すべてのj=1,2,\cdots及び、非負整数t_1 <t_2 <\cdots <t_jについて\begin{equation}X_{t_1},X_{t_2}-X_{t_1} , \cdots,X_{t_j}-X_{t_{j-1}}\end{equation}は独立

ここで、{\sigma^2}_{ \text{noise}} >0は定数である。□

条件(BM3)は、独立増分性と呼ばれている条件です。

さて、まず時刻tにおける分布について調べましょう。

proposition2.3.

{X_t}_{t=0,1,2,\cdots}がdefinition2.1.で与えられるブラウン運動であるとき、

X_tは平均0,分散t\sigma_{ \text{noise}}^2正規分布に従う。

このことを見るには、次のようにします

X_t=X_0+(X_1-X_0)+(X_2-X_1)+\cdots+(X_t-X_{t-1})

と書いたときX_j-X_{j-1}=w_j (j=1,2,\cdots,t)正規分布N(0,\sigma^2_ \text{noise})に従い(BM2)ます。さらにこれらは独立(BM3)なので、正規分布の再生性よりX_tは平均t \times 0,t \times \sigma^2_{ \text{noise}}正規分布であることが言えます。

proposition2.3.より、ブラウン運動は平均定常性は持つが、分散定常性は持たないことがわかりました。うえと同様に、時刻t,s (t>s) のとき,

X_t-X_s=(X_{s+1}-X_s)+\cdots+(X_t-X_{t-1})

に着目してX_t-X_sは平均0,分散(t-s)\sigma^2_{ \text{noise}}正規分布であることがわかります。これは、X_{t-s}と同分布であるためブラウン運動定常増分を持つことになります。

さて次にブラウン運動の自己共分散\phi(t,t') t>t'を求めるには次のようにします。
(BM3)よりX_{t'}X_t-X_{t'}は独立なので
 \text{Cov}(X_t-X_{t'},X_{t'})は0になることから、
\phi(t,t')

= \text{Cov}(X_{t},X_{t'})

= \text{Cov}(X_{t}-X_{t'},X_{t'})+ ext{Cov}(X_{t'},X_{t'})

={\mathbb V}[ X_{t'}]

=t' \sigma_{ \text{noise}}^2.

けっか\phi(t,t')= \text{min}{t,t'}\sigma_{ \text{noise}}^2
となりますが、これは|t-t'|の関数ではないため、ブラウン運動は自己共分散に
関して定常性を持ちません。

AR(1)過程

definition 

時系列データX_t (t=0,\pm1 ,\pm2,\cdots)がparameter \alpha,\beta,\sigma_{ \text{noise}}^2をもつAR(1)過程(AR(1) process)に従うとは、次の2条件

  • (AR1)\begin{equation}w_t:=X_t-\alpha X_{t-1}\end{equation}は、平均0分散\sigma_{ \text{noise}}^2正規分布に従う
  • (AR2)t \neq sならば,w_t,w_sは独立。

が成り立つ時をいう。

さてここで例えば、parameterを\alpha=1,\beta=0と選ぶと、AR方程式はブラウン運動を定める式(BM2)に一致してしまいます。するとそのparameterで定まるAR(1)過程も、ブラウン運動が分散定常性を持たないのと同様に、分散定常性を失ってしまいます。じつは\textcolor{red}{\alpha=\pm1}を境にAR(1)過程は定常性に関して異なる性質を持つのを述べたのが、つぎの命題になります

proposition 

AR(1)過程が分散定常性を持つならば、|\alpha|<1

 

(証明)
{\mathbb{V}}[X_t]=\sigma^2>0 と仮定する。(AR1)を繰り返し用いて、[h=1,2,\cdots]のとき

X_t=\alpha^h X_{t-h}+\beta \sum_{j=0}^{h-1}\alpha^j+\sum_{j=0}^{h-1}\alpha^j w_{t-h}・・・①

となる。この等式の右辺の分散を計算すると、\alpha^{2h}{\mathbb{V}}[X_{t-h}]+\sigma_{ \text{noise}}^2\sum_{h=0}^{j-1}\alpha^{2h}これが,つねに{\mathbb{V}}[X_t]=\sigma^2となるには、[|\alpha|]<1が必要である。

 

さて以降、|\alpha|<1を条件に課すことにします。

proposition

parameter\alpha,\beta,\sigma_{ \text{\noise}}^2 (|\alpha|<1)で与えられる定常AR(1)過程が、平均[\mu],分散\sigma^2 をもつならば、関係式

\mu=\frac{\beta}{1-\alpha}

\sigma^2=\frac{\sigma_{ \text{noise}}^2}{1-\alpha^2}

が成り立つ。

これを導くには、式①の両辺の平均、分散をとり

\mu=\alpha^h \mu+\beta(1+\alpha+\cdots+\alpha^{h-1})\sigma^2=\alpha^{2h}\sigma^2 +\sigma_{ \bext{noise}}(1+\alpha^2+\cdots+\alpha^{2(h-1)})

を導いてh \to \infty極限を取れば得られます。


さいごに、定常AR(1)過程の自己共分散を調べてこのセクションの締めとします。それを求めるには、まず次の事実を確かめなければなりません
事実:定常AR(1)過程X_t 時刻t,t' (t<t')について

 \text{Cov}(X_t,w_{t'})=0

 

上の事実は、まず式①,(AR2)を使って

 \text{Cov}(X_t,w_{t'})=\alpha^{2h} \text{Cov}(X_{t-h},w_{t'})

となります。

この量は、コーシー・主ワルツ不等式より絶対値でもって

\alpha^{2h} \sqrt{{\mathbb V}[X_{t-h}]{\mathbb V}[w_{t'}]}=\sqrt{\frac{\alpha^{2h}\sigma^4_{\text{noise}}}{1-\alpha^2} }で抑えられているので、h \to \inftyの極限をとれば導ける。

proposotion

parameter\alpha,\beta,\sigma_{ \text{noise}}^2 (|\alpha|<1)で与えられる定常AR(1)過程X_t、および時刻 t,t' (t>t')にたいし自己共分散は次のようになる

 \text{Cov}(X_t,X_{t'})=\frac{\alpha^{2(t-t')}}{1-\alpha^2}\sigma^2_{ \text{noise}}

とくに、AR(1)の自己共分散は定常、かつ時刻差に関して指数減衰する。

(証明)

式でh=t-t'として、 \text{Cov}(X_t,X_{t'})=\alpha^{2h} \text{Cov}(X_{t'},X_{t'})=\frac{\alpha^{2(t-t')}}{1-\alpha^2}\sigma^2_{ \text{noise}}よりわかる。

 

今回導入したAR(1)過程は、さらにAR(p)過程 p:2以上の自然数へと一般化されます。

それらについて調べるには、フーリエ変換が必要なのですが、今回見たように

AR(1)についてはフーリエ変換を用いずとも、初等的に性質を導くことができます。

もちろん、フーリエ変換を使ってAR(1)を調べることもでき、それから自己共分散の

指数減衰則を導くことは、統計力学で有名な ウィーナー・ヒンチンの関係式の再証明であるといえます。

 

以上わかったことを表にまとめると次のようになります。

  ホワイトノイズ ブラウン運動 定常AR(1)過程
定常平均
定常分散
定常増分
定常自己共分散
 \text{Cov}(X_t,X_s)の振る舞い t=sを除き0;  \text{min}{t,s}に関し、線形に増大 |t-s|に関し指数的減衰

トーラスのde Rham cohomology(part1)

 \require{AMScd}
こんにちは。今日からトーラスのde Rhamコホモロジーの話をしようと思います。内容としては,(1)de Rham コホモロジーを導入し次に、(2)Fourie級数を復習、最後に(3)2次元トーラスのde Rham コホモロジーを計算という流れとなります。part1では(1),(2)を扱っていきます.

1.de Rham cohomology群
{\displaystyle X}を可微分多様体とします.
{\displaystyle k}を非負整数とし,{\displaystyle X}上の滑らかな{\displaystyle {\mathbb C}}{\displaystyle k}-formのなす線形空間{\displaystyle {\mathcal A}^k_X}で表します.外微分

{\displaystyle d^k: {\mathcal A}^k_X \rightarrow  {\mathcal A}^{k+1}_X}

{\displaystyle d^{k+1} \circ d^k=0}を満たしているので{\displaystyle 
\{  {\mathcal A}^k_X , d^k\}}は複体になります.
これのk次コホモロジー{\displaystyle \frac{{\rm Ker} \ d^k}{{\rm Im} \ d^{k-1}}}
{\displaystyle X}{\displaystyle {\sf de-Rham \ cohomology}}群といい
{\displaystyle H^k_{\rm dR}(X)}で表します.de Rham cohomology群を素直に計算するには,微分方程式を解けば良いことになります.

例えば,{\displaystyle X={\mathbb R}}の0次de Rhamコホモロジーならば、次のように計算できます.{\phi \in \displaystyle {\rm Ker }\ d^0}微分方程式{\displaystyle \frac{d\phi}{dx}=0}と同値で、その解空間は積分定数だけの自由度分あるから{\displaystyle {\mathbb C}}と同型.{\displaystyle {\rm Im} \ d^{-1}=0}と併せれば{\displaystyle H^0_{\rm dR}({\mathbb R})\cong {\mathbb C}}となります.さらに1次のde Rham cohomology群について,{\displaystyle {\rm Ker}\  d^1={\mathcal A}^1_X}であることがわかります.微分積分学の基本定理から,{\displaystyle {\rm Im} \ d^0={\mathcal A}_X^1}がわかりますから,{\displaystyle H^1_{\rm dR}({\mathbb R})\cong 0}が言えるのです.

次に{\displaystyle X={\mathbb T}^1}(円周)の場合を考えてみましょう.{\displaystyle {\mathbb T}^1={\mathbb R}/{\mathbb Z}}で,特に1次元トーラスとなっています.そのde Rhamコホモロジーを求めるのですから,トーラス上の微分方程式を解くことになります.それがFourie級数を導入するモチベーションとなります.

2.Fourie級数
トーラス上の関数を表示するツールとしてしばしば使われるのがFourie級数です.
{\displaystyle {\mathbb T}^1}上の滑らかな{\displaystyle {\mathbb C}}値関数{\displaystyle \phi(t)}に対する{\displaystyle {\sf Fourie} }級数を次で定める:

\begin{equation} \sum_{m \in {\mathbb Z}} \hat{\phi}(m)e^{2\pi \sqrt{-1}mt}\end{equation}

ここで{\displaystyle \hat{\phi}(m)}{\displaystyle {\sf Fourie}}係数というもので

{\displaystyle \hat{\phi}(m):=\int_{{\mathbb T}^1} dt \  \phi(t)e^{-2\pi \sqrt{-1}mt}}で与えられるものです.

定理:
{\displaystyle {\mathbb T}^1}上の滑らかな{\displaystyle {\mathbb C}}値関数{\displaystyle \phi(t)},および{\displaystyle r=0,1,2,\cdots }に対して

(1)\begin{equation} \phi^{(r)}(t)=\sum_{m \in {\mathbb Z}} (2\pi \sqrt{-1}m)^r\hat{\phi}(m)e^{2\pi \sqrt{-1}mt}\end{equation}
が成り立つ.

(2)対応{\displaystyle {\mathcal F}:\phi \mapsto \{ \hat{\phi} (m)\}_{m \in {\mathbb Z}}},は
{\displaystyle {\mathbb T}^1}上の滑らかな{\displaystyle {\mathbb C}}値関数全体のなす{\displaystyle {\mathbb C}}線形空間
{\displaystyle {\mathcal V}:=\{ \{ a_m \}_{m \in {\mathbb Z}} | a_m \in {\mathbb C}\ ,  a_m={\mathcal O}(|m|^{-h} )(^\forall h >0).\} }
との間の同型を導く.

上の主張は言い換えると次のようになります:
{\mathcal A}_{{\mathbb T}^1} 上の微分作用素\frac{d^r}{dt^r}はFourie変換を通してみると
{\mathcal V}上の対角行列

D:=\begin{bmatrix}
\ddots&&&&&&\\
&(2 \pi \sqrt{-1} (-2) )^r&&&&&\\
&&(2 \pi \sqrt{-1} (-1) )^r&&&&\\
&&&0&&&\\
&&&&(2 \pi \sqrt{-1} )^r&&\\
&&&&&(2 \pi \sqrt{-1} (2) )^r&\\
&&&&&&
\ddots
\end{bmatrix}
に対応する.

業界用語(?)で「Fourie変換で,微分が掛け算に変わる」ってやつです.それでは次回,上の主張を使って1次元トーラスのde Rham cohomology群を計算してみましょう.


同型{\displaystyle \nu:{\mathcal A}^0_{{\mathbb T}^1} \cong {\mathcal A}^1_{{\mathbb T}^1} , \phi \mapsto \phi dt}に注意しましょう.Fourie変換によって次の可換図式が得られます.

\begin{CD}
{\mathcal A}_{{\mathbb T}^1} ^0                 @>{d^0}>> {\mathcal A}_{{\mathbb T}^1} ^1@>{d^1}>>0\\
@V{\mathcal F}VV @VV{{\mathcal F}\circ \nu^{-1}}V 
@VV0V\\
{\mathcal V}@>D>> {\mathcal V}@>0>> 0

\end{CD}

\underline{H^0_{\rm dR}({\mathbb T}^1})
上の可換図式において,垂直方向は同型ですから,
{\displaystyle H^0_{\rm dR}({\mathbb T}^1) \cong {\rm Ker}\ D \cong {\mathbb C}}となります.

\underline{H^1_{\rm dR}({\mathbb T}^1})
これも可換図式から
{\displaystyle H^1_{\rm dR}({\mathbb T}^1) \cong {\rm Coker}\ D \cong {\mathbb C}}となります.

1次元トーラスのde Rham cohomologyが以上のようにFourie級数を使ってできました.次回は2次元トーラスの場合を実行します(続く).