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

常微分方程式を勉強する　２

MIKU R 常微分方程式 

library(rgl)

v1<-c(0,0,0)
v2<-c(1,0,0)
v3<-c(1,1,0)
v4<-c(0,1,0)

alpha<-0.3
beta1<-0.2
beta2<-0.7


v5<-c(alpha,beta1,0)
v6<-c(alpha,beta2,0)

L<-1.5

v7<-c(alpha,beta2,L)

Vs<-rbind(v1,v2,v3,v4,v5,v6,v7)

plot3d(Vs,xlim=c(-0.5,1.5),ylim=c(-0.5,1.5),zlim=c(-0.5,1.5))



segments3d(Vs[c(1,2,2,3,3,4,4,1,5,6,6,7,7,5),1],Vs[c(1,2,2,3,3,4,4,1,5,6,6,7,7,5),2],Vs[c(1,2,2,3,3,4,4,1,5,6,6,7,7,5),3],col=2)


Ntri<-20

for(i in 1:Ntri){
	tmpalpha<-runif(1)
	tmpbeta1<-runif(1)
	tmpbeta2<-runif(1)
	tmpbeta1.2<-min(tmpbeta1,tmpbeta2)
	tmpbeta2.2<-max(tmpbeta1,tmpbeta2)
	tmpL<-runif(1)*L
	
	tmpv5<-c(tmpalpha,tmpbeta1.2,0)
	tmpv6<-c(tmpalpha,tmpbeta2.2,0)
	tmpv7<-c(tmpalpha,tmpbeta2.2,tmpL)
	tri<-rbind(tmpv5,tmpv6,tmpv7)
	triangles3d(tri)
}

常微分方程式を勉強する

MIKU R 常微分方程式

講座 数学の考え方〈7〉常微分方程式論

講座 数学の考え方〈7〉常微分方程式論


C1<-runif(1)
C2<-runif(1)
as<-seq(from=-1,to=1,length=11)

xs<-seq(from=-10,to=10,length=101)

# ys1,ys2は第１、第２導関数値
ys<-ys1<-ys2<-matrix(0,length(as),length(xs))

for(i in 1:length(ys[,1])){
	a<-as[i]
	absa<-abs(a)
	if(a<0){
		ys[i,]<-C1*sinh(xs*sqrt(absa))+C2*cosh(xs*sqrt(absa))
		# 微分しよう
		dx<-deriv(~C1*sinh(xs*sqrt(absa))+C2*cosh(xs*sqrt(absa)),"xs")
		# 値の評価
		ys1[i,]<-eval(dx)
	}else if(a>0){
		ys[i,]<-C1*sin(xs*sqrt(absa))+C2*cos(xs*sqrt(absa))
		dx<-deriv(~C1*sin(xs*sqrt(absa))+C2*cos(xs*sqrt(absa)),"xs")
		ys1[i,]<-eval(dx)
	}else{
		#ys[i,]<-C1+C2*xs
		ys[i,]<-C1*xs+C2
		dx<-deriv(~C1*xs+C2,"xs")
		ys1[i,]<-eval(dx)
	}
	ys2[i,]<-(-a)*ys[i,]
}

matplot(t(ys),type="l",ylim=range(ys[11,])*5)

library(rgl)
plot3d(xs,ys[1,],ys1[1,])

# a>0 の範囲では
# ysとその第１導関数は周期が半周ずれていて、第２導関数の周期はもとに戻っている

xys<-cbind(xs,ys[8,],ys1[8,],ys2[8,])
plot(as.data.frame(xys))
matplot(xys[,2:4],type="l")

# a<0 の範囲では
# 発散する
xys<-cbind(xs,ys[3,],ys1[3,],ys2[3,])
plot(as.data.frame(xys))
matplot(xys[,2:4],type="l")

与えられた条件から、推定する〜信頼区間？〜

MIKU 統計 R

数学いらずの医科統計学 第2版

数学いらずの医科統計学 第2版

  • 関連記事など(こちら)
  • 条件とデータが与えられたときに、知りたい変数の値を推定したい
  • 最尤推定
    • 最も尤もらしい値を推定する(最尤推定値)
    • 推定したい変数の尤度関数
    • 観察データから尤度関数の値を求める
    • 最尤推定値は、尤度関数の最大値を与える変数値。微分できる？→微分する
  • 尤もらしい範囲を推定する
    • 尤度に基づくありそうな値とありそうもない値とがわかる
    • 尤もらしい値を範囲で選ぶとしたらどうしたいか
    • Rの信頼区間(こちら)
  • メモ

p<- seq(0,1,0.00001)
pp<- rep(0,length(p))
for(i in 1:length(p)){
	pp[i]<- dbinom(4,15,p[i])
}
plot(pp)


#plot(pp/sum(pp),cex=0.1,axes=FALSE,ylim=c(0,0.0005),
#	 xlab="Ratio of Red balls",ylab="Likelyhood")
#axis(1,
#	 c(0,which(pp==max(pp)),length(p)),
#	 c(0,round(which(pp==max(pp))/length(p),3),1))
#axis(2)

a<- pp/sum(pp)
a_order<- order(a,decreasing=TRUE)

a_ordered<-a[a_order]
cumsum_a_ordered<-cumsum(a_ordered)
plot(cumsum_a_ordered)

closest<-which(abs(cumsum_a_ordered-0.95)==min(abs(cumsum_a_ordered-0.95)))
abline(v=closest)
abline(h=0.95)

plot(pp,type="l")
#abline(h=a[a_order[closest]])

thres<-a[a_order[closest]]
mores<-which(a>thres)
lessANDeq<-which(a<=thres)

sum(a[mores])
sum(a[lessANDeq])

CI<-qbeta(c(1-0.95)/2,1-(1-cirange)/2),4+1,15-4+1)
abline(v=CI*length(p),col=2)
abline(v=a_order[closest])
abline(v=a_order[closest-1])


N<-15
k<-4
p<-seq(from=0,to=1,by=0.001)
v<-dbeta(p,k+1,N-k+1)
v<-v/sum(v)
plot(p,v,type="l")
abline(v=k/N)
cirange<-0.95
ci1<-qbeta(c((1-cirange)/2,1-(1-cirange)/2),k+1,N-k+1)
abline(v=ci,col=2)

lesses<-which(p<ci1[1])
mores<-which(p>ci1[2])


a_order<- order(v,decreasing=TRUE)

a_ordered<-v[a_order]
cumsum_a_ordered<-cumsum(a_ordered)

closest<-which(abs(cumsum_a_ordered-0.95)==min(abs(cumsum_a_ordered-0.95)))
thres<-v[a_order[closest]]
mores2<-which(v>thres)
lesss2<-which(v<thres)
col<-rep(0,length(p))
col[mores2]<-2
col[lesss2]<-3
par(new=TRUE)
plot(p,v,col=col,cex=0.1)
abline(h=thres,col=3)
sum(v[c(lesses,mores)])
sum(v[-c(lesses,mores)])
sum(v[c(lesss2)])
sum(v[c(mores2)])
  • dbinom()はp引数がベクトル対応なのでループがはしょれるらしい(こちら)
    • for()は、好きかもしれない