9月5日（金）午後

ガリガリ計算して確率の対数をとって、もっともらしいpとqの値をそれぞれ0.01～0.05まで0.0050刻みで選べる関数をつくった。

p<-runif(1,0.01,0.05)
q<-runif(1,0.01,0.05)

Surv.time<-function(p,q,l){
cnt<-1
repeat{ 
if (runif(1)< p) {x<- c(cnt,1);break}
else if (runif(1)< q) {x<- c(cnt,0);break}
else if (cnt==l){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人を対象にしたデータを作成

d<-data.sample(p,q,20,100)

pt<-function(p,q,t){
log(prod((((1-p)*(1-q))^((1:t)-1))*p))
}# tに死亡が起こる確率

qt<-function(p,q,t){
if(t==20)x<-log((((1-p)*(1-q))^((1:t)-1))*(1-p))
else x<-log(prod((((1-p)*(1-q))^((1:t)-1))*(1-p)*q))
return(x)
}# tに脱落が起こる確率

dp<-function(d){
subset(d,d[,2]==1)
}

dq<-function(d){
subset(d,d[,2]==0)
}

pa<-function(p,q,d){
x<-c()
for(i in 1:length(dp(d)[,1]))x<-c(x,pt(p,q,dp(d)[i,1]))
return(sum(x))
}

qa<-function(p,q,d){
x<-c()
for(i in 1:length(dq(d)[,1]))x<-c(x,qt(p,q,dq(d)[i,1]))
return(sum(x))
}

A<-function(p,q,d){
pa(p,q,dp(d))+qa(p,q,dq(d))
}

mat<-matrix(nrow=9,ncol=9)
ps<-(seq(0.01,0.05,0.005))
qs<-(seq(0.01,0.05,0.005))

for(i in 1:9)
for(j in 1:9)
mat[i,j]<-A(ps[i],qs[j],d)

ex<-function(m){
{if(which(mat==max(mat))%%9==0)x<-ps[9]
else x<-ps[which(mat==max(mat))%%9]}
y<-qs[(which(mat==max(mat))%/%9)+1]
return(c(x,y))
}

p

q

mat

ex(mat)

> p
[1] 0.02417261
> 
> q
[1] 0.0497098
> 
> mat
           [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
 [1,] -3299.853 -3150.888 -3056.414 -2991.888 -2946.378 -2914.053 -2891.436
 [2,] -3252.458 -3103.493 -3009.019 -2944.493 -2898.983 -2866.658 -2844.041
 [3,] -3231.308 -3082.343 -2987.869 -2923.343 -2877.832 -2845.508 -2822.890
 [4,] -3224.638 -3075.673 -2981.199 -2916.673 -2871.162 -2838.838 -2816.220
 [5,] -3227.208 -3078.244 -2983.770 -2919.243 -2873.733 -2841.409 -2818.791
 [6,] -3236.227 -3087.262 -2992.788 -2928.261 -2882.751 -2850.427 -2827.809
 [7,] -3250.025 -3101.060 -3006.587 -2942.060 -2896.550 -2864.226 -2841.608
 [8,] -3267.531 -3118.566 -3024.092 -2959.566 -2914.055 -2881.731 -2859.113
 [9,] -3288.010 -3139.045 -3044.571 -2980.045 -2934.534 -2902.210 -2879.592
           [,8]      [,9]
 [1,] -2876.281 -2867.059
 [2,] -2828.886 -2819.664
 [3,] -2807.736 -2798.514
 [4,] -2801.066 -2791.844
 [5,] -2803.637 -2794.415
 [6,] -2812.655 -2803.433
 [7,] -2826.454 -2817.232
 [8,] -2843.959 -2834.737
 [9,] -2864.438 -2855.216
> 
> ex(mat)
[1] 0.025 0.050

できるだけfor-loopを使わないようにしていたら、なんだかすごくつくるのに時間がかかった気がする。
処理を早くできたのかはわからないが…。
（前回のコードも書き直そうかと思ったが、for-loopをなくせそうなところがすぐには見つからなかったので、追究するのはやめにした。）

ちなみにこれ、

> exp(sum(mat))
[1] 0

となってしまうので、尤度そのものを出すことは出来なかった。
しかし、尤度の計算の要領はつかめたと思うので、今回はこれでよしとしようと思う。

次は、ケース・コントロール研究のシミュレーションをしていく予定である。

以上