9月5日（金）午前

生存曲線のシミュレーション 4

昨日のコード、ところどころ修正したミスを書き換えてなかった部分があったので、ちゃんとしたものを記録として残しておく。

Surv.time<-function(p,q,l){
cnt<-1
repeat{ 
if (cnt==l){x<- c(cnt,0);break}
else if (runif(1)< p) {x<- c(cnt,1);break}
else if (runif(1)< q) {x<- c(cnt,0);break}
else cnt<-cnt+1
}
return(x)
}# 生存時間のシミュレーション

A.St<-function(p,q,l){
Surv.time(p,q,l)
}

data.sample<-function(p,q,l,N){
mat<-numeric(0)
for(i in 1:N)
mat<-rbind(mat,A.St(p,q,l))
return(mat)
}# N人を対象にしたデータを作成

ri<-function(d,t){
if (length(d[,1][d[,1]>=t])==0)x<-1
else x<-length(d[,1][d[,1]>=t])
return(x)
}

di<-function(d,t){
cnt<-0
for(i in 1:nrow(d)){
if(d[i,1]==t&&d[i,2]==1)cnt<-cnt+1
else cnt<-cnt
}
return(cnt)
}#diを出す関数

S.hat<-function(d,t){
x<-1
for(i in 0:t){
x<-x*(1-di(d,i)/ri(d,i))
}
return(x)
}#カプラン・マイヤー法に適用

S.true<-function(p,t){
(1-p)^t
}

E<-function(p,q,l,N){
d<-data.sample(p,q,l,N)
x<-0
for(i in 1:l)x<-x+(S.hat(d,i)-S.true(p,i))^2
return(x/l)
}

Ev<-function(p,q,l,N){
x<-c()
for(i in 1:1000) x<-c(x,E(p,q,l,N))
return(var(x)*99/100)
}

xq<-seq(0,0.5,0.05)

yq<-function(j){
yq<-c()
for(i in xq)yq<-c(yq,Ev(0.01*j,i,20,50))
return(yq)
}

yq1<-yq(1)
yq2<-yq(2)
yq3<-yq(3)
yq4<-yq(4)
yq5<-yq(5)

plot(xq,yq1,xlim=c(0,0.5),ylim=c(0,0.01), 
ann=F,type="l", lty=1,col=1)

par(new=T)
plot(xq,yq2,xlim=c(0,0.5),ylim=c(0,0.01), 
ann=F,type="l", lty=1,col=2) 

par(new=T) 
plot(xq,yq3,xlim=c(0,0.5),ylim=c(0,0.01), 
ann=F,type="l", lty=1,col=3)

par(new=T) 
plot(xq,yq4,xlim=c(0,0.5),ylim=c(0,0.01), 
ann=F,type="l", lty=1,col=4)

par(new=T) 
plot(xq,yq5,xlim=c(0,0.5),ylim=c(0,0.01), 
,type="l", lty=1,col=5,xlab="q 脱落率" , ylab="Eの分散値" ,main="p,qによるEの分散")

legend(0.01,0.0095,c("p=0.01","p=0.02","p=0.03","p=0.04","p=0.05"),lty=1,col=1:5)

さすがにあの結果はブレがありすぎるということで、先生の助言通り、サンプル数を100から1000にしてみた。（パソコンの限界が近い。）

結果は下の通り。

N=50のとき

N=100のとき

qが大きくなるとEの分散は大きくなるのに、pが大きくなるとEの分散は小さくなるようである。

pとqの変域を対称にしたいところだが、そうするとの値が小さくなりすぎるので、難しい。

Eという指標に問題があるという可能性もあり、あまり本質的ではない問題な気もする。
他にしたいこともあるので、平均の結果が得られただけよしとして、とりあえず一旦保留にしようと思う。

尤度の計算

今日のミーティングで「尤度」の概念は理解した方がよいという話が出たので、この機会に計算してみることにする。

すなわち、データからモデルの復元を試みる。

まず、p、qを適当な乱数で設定し、データをとる。

そこから、最も確からしいp、qが何かを推定してみる。

成果は午後にまとめる。

以上