アレル頻度を周期関数で

アレル頻度を空間中に関数で表現するための第１歩として、三角関数を用いてみる。
そのごく初期段階のお試し。

df<-1 # 空間の次元


Nd<-4 # サンプルセット数

Nm<-1 # マーカー数

#Means<-matrix(rep(0,Nd*df),Nd,df) # サンプルセットのサンプリングの中心
Means<-matrix(rnorm(Nd*df),Nd,df)
SDs<-matrix(rep(1,Nd*df),Nd,df) # サンプルセットのサンプリングのSD

# マーカーのアレル頻度は各軸について、1/2*(sin((x-t)/k) +1) で0-1の範囲をとる周期関数とする
# kがサンプリングの分布の中心とSDに比べて、十分に大きいときは、１峰性の関数とみなすことができる

AfCalc<-function(k,t,x){
 1/2*(sin((x-t)/k)+1)	
}

Ks<-matrix(rep(1,Nm*df),Nm,df)
Ts<-matrix(rep(0,Nm*df),Nm,df)
L<-1000
Xs<-seq(from=-10,to=10,length.out=L)

res<-c()
probs1<-probs2<-matrix(0,Nd,Nm)
for(im in 1:Nm){
for(id in 1:Nd){
 DfCounter<-rep(0,df)
 loop<-TRUE
 while(loop){
  DfCounter[1]<-DfCounter[1]+1
  if(df>1){
  for(i in 1:(df-1)){
   if(DfCounter[i]==(L+1)){
    DfCounter[i]<-0
    DfCounter[i+1]<-DfCounter[i+1]+1
   }
  }	
  }
 af<-1
 sp<-1
 for(i in 1:df){
  af<-af*AfCalc(k=Ks[im,i],t=Ts[im,i],x=Xs[DfCounter[i]])
  sp<-sp*dnorm(x=Xs[DfCounter[i]],mean=Means[id,i],sd=SDs[id,i])
  
 }
 probs1[id,im]<-probs1[id,im]+af*sp
  probs2[id,im]<-probs2[id,im]+(1-af)*sp
  res<-c(res,af*sp)
  if(sum(DfCounter)==L*df){
   loop<-FALSE
  }	
 }


}
}

plot(res)

probFinal<-probs1/(probs1+probs2)

probFinal

plot(Xs,AfCalc(k=Ks[1,1],t=Ts[1,1],x=Xs))