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

関数間のコサインを使って位相ずれを見つけて補正する

時系列 相関関数 フーリエ変換 畳み込み R Haskell 

# 離散データ
# 補間する(補間方法はいろいろあり得る)
# 保管した結果、regularly sampled discrete time seriesとなったので
# それらについて、離散フーリエ展開を介して、相関関数を求める
Np<-15
ts<-sort(runif(Np))
vs<-sort(runif(Np))
#vs<-sin((1:Np)/10)
vs<-vs-min(vs)
library(pracma)
my.convolve<-function(x,y,type="open"){
	i.p<-(convolve(x,y,type=type))
	#print(i.p)
	x.x<-max(convolve(x,x,type=type))
	#print(x.x)
	y.y<-max(convolve(y,y,type=type))
	#print(y.y)
	return(list(corr=max(i.p)/sqrt(x.x[1]*y.y[1]),vs=i.p/sqrt(x.x[1]*y.y[1])))
}
# vt1, vt2 は値と時刻の２列行列
# n.t は２つの時刻列の重複時間帯を何等分するか
corr.fx<-function(vt1,vt2,n.t=100,method="linear",type="open",permutation=FALSE,n.perm=1000){
	len1<-length(vt1[,1])
	len2<-length(vt2[,1])
	ss<-seq(from=max(vt1[1,2],vt2[1,2]),to=min(vt1[len1,2],vt2[len2,2]),length=n.t)
	x1<-interp1(vt1[,2],vt1[,1],ss,method=method)
	x2<-interp1(vt2[,2],vt2[,1],ss,method=method)
	reserve.vs<-vt2[,1]
	perm.max.id<-0
	perm.max.time<-0
	i.p<-(convolve(x1,x2,type=type))
	x1.x1<-max(convolve(x1,x1,type=type))
	x2.x2<-max(convolve(x2,x2,type=type))
	stat.ori<-max(i.p)/sqrt(x1.x1[1]*x2.x2[1])
	max.id<-which(i.p==max(i.p))
	max.time<-(max.id-n.t)*(ss[2]-ss[1])
	if(permutation){
		stat.perm<-rep(0,n.perm)
		for(i in 1:n.perm){
			shuffle<-sample(1:len2)
			tmp.x2<-interp1(vt2[,2],vt2[shuffle,1],ss,method=method)
			tmp.x2.x2<-max(convolve(tmp.x2,tmp.x2,type=type))
			tmp.i.p<-(convolve(x1,tmp.x2,type=type))
			stat.perm[i]<-max(tmp.i.p)/sqrt(x1.x1[1]*tmp.x2.x2[1])
			tmp.max.id<-which(tmp.i.p==max(tmp.i.p))
			if(stat.perm[i]==max(stat.perm)){
				reserve.vs<-vt2[shuffle,1]
				perm.max.id<-tmp.max.id
				
			}
		}
		p.perm<-length(which(stat.perm>=stat.ori))/n.perm
		perm.max.time<-(perm.max.id-n.t)*(ss[2]-ss[1])
	}else{
		stat.perm=NULL
		p.perm=NULL
		perm.max.vt2=NULL
		perm.max.id=NULL
		perm.max.time=NULL
	}
	return(list(vt1=vt1,vt2=vt2,ss=ss,max.id=max.id,max.time=max.time,n.t=n.t,method=method,type=type,corr=stat.ori,vs=i.p/sqrt(x1.x1[1]*x2.x2[1]),stat.perm=stat.perm,p.perm=p.perm,perm.max.vt2=cbind(reserve.vs,vt2[,2]),perm.max.id=perm.max.id,perm.max.time=perm.max.time,n.t=n.t,permutation=permutation,n.perm=n.perm))
}

ss<-seq(from=0, to=1,length=1000)

ss<-ss[which(ss>(min(ts)) & ss< max(ts))]
yc<-interp1(ts,vs,ss,method="constant")
yl<-interp1(ts,vs,ss,method="linear")
yn<-interp1(ts,vs,ss,method="nearest")
ys<-interp1(ts,vs,ss,method="spline")
ys<-interp1(ts,vs,ss,method="cubic")


plot(ts,vs)
grid()
lines(ss,yc,col=2)
lines(ss,yl,col=3)
lines(ss,yn,col=4)
lines(ss,ys,col=5)
lines(ss,yc,col=6)

#Fyc<-fft(yc)
#Fyl<-fft(yc)
Niter<-1000
ori<-my.convolve(yc,yl)$corr
perm.stat<-rep(0,Niter)
reserve.vs<-vs
for(i in 1:Niter){
	tmpvs<-sample(vs)
	#plot(ts,tmpvs)
	tmpyc<-interp1(ts,tmpvs,ss,method="constant")
	#lines(ss,tmpyc,col=2)
	#lines(ss,yl,col=3)
	perm.stat[i]<-my.convolve(tmpyc,yl)$corr
	if(perm.stat[i]==max(perm.stat,ori)){
		reserve.vs<-tmpvs
	}
}
par(mfcol=c(1,2))
plot(ts,vs)
grid()
tmpyc<-interp1(ts,reserve.vs,ss,method="constant")
lines(ss,tmpyc,col=2)
lines(ss,yl,col=3)
lines(ss,yc,col=4)

plot(sort(perm.stat),ylim=range(c(ori,perm.stat)))
abline(h=ori,col=2)
par(mfcol=c(1,1))


Np1<-100
Np2<-80
ts1<-sort(runif(Np1))*1.5*pi
#vs1<-sort(runif(Np1))
vs1<-sin(ts1)+rnorm(Np1)*0.1
ts2<-sort(runif(Np2))*1.5*pi
#vs2<-sort(runif(Np2))
vs2<-sin(ts2+pi/2)+rnorm(Np2)*0.2
#vs2<-vs2+rnorm(Np2)*0.1
xlim<-range(c(ts1,ts2))
ylim<-range(c(vs1,vs2))

vt1<-cbind(vs1,ts1)
vt2<-cbind(vs2,ts2)

out.perm<-corr.fx(vt1,vt2,permutation=TRUE,n.perm=1000)
par(mfcol=c(2,2))
plot(ts1,vs1,type="l",col=2,xlim=xlim,ylim=ylim,main="Raw Plot")
par(new=TRUE)
plot(ts2,vs2,type="l",col=3,xlim=xlim,ylim=ylim)

plot(ts1,vs1,type="l",col=2,xlim=xlim,ylim=ylim,main="Best-shifted")
par(new=TRUE)
plot(ts2+out.perm$max.time,vs2,type="l",col=3,xlim=xlim,ylim=ylim)

plot(ts1,vs1,type="l",col=2,xlim=xlim,ylim=ylim,main="Best_among_permutations")
par(new=TRUE)
plot(ts2+out.perm$max.time,vs2,type="l",col=3,xlim=xlim,ylim=ylim)
par(new=TRUE)
plot(out.perm$perm.max.vt2[,2]+out.perm$perm.max.time,out.perm$perm.max.vt2[,1],type="l",col=4,xlim=xlim,ylim=ylim)

plot(sort(out.perm$stat.perm),ylim=range(c(out.perm$stat.perm,out.perm$corr)))
abline(h=out.perm$corr,col=2,main="Distribution_of_stats")

par(mfcol=c(1,1))