# xk=sin(vk) # xk-1=cos(vk)*sin(vk-1) # xk-2=cos(vk)*cos(vk-1)*sin(vk-2) # ... # x2=cos(vk)*cos(vk-1)*cos(vk-2)*...*sin(v1) # x1=cos(vk)*cos(vk-1)*cos(vk-2)*...*cos(v1) Niter<-1000 n<-3 xs<-matrix(0,Niter,n) xs[1,]<-runif(n) ps<-runif(n-1)*10/(2*pi) dt<-0.01 for(i in 2:Niter){ tmp<-Cartesian2Angular(xs[i-1,]) newang<-tmp$t+dt*ps tmp2<-Angular2Cartesian(tmp$r,newang) xs[i,]<-tmp2 } matplot(as.data.frame(xs)) t<-seq(from=0,to=2,length=Niter)*2*pi xs[,1]<-cos(t) xs[,2]<-sin(t) xs[,3]<-0 rot<-NormalBase(3) #rot<-diag(rep(1,3)) ts<-matrix(0,Niter,2) rs<-rep(0,Niter) newxs<-xs for(i in 1:Niter){ xs[i,]<-rot%*%xs[i,] tmp<-Cartesian2Angular(xs[i,]) ts[i,]<-tmp$t rs[i]<-tmp$r newxs[i,]<-Angular2Cartesian(rs[i],ts[i,]) } plot(ts) plot(xs[,1],newxs[,1]) plot(xs[,2],newxs[,2]) plot(xs[,3],newxs[,3]) Cartesian2Angular(c(0,1/sqrt(2),1/sqrt(2))) # 各座標からデカルト座標へ Angular2Cartesian<-function(r,t){ n<-length(t)+1 #tmpt<-(t+pi)%%(2*pi)-pi C<-abs(cos(t)) S<-sin(t) C2<-cumprod(C[length(C):1]) x<-c(1,S) x[1:(n-1)]<-x[1:(n-1)]*C2[length(C):1] x[1]<-x[1]*sign(cos(t[1])) x<-x*r x } # 検算する x<-runif(4) a<-Cartesian2Angular(x) x2<-Angular2Cartesian(a$r,a$t) x a x2 # n-th座標の正負を無視してもよければ、こんな方法も c Theta<-function(x,v){ acos(sum(x*v)/sqrt(sum(x^2)*sum(v^2))) } n<-6 Es<-diag(rep(1,n)) x<-c(rep(0,n-1),1) x<-x[sample(1:n)] x<-rnorm(n) thetas<-rep(0,n) for(i in 1:n){ thetas[i]<-Theta(x,Es[i,]) } Angular2Cartesian2<-function(r,t){ x<-c(cos(t),0) x[length(x)]<-sqrt(1-sum(x^2)) x<-r*x x } x2<-Angular2Cartesian2(sqrt(sum(x^2)),thetas[1:(n-1)]) x x2
遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典