修理する６　指数行列による初期値問題の解 - ryamadaの遺伝学・遺伝統計学メモ

微分行列が以下のようになった
$A=\begin{pmatrix}-p & q & 0 & 0 & 0 &...\\ p & -(p+q) & q & 0 & 0 ... \\0 & p & -(p+q) & q & 0 ... \\ 0 & 0 & p & -(p+q) & q & ... \\ 0& 0& 0 & p & -(p+q) & ...\end{pmatrix}$
実際、これは $\infty \times \infty$ 行列である
現実的には、逆イベントがかなりの確率で起きるなら、vには現実的な上限があるだろうし
vが、「複数個所に起きうる」イベントで、「起きた箇所数」を数えているとしたら、vには真の上限がある
その場合には、行列の最終行・最終列のところが書き換えられて、以下のようにできる
$A=\begin{pmatrix}-p & q & 0 & 0 &...& 0 & 0\\ p & -(p+q) & q & 0 & ... & 0& 0 \\0 & p & -(p+q) & q & ... & 0& 0 \\ 0 & 0 & p & -(p+q) & ... & 0& 0\\ 0& 0& 0 & 0&... & p & -q \end{pmatrix}$
このような行列は、固有値分解を通じて得られる、指数行列を用いて、初期値問題に利用できる(こちら)

# k<-p, h<-q
k<-5
h<-7.5

N<-100
M<-matrix(0,N,N)

M[1,1]<--k
M[1,2]<-h
M[2,1]<-k
M[2,2]<--(k+h)
M[2,3]<-h

for(i in 2:(N-1)){
	M[i,i-1]<-k
	M[i,i]<--(k+h)
	M[i,i+1]<-h
}

M[N,N-1]<-k
M[N,N]<--h

M

n<-N
A<-M
# 固有値分解
eigen.out<-eigen(A)
V<-eigen.out[[2]]
U<-solve(V)
B<-diag(eigen.out[[1]])
# 計算する
Xmax<-1
x<-seq(from=0,to=1,length=101)*Xmax
Ys<-matrix(0,length(x),n)
# 初期値も適当に
Ys0<-rep(0,n)
Ys0[1]<-1
# 複素成分の自乗の総和がほぼゼロであることを確認
IMS<-rep(0,length(x))
for(i in 1:(length(x))){
	Bex<-diag(exp(eigen.out[[1]]*x[i]))
	Aex<-V%*%Bex%*%U
	Ys[i,]<-Aex%*%Ys0
	ims<-sum(Im(Ys[i,])^2)
	IMS[i]<-ims
# 複素成分が大きかったら表示させるようにして確認
	if(ims>0.000001)print(ims)
}
plot(IMS)
par(mfcol=c(1,2))
matplot(Ys,type="l")

Ys1<-Ys


# ただのポアッソン過程と比べてみよう
# k<-p, h<-q
v<-0:(N-1)
t<-x
lambda<-max(k-h,0.00001)
Ypois<-matrix(0,length(t),length(v))
for(i in 1:length(t)){
	Ypois[i,]<-dpois(v,lambda*t[i])
}

matplot(Ypois,type="l")
par(mfcol=c(1,1))