確率論的リスク解析の数理と方法 コロナ社 金野秀敏先生

#p.65 多重ワイブル分布　最尤法
library(survival)

data=data.frame(lung)
data=na.omit(data)

t=data$time
t0=min(t)
status=data$status-1

x=cbind(rep(1,nrow(data)),data$age,ifelse(data$sex==1,1,0),ifelse(data$sex==2,1,0),data$ph.karno,data$pat.karno,data$meal.cal/100,data$wt.loss)

ite=5*10^4
eta=10^(-6)
beta=rep(10^(-9),ncol(x))
m=0.5

for(l in 1:ite){
  
lam=log(t^m/t0)+x%*%beta  
ptx=1-exp(-exp(lam))  

dptx1=t(x)%*%c(exp(lam)*exp(-exp(lam))*status/ptx)
dptx2=-t(x)%*%c(exp(lam)*exp(-exp(lam))*(1-status)/(1-ptx))
dptx=dptx1+dptx2
dm=sum(c(log(t)*exp(lam)*exp(-exp(lam))*status/ptx)-c(log(t)*exp(lam)*exp(-exp(lam))*(1-status)/(1-ptx)))
m=m+eta*dm
beta=beta+eta*dptx
loglik=sum(status*log(ptx))+sum((1-status)*log(1-ptx))

#print(loglik)
}

avet=(1/m)*(log(t0)-x%*%beta)+log(gamma(1+1/m))

accuracyvec=c()
values=seq(0.1,0.9,0.01)

for(j in 1:length(values)){
accuracy=sum(ifelse(ptx>values[j],1,0)==status)/length(status)
accuracyvec=c(accuracyvec,accuracy)
}

plot(values,accuracyvec)

#SIRの離散近似計算　p.97
#確認中
t=c(1:20)
I=c(2,3,4,5,6,7,8,9,10,9,7,5,4,3,2,1,1,1,1,1)
S=c(18,17,16,15,14,13,12,11,10,10,rep(10,10))
R=20-(I+S)
r0=1

N=S[1]-S

beta=N[length(N)]/sum(S*I)
rinv=sum(I)/(R[length(R)]-r0)
r=1/rinv

#Pr(dN(t)=1,dR(t)=0|St-)
beta*S*I

#Pr(dN(t)=0,dR(t)=1|St-)
r*I

#sig(beta)
beta/sqrt(N[length(N)])

#sig(rinv)
rinv/sqrt(R[length(R)]-r0)