こちらでperspを使って描いているけど、plot3dでグリグリしたいので
N<-50 #集団のサイズ k<-3 #残せる遺伝子の最大数 #k<-"infinity" #残せる遺伝子が無限大の場合 T<-10 #世代数 initial <- 5 #最初に存在する変異アレル数 #推移確率の行列を作る A<-matrix(0,N+1,N+1) for (i in 0:N){ for(j in 0:N){ if (k == "infinity"){ A[j+1,i+1]<-choose(N,j)*i^j*(N-i)^(N-j)/N^N }else{ A[j+1,i+1]<- choose(i*k,j)*choose((N-i)*k,N-j)/choose(N*k,N) } } } #apply(A,1,sum) #確認 #各世代ごとの変異アレル数の存在確率のベクトル v<-c(rep(0,N+1)) v[initial+1]<-1 #第１世代の分布 #各世代のベクトルvを並べて行列Dを作る D<-c() D<-cbind(D,v) for (t in 1:T){ v<-A%*%v #次世代のベクトル D<-cbind(D,v) } #apply(D,2,sum) #確認 persp(D,theta=150,phi=30,shade=0.2,xlab="mutatns",ylab="generations",zlab="probability",col="green")
で、ここからplot3d用
x<- 1:(N+1) # x点を用意 y<- 1:(T+1) # y点を用意 library(rgl) plot3d(x,y,D,xlab="Mutatns",ylab="Generations",zlab="Probability")
こんな感じで大丈夫か?
なんか違うのであとで修正する。
やっぱりな。
open3d() plot3d(x,y,D[,1])
でなんか傾いていることが判明。
open3d() plot3d(x,1,D[,1])
ならうれしい。
行列を強引に(x,y,z)に持っていくことにしよう。
今、上でN,T,Dが生成されたとして
x<- t(t(c(rep(1:(N+1),T+1)))) # x点の生成 y<- c() for(i in 1:(T+1)){ y<- c(y,rep(i,N+1)) } y<- t(t(y)) # y点の生成 z<- c(D) # z点の生成 library(rgl) plot3d(x,y,z,xlab="Mutatns",ylab="Generations",zlab="Probability",col=rainbow((T+1)*(N+1)))
