遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

周期データのパラメタ推定

R 時系列解析 線形代数 複素平面 球 多次元球
  • こちらの続き
  • 複素関数を用いた周期データの変数表現があった
  • 複素数は２つの実数からなる(実部と虚部)
  • 複数の実数に関する回帰・最適化問題は、いろいろなやり方がある
  • 掲載図はoptim()関数を使ってみた例
  • 緑が理論値、黒が実測値(理論値からのずれを正規分布で与えた)、赤が実測値に基づいた推定値
#######################
# 関数
# 多次元単位球面の角座標から、直交座標の座標を作る
sphereCoords<-function(v){
	C<-cos(v)
	S<-sin(v)
	ret<-c(1,S)

	for(i in 1:length(v)){
		for(j in i:length(v)){
			ret[i]<-ret[i]*C[j]
		}
	}
	return(ret)
}
# 各座標から推移行列を作る
TransitionMat<-function(e.val.theta,e.vec.theta,e.vec.psi){
	k<-length(e.val.theta)
	S<-diag(complex(real=cos(e.val.theta),imaginary=sin(e.val.theta)))
	V<-matrix(0,k,k)
	Co<-cos(e.vec.theta)
	Si<-sin(e.vec.theta)
	Comp<-matrix(complex(real=Co,imaginary=Si),k,k)
	V<-Comp
	for(i in 1:k){
		A<-sphereCoords(e.vec.psi[i,])
		V[i,]<-A*Comp[i,]
	}
	return(V%*%S%*%solve(V))
}

# 周期的データを作る
k<-20
e.val.theta<-runif(k)*2*pi
e.vec.theta<-matrix(runif(k^2)*2*pi,k,k)
e.vec.psi<-matrix(runif(k*(k-1))*2*pi,k,k-1)

e.vector.norm<-matrix(runif(k*(k-1))*2*pi,k,k-1)

v<-runif(k)

ret<-sphereCoords(v)

sum(ret^2)

m<-TransitionMat(e.val.theta,e.vec.theta,e.vec.psi)
# 時間を経過させる
Niter<-10

dseries<-matrix(0,Niter,k)
dseries[1,]<-runif(k)+rpois(k,3)
for(i in 2:Niter){
	dseries[i,]<-m%*%dseries[i-1,]
}
# 一部の状態量の時系列変化をプロットする
tobedrawn<-3
matplot(Re(dseries[,1:tobedrawn]),type="l")

#
#initdata<-c(k,e.val.theta,e.vec.theta,e.vec.psi,d.init)
# 完璧周期的データから逸脱させた「実測データ」を作る
dseriesOriginal<-dseries
dseries<-dseriesOriginal + matrix(rnorm(Niter*k),Niter,k)*10
plot(dseriesOriginal[,1],dseries[,1])

# 初期値は「実測データ」に基づかせる
d.init<-Re(dseries[1,])
d.init.theta<-Im(dseries[1,])
#dseries2<-Re(dseries)

# 最適化関数
fn2<-function(x){
	e.val.theta<-x[1:(k)]
	e.vec.theta<-matrix(x[(k+1):(k+k^2)],k,k)
	e.vec.psi<-matrix(x[(k+1+k^2):(k+k^2+k*(k-1))],k,(k-1))
	d.init<-x[(k+1+k^2+k*(k-1)):(k+k^2+k*(k-1)+k)]
	d.init.theta<-x[(k+1+k^2+k*(k-1)+k):(k+k^2+k*(k-1)+k+k)]
	m<-TransitionMat(e.val.theta,e.vec.theta,e.vec.psi)
	est.data<-matrix(0,nrow(dseries),k)
	est.data[1,]<-complex(real=d.init,imaginary=d.init.theta)
	for(i in 2:length(est.data[,1])){
		est.data[i,]<-m%*%est.data[i-1,]
	}
	sum((Re(est.data)-Re(dseries))^2)
}


# 最適化処理の初期データベクトル
x<-c(e.val.theta,e.vec.theta,e.vec.psi,d.init,d.init.theta)

# 最適化
outOptim<-optim(x,fn2)
#outNlm<-nlm(fn2,x)

# 最適化推定結果取り出し

E.x<-outOptim$par
E.x2<-outNlm$estimate
#E.x2<-E.x
E.e.val.theta<-E.x[1:k]
E.e.vec.theta<-matrix(E.x[(k+1):(k+k^2)],k,k)
E.e.vec.psi<-matrix(E.x[(k+1+k^2):(k+k^2+k*(k-1))],k,(k-1))
E.d.init<-E.x[(k+1+k^2+k*(k-1)):(k+k^2+k*(k-1)+k)]
E.d.init.theta<-E.x[(k+1+k^2+k*(k-1)+k):(k+k^2+k*(k-1)+k+k)]
E.m<-TransitionMat(E.e.val.theta,E.e.vec.theta,E.e.vec.psi)
E.dseries<-dseries
E.dseries[1,]<-E.d.init
for(i in 2:length(dseries[,1])){
	E.dseries[i,]<-E.m%*%E.dseries[i-1,]
}

plot(Re(dseries[,1]),Re(E.dseries[,1]))
plot(Re(dseriesOriginal[,1]),Re(E.dseries[,1]))

ylim<-c(min(Re(dseries[,1]),Re(E.dseries[,1])),max(Re(dseries[,1]),Re(E.dseries[,1])))
plot(Re(dseries[,3]),type="l",ylim=ylim)
par(new=TRUE)
plot(Re(E.dseries[,3]),type="l",col=2,ylim=ylim)
par(new=TRUE)
plot(Re(dseriesOriginal[,3]),type="l",col=3,ylim=ylim)