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