一般化線形モデルと生存分析（朝倉書店）・サービスサイエンスの事訳　ーデータサイエンスと数理科学の融合ー　（筑波大学出版）

```{r}

#packageを使わないでやってみた(生存時間解析)

m=80
lambda=function(t){
z<-exp(-((t-m)^2)/100)*0.2
return(z)}

#plot(lambda(seq(0,50,0.1)))

#exp(-integrate(lambda,0,10)$value)

#plot(exp(-seq(0,20,0.1)))

C=exp(-integrate(lambda,0,Inf)$value)

S=function(x){

func=function(t){
  y<-lambda(t)
  return(y)
}

a=0;H=x;n=100000

Y=(H/3)*(1/(n-1))*(func(a)+func(a+H));

for(j in 1:(n-2)){
  coefficient=ifelse(j-floor(j/2)*2==0,4,2)
  Y=Y+((H*coefficient)/(3*(n-1)))*func(a+(j*H)/(n-1))

}
return(exp(-Y)-C)}

#discrete_times
t=seq(0,100,1)

hazard_randoms=data.frame(t=t,prob=0,dis=0)

for(j in 2:nrow(hazard_randoms)){

hazard_randoms$prob[j]=S(hazard_randoms$t[j-1])-S(hazard_randoms$t[j])

hazard_randoms$dis[j]=S(hazard_randoms$t[j])

}

fun=function(s){
z<-log(lambda(s))+log(S(s))
return(z)}

fun_data=data.frame(t=t,value=0)

for(j in 1:nrow(fun_data)){
fun_data$value[j]=fun(fun_data$t[j])
}

times=100
newton_raphson=data.frame(times=1:times,value=0,f=0,df=0)
value=0
for(j in 1:nrow(newton_raphson)){
f<-(-(value-80)-exp(-((value-80)^2)/100)*10)
df<-((value-80)/5)*exp(-((value-80)^2)/100)-1

  newton_raphson$f[j]=f
  newton_raphson$df[j]=df
  value=-f/df+value

  newton_raphson$value[j]=value
}

solution=tail(newton_raphson$value,1)

lambda_hat=lambda(solution)

S_hat=S(solution)

likelihood_data=data.frame(t=seq(0,150,1),ratio=0)

for(j in 1:nrow(likelihood_data)){
  
likelihood_data$ratio[j]=(lambda(likelihood_data$t[j])*S(likelihood_data$t[j]))/(lambda_hat*S_hat)  
  
}

plot(likelihood_data$t,likelihood_data$ratio,type="o",col=3,xlab="t",ylab="ratio",main="ratio_values")

alpha=0.95

num=round((nrow(likelihood_data)*(alpha)))

t_max=likelihood_data$t[likelihood_data$ratio==max(likelihood_data$ratio)]

likelihood_data_sub=likelihood_data

for(j in 1:num){
  
likelihood_data_sub=likelihood_data_sub %>% dplyr::filter(ratio!=min(likelihood_data_sub$ratio))  
  
  
}

under=min(likelihood_data_sub$ratio)

print(paste0("精度",alpha*100,"%の尤度比の値は",under))

sample_t=sample(60:80,100,replace=T)

t_data=data.frame(t=unique(sample_t),freq=0)

S_value=c();lambda_value=c()

for(j in 1:nrow(t_data)){


t_data$freq[j]=length(sample_t[sample_t==t_data$t[j]])/length(sample_t)

S_value=c(S_value,S(t_data$t[j]))

lambda_value=c(lambda_value,lambda(t_data$t[j]))

}

KL=t_data$freq*log(t_data$freq/(S_value*lambda_value))

#KLは0に近いほどいい

KL=sum(KL)

N=100

X=cbind(rnorm(N,mean=0.15,sd=0.01),rnorm(N,mean=20,sd=3))

Y=c()

for(j in 1:N){

y_value=ifelse(X[j,1]>=0.15,ifelse(X[j,2]>=20,rnorm(1,mean=0.5,sd=0.01),rnorm(1,mean=0.25,sd=0.01)),ifelse(X[j,2]>=20,rnorm(1,mean=0.75,sd=0.01),rnorm(1,mean=0.25,sd=0.01)))

Y=c(Y,y_value)

}

log_Y=log(Y)

beta=solve(t(X)%*%X)%*%t(X)%*%log_Y

alpha=mean(log_Y)-sum(apply(X,2,mean)*beta)

predict=c()

for(j in 1:length(Y)){

predict=c(predict,exp(alpha+sum(X[j,]*beta)))

}

plot(c(1:length(Y)),Y,type="o",col=2,ylab="values",xlab="num",ylim=c(0.1,1.0),main="cox partial likelihood")

par(new=T)

plot(c(1:length(Y)),predict,type="o",col=3,ylab="values",xlab="num",ylim=c(0.1,1.0),main="cox partial likelihood")

print(paste0("相関は",cor(log(predict),log_Y)))

#baseline survivor function

S_t_x=function(x,t){

z<-S(t)

if(z>((10^(-3)))){

z<-z^exp(alpha+sum(beta*x))

}else{

z<-0

}

return(z)

}

S_t_x_data=data.frame(t=seq(0,150,1),x1=0.3,x2=20,S_t_x=0,S=0)

for(j in 1:nrow(S_t_x_data)){

S_t_x_data$S_t_x[j]=S_t_x(c(S_t_x_data$x1[j],S_t_x_data$x2[j]),S_t_x_data$t[j])

S_t_x_data$S[j]=S(S_t_x_data$t[j])

}

plot(S_t_x_data$t,S_t_x_data$S_t_x,type="o",col=2,xlab="t",ylab="S",main="baseline survivor function")

plot(S_t_x_data$t,S_t_x_data$S,type="o",col=3,xlab="t",ylab="S",main="baseline survivor function")




```

```{r}
#指数分布;S(t)

data=data.frame(i=c(1:21),t=c(6,6,6,6,7,9,10,10,11,13,16,17,19,20,22,23,25,32,32,34,35),delta=c(1,1,1,0,1,0,1,0,0,1,1,0,0,0,1,1,0,0,0,0,0))

p_hat=sum(data$delta)/sum(data$t)

#weibull分布:S(t)

d=sum(data$delta)

delta=data$delta

t=data$t

times=100
newton_raphson=data.frame(times=1:times,value=0,f=0,df=0)
value=10
for(j in 1:nrow(newton_raphson)){
f<-(d/value+sum(delta*log(t))-d*sum((t^value)*log(t))/sum(t^value))

df<-(-d/(value^2)-d*(sum(value*t^(value-1))/sum(t^value)+sum(log(t)*t^value)*sum(value*t^(value-1))/(sum(t^value)^2)))
  
  newton_raphson$f[j]=f
  newton_raphson$df[j]=df
  value=-f/df+value

  newton_raphson$value[j]=value
}

value=tail(newton_raphson$value,1)

k=value;p=(d/sum(t^k))^(1/k)

#kaplan-meier estimator

data=data[order(data$t),]

t_data=data %>% group_by(t) %>% dplyr::summarise(d=sum(delta))

t_data=t_data %>% dplyr::filter(d!=0)

t_data=rbind(data.frame(t=0,d=0),t_data)

t_data=t_data %>% dplyr::mutate(r=0,c=0)

t_data=t_data %>% dplyr::mutate(t1=c(t_data$t),t2=c(t_data$t[2:nrow(t_data)],NA))



for(j in 1:nrow(t_data)){
  


if(j==nrow(t_data)){
  
data_sub=data %>% dplyr::filter(t<=max(data$t) & t>=t_data$t1[j])  
  
}else{
  
data_sub=data %>% dplyr::filter(t<t_data$t2[j] & t>=t_data$t1[j])  

}  
  
t_data$c[j]=length(data_sub$delta[data_sub$delta==0])

}

for(j in 1:nrow(t_data)){
  
t_data$r[j]=nrow(data)-sum(c(t_data$d[1:j],t_data$c[1:j]))  
  
}

t_data$r=c(t_data$r[1],lag(t_data$r)[2:nrow(t_data)])

t_data=t_data %>% dplyr::mutate(ratio=1-d/r) %>% dplyr::mutate(S=cumprod(ratio))



```
```{r}

#Cox部分尤度法 勾配上昇法 p.157

x1=c(rep(0,10),1,0,1,1,1,0,0,1,0,1,1,1,0,0,0,1,1,0,1,1,1,0,0,0,0,0,1,1,0,0,0,0,0,0,1,1,0,0,0,0,1,0,0,1,0,0,1,0,0,1,0,1,1,1,0,1,0,1,1,0,1,1,0,1,0,0,0,1,0,0)

x2=c(2,3,0,3,2,1,1,1,3,0,2,0,2,0,2,3,1,2,2,1,1,2,1,3,3,1,1,2,3,1,1,3,3,3,1,3,0,1,0,2,3,2,0,2,1,2,0,3,1,1,3,3,3,2,0,3,3,2,3,3,0,3,0,2,2,2,3,1,3,0,0,3,3,3,0,2,0,1,3,0)

x3=c(1,0,1,1,1,1,0,1,0,1,0,0,1,0,1,0,1,1,1,0,1,0,1,0,0,1,1,1,0,0,1,1,1,1,1,1,1,0,1,1,1,0,0,1,0,1,1,1,0,0,0,1,0,0,1,0,1,1,1,1,1,1,0,1,0,1,1,0,1,1,0,1,0,0,1,1,0,0,0,1)

t=c(15,13,15,10,3,2,2,2,1,18,4,5,2,6,1,1,22,11,19,1,23,1,22,22,10,1,31,1,3,4,10,20,2,6,9,1,2,1,3,8,31,1,30,2,5,2,6,5,5,10,3,7,1,1,5,1,2,11,21,14,3,1,2,6,31,4,14,25,6,31,11,7,20,3,8,13,8,1,14,31)

delta=c(rep(1,9),0,rep(1,6),0,0,0,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0,1,1,0,1,0,1,0,rep(1,9),0,1,0,1,1,1,0,1,0,1,0,0,0,1,1,1,1,0,1,1,1,0,0,1,1,1,1,0)

times=sort(unique(t))

X=cbind(x1,x2,x3)

beta=rep(0.001,ncol(X))

ite=1000

eta=0.001

for(l in 1:ite){
  
dbeta=rep(0,length(beta))  

loglik=0
  
for(i in 1:nrow(X)){
  
x=X[t>=t[i],]  
  
dbeta=dbeta+delta[i]*(X[i,]-t(x)%*%exp(x%*%beta)/sum(exp(x%*%beta)))

loglik=loglik+delta[i]*(sum(X[i,]*beta)-log(sum(exp(x%*%beta))))

}  
  
beta=beta+eta*dbeta

print(loglik)
  
}

```



```{r}

#Cox部分尤度法　ニュートンラフソン p.157

x1=c(rep(0,10),1,0,1,1,1,0,0,1,0,1,1,1,0,0,0,1,1,0,1,1,1,0,0,0,0,0,1,1,0,0,0,0,0,0,1,1,0,0,0,0,1,0,0,1,0,0,1,0,0,1,0,1,1,1,0,1,0,1,1,0,1,1,0,1,0,0,0,1,0,0)

x2=c(2,3,0,3,2,1,1,1,3,0,2,0,2,0,2,3,1,2,2,1,1,2,1,3,3,1,1,2,3,1,1,3,3,3,1,3,0,1,0,2,3,2,0,2,1,2,0,3,1,1,3,3,3,2,0,3,3,2,3,3,0,3,0,2,2,2,3,1,3,0,0,3,3,3,0,2,0,1,3,0)

x3=c(1,0,1,1,1,1,0,1,0,1,0,0,1,0,1,0,1,1,1,0,1,0,1,0,0,1,1,1,0,0,1,1,1,1,1,1,1,0,1,1,1,0,0,1,0,1,1,1,0,0,0,1,0,0,1,0,1,1,1,1,1,1,0,1,0,1,1,0,1,1,0,1,0,0,1,1,0,0,0,1)

t=c(15,13,15,10,3,2,2,2,1,18,4,5,2,6,1,1,22,11,19,1,23,1,22,22,10,1,31,1,3,4,10,20,2,6,9,1,2,1,3,8,31,1,30,2,5,2,6,5,5,10,3,7,1,1,5,1,2,11,21,14,3,1,2,6,31,4,14,25,6,31,11,7,20,3,8,13,8,1,14,31)

delta=c(rep(1,9),0,rep(1,6),0,0,0,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0,1,1,0,1,0,1,0,rep(1,9),0,1,0,1,1,1,0,1,0,1,0,0,0,1,1,1,1,0,1,1,1,0,0,1,1,1,1,0)

times=sort(unique(t))

X=cbind(x1,x2,x3)

beta=rep(0.001,ncol(X))

ite=10000

eta=0.01

for(l in 1:ite){
  
dbeta=rep(0,length(beta)) 

ddbeta=array(0,dim=c(length(beta),length(beta)))

loglik=0
  
for(i in 1:nrow(X)){
  
x=X[t>=t[i],]  
  
dbeta=dbeta+delta[i]*(X[i,]-t(x)%*%exp(x%*%beta)/sum(exp(x%*%beta)))

H1=t(x)%*%diag(c(exp(x%*%beta)/sum(exp(x%*%beta))))%*%x

H2=t(x*c(exp(x%*%beta)))%*%diag(rep(1/sum(exp(x%*%beta)),nrow(x)))%*%(x*c(exp(x%*%beta)))

ddbeta=ddbeta-delta[i]*(H1-H2)

loglik=loglik+delta[i]*(sum(X[i,]*beta)-log(sum(exp(x%*%beta))))

}  
  
beta=beta+eta*solve(ddbeta)%*%dbeta

print(loglik)
  
}

```

```{r}

library(dplyr)

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56))

X=as.matrix(data[,colnames(data) %in% c("x1","x2")])

Y=data$y

beta=solve(t(X)%*%X)%*%t(X)%*%Y

e=Y-X%*%beta

sigma2=mean(e^2)

n=nrow(X);k=ncol(X)

#H0:R%*%beta=0

R=diag(1,ncol(X));r=rep(0,ncol(X))

H=X%*%solve(t(X)%*%X)%*%t(X)

M=diag(1,nrow(X))-H

s2=t(Y)%*%M%*%t(t(Y))/(n-k)

q=nrow(R)

f=t(R%*%beta-r)%*%solve(R%*%solve(t(X)%*%X)%*%t(R))%*%(R%*%beta-r)/q

f=f/s2

alpha=0.01

qf(1-alpha,q,n-k)

#誤差データの正規性の検定

m2=mean((e-mean(e))^2);m3=mean((e-mean(e))^3);m4=mean((e-mean(e))^4)

#標本歪度

m3/(m2^(3/2))

#標本尖度

m4/(m2^2)

#plot

hist(e)

JB=(m3/(m2^(3/2))+(m4/(m2^2)-3)^2/4)*n/6

qchisq(1-alpha,2)






#p.68

sigma=3;

E_uu=diag(3,n)

omega=E_uu

beta=solve(t(X)%*%solve(omega)%*%X)%*%t(X)%*%solve(omega)%*%Y

P=eigen(omega)$vector

var_beta=solve(t(X)%*%solve(omega)%*%X)


f=function(x1,x2,l){
  
if(l==1){  
  
return(exp(-(x1^2+x2^2)))  
  
}else{

return(exp(-(x1^2-x2^2)))  

}
  
}  

h=0.01;C=array(0,dim=c(k,k))

f_vec=c()

for(j in 1:k){

C[j,1]=(f(beta[1]+h,beta[2],j)-f(beta[1],beta[2],j))/h  
  
C[j,2]=(f(beta[1],beta[2]+h,j)-f(beta[1],beta[2],j))/h  

f_vec=c(f_vec,f(beta[1],beta[2],j))

}  

var=array(s2,dim=c(k,k))*C%*%solve(t(X)%*%X)%*%t(C)

library(mvtnorm)

qmvnorm(1-alpha,mean=f_vec,sigma=var)

plot_dat=data.frame(x1=seq(-1,1,0.01),x2=rep(seq(-1,1,0.01),length(seq(-1,1,0.01))),density=0)

for(j in 1:nrow(plot_dat)){

plot_dat$density[j]=dmvnorm(c(plot_dat$x1[j],plot_dat$x2[j]),mean=f_vec,sigma=var)

}

library(rgl)

plot3d(plot_dat$x1,plot_dat$x2,plot_dat$density)

```

```{r}

#ハット行列 p.74

library(dplyr)

data=data.frame(num=1:20,y=c(167,167.5,168.4,172,155.3,151.4,163,174,168,160.4,164.7,171,162.6,164.8,163.3,167.6,169.2,168,167.4,172),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56))

X=as.matrix(data[,colnames(data) %in% c("x1","x2")])

X=cbind(rep(1,nrow(X)),X)

Y=data$y

H=X%*%solve(t(X)%*%X)%*%t(X)

e=(diag(1,nrow(H))-H)%*%Y


#スチューデント化残差

#データと予測値に対して、標準化された残差

s2=sum(e^2)/(n-k)

sd_ei=sqrt(s2)*sqrt(1-diag(H))

ri=e/(sd_ei)

plot(ri)

#データの欠損におけるデータへの影響力の大きさ

X=as.matrix(data[,colnames(data) %in% c("x1","x2")])

beta=solve(t(X)%*%X)%*%t(X)%*%Y

beta=c(mean(Y)-sum(beta*apply(X,2,mean)),beta)

DFBETAS=array(0,dim=c(n,k+1))

s2_vec=(n-k)*s2/(n-k-1)-e^2/((1-diag(H))*(n-k-1))

for(j in 1:nrow(X)){
  
X_sub=X[-j,]

beta_sub=solve(t(X_sub)%*%X_sub)%*%t(X_sub)%*%Y[-j]

beta_sub=c(mean(Y[-j])-sum(apply(X_sub,2,mean)*beta_sub),beta_sub)

s2_sub=s2_vec[j]

X_sub=cbind(rep(1,n-1),X_sub)

DFBETAS[j,]=((beta-beta_sub)/(sqrt(s2_sub)*sqrt(diag(solve(t(X_sub)%*%X_sub)))))
  
}

#データに極めておおきな影響を与える大きさ

2/sqrt(n)

#p.81 スチューデント化残差

ti=e/(sqrt(s2_vec)*sqrt(1-diag(H)))

alpha=0.01

qt(1-alpha,n-k-1)

#ti^2=((n-k-1)/(n-k-ri^2))*ri^2

H_mat=H

for(j in 1:ncol(H)){
for(i in 1:nrow(H)){
  
H_mat[i,j]=e[i]*H_mat[i,j]/(sqrt(s2_vec[i])*sqrt(diag(H)[j])*(1-diag(H)[i]))  
  
}
}

DFFITS=H_mat

abs(DFFITS)>sqrt(k/(n-k))*2

#クックのD

Di=(ri^2)*(diag(H)/(1-diag(H)))/k

alpha=0.01

qf(1-alpha,1,n-k-1)/(n-k-1+qf(1-alpha,1,n-k-1))

#L-R plot (leverage,residual)

Z=cbind(X,Y)

H_Z=Z%*%solve(t(Z)%*%Z)%*%t(Z)

plot(e^2/sum(e^2),diag(H))



#p.99 

H_childer=H-array(1/n,dim=c(n,n))



```

```{r}

#bernoulli-logistic regression p.105

#multi dimension

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)

for(j in 1:nrow(data)){

data$sign[j]=ifelse(data$y[j]>0.49,1,0)

}

X=as.matrix(cbind(rep(1,nrow(data)),data[,colnames(data) %in% c("x1","x2")]))

pis=data$y

Y=as.numeric(data$sign)

alpha=rep(0.01,(ncol(X)))

ite=10

alpha_beta_data=array(0,dim=c(nrow(X),ncol(X)))

for(j in 1:ite){

V=array(0,dim=c(nrow(data),nrow(data)))

diag(V)=((1/(1+exp(-(X%*%alpha))))*(1-1/(1+exp(-(X%*%alpha)))))

print(sum(V))

mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%(pis-(1/(1+exp(-(X%*%alpha)))))

alpha=alpha+mat

alpha_beta_data[j,]=alpha

}

data=data %>% dplyr::mutate(predict=1/(1+exp(-(X%*%alpha))))

I=t(X)%*%solve(V)%*%X

#回帰係数の検定 alpha=0 ~N(0,1)

Z=alpha/sqrt(diag(solve(I)))

#confidence interval

p_hat=1/(1+exp(-(X%*%alpha)))

confidence=data.frame(num=1:nrow(data),value=0)

for(j in 1:nrow(confidence)){
  
vec=X[j,]  
  
confidence$value[j]=as.numeric(t(vec)%*%solve(I)%*%t(t(vec)))

}

p=0.01

confidence=data.frame(confidence,upper=p_hat+qnorm(1-p/2)*p_hat*(1-p_hat)*confidence$value,under=p_hat-qnorm(1-p/2)*p_hat*(1-p_hat)*confidence$value)



```



```{r}

#pois regression

library(dplyr)

Y <- c( 6,  9, 11, 11,
       21, 14, 13, 20,
       38, 36, 23, 29,
       30, 33, 46, 40,
       61, 48, 54, 44)

X <- cbind(rep(1, length(Y)),
       c(1, 1, 1, 1,
         2, 2, 2, 2,
         3, 3, 3, 3,
         4, 4, 4, 4, 
         5, 5, 5, 5)
       )

beta <- matrix(0, ncol = 2, nrow = 100)

beta[1, ] <- c(1, 10)

for(m in 2:100){
  lambda <- exp(X %*% beta[m - 1, ])
  W <- diag(lambda[, 1])
  XtWX <- t(X) %*% W %*% X
  U <- t(Y - lambda) %*% X
  beta[m, ] <- beta[m - 1, ] + solve(XtWX) %*% t(U)
}

params=beta[nrow(beta),]

predict=exp(X%*%params)

plot(c(1:length(Y)),(Y),type="p",col=2,xlab="num",ylab="values",ylim=c(min((Y)),max((Y))))

par(new=T)

plot(c(1:length(Y)),predict,type="p",col=3,xlab="num",ylab="values",ylim=c(min((Y)),max((Y))))

#尤離度

D=2*sum(Y*log(Y/predict)-(Y-predict))

qchisq(1-0.01,length(Y)-1)





```

```{r}

#冠状動脈性心疾患による死亡者数

#最小二乗法：

#死亡者数
DEATH=c(32,104,206,186,102,2,12,28,28,31)

#観測者数
PERSON=c(52407,43248,28612,12663,5317,18790,10673,5710,2585,1462)

MR=100000*DEATH/PERSON


#年齢階層(0:35歳未満、1:35~44歳、2:45~54歳、3:55～64歳、4:65~74歳、5：75~84歳
X3=c(1:5,1:5)

X4=X3^2

#喫煙X3、非喫煙０
X5=c(1:5,rep(0,5))

#喫煙1、非喫煙０
X2=c(rep(1,5),rep(0,5))

#特徴量行列
X=cbind(X2,X3,X4,X5)

#最小二乗法
beta=solve(t(X)%*%X)%*%t(X)%*%log(MR)

#予測値
pre_MR=exp(X%*%beta)

#喫煙者
exp(sum(beta))

#非喫煙者
exp(sum(beta*c(0,1,1,0)))




#ポアソン回帰:

#死亡者数
Y=DEATH

#特徴量行列
X=cbind(rep(1,length(X2)),X)

#反復計算回数
times=1000

#回帰係数の計算履歴を保存するデータを作成
beta <- matrix(0, ncol = ncol(X), nrow = times)

#初期値
beta[1, ] <- rep(1,ncol(X))

for(m in 2:times){
  
  lambda <- exp(X %*% beta[m - 1, ])
  
  W <- diag(lambda[, 1])
  
  XtWX <- t(X) %*% W %*% X
  
  D=diag(1/lambda[,1])
  
  U <- t(X)%*%W%*%D%*%(Y-lambda[,1])
  
  beta[m, ] <- beta[m - 1, ] + solve(XtWX) %*% (U)
  
}

#回帰係数
params1=beta[nrow(beta),]

#予測値
predict1=exp(X%*%params1)

#実測値のプロット
plot(c(1:length(Y)),(Y),type="p",col=2,xlab="num",ylab="values",ylim=c(min((Y)),max((Y))))

par(new=T)

#予測値のプロット
plot(c(1:length(Y)),predict1,type="p",col=3,xlab="num",ylab="values",ylim=c(min((Y)),max((Y))))

#H：回帰係数=0の統計量

Z=params1/sqrt(diag(solve(t(X)%*%W%*%X)))



X_star=sqrt(W)%*%X

H_star=X_star%*%solve(t(X_star)%*%X_star)%*%t(X_star)

d=(Y*log(Y/predict1)-(Y-predict1))

L_y_y=(Y*log(Y/mean(Y))-(Y-mean(Y)))


#尤離度決定係数
R2_dev=1-sum(d)/sum(L_y_y)


#ピアソン残差
r_pi=(Y-predict1)/sqrt(predict1)

#標準ピアソン残差
S_pi=r_pi/sqrt(1-diag(H_star))

#尤離度残差
r_di=sign(Y-predict1)*sqrt(d)

#標準尤離度残差
S_di=r_di/sqrt(1-diag(H_star))

#結果まとめ
res=data.frame(Y=Y,Y_hat=predict1,r_pi=r_pi,r_di=r_di,S_pi=S_pi,S_di=S_di)


```



```{r}



#p.137

library(dplyr)

#データ
data=data.frame(WAIS=c(4:20),n=c(2,1,2,3,2,6,6,6,2,6,7,3,4,1,1,1,1),Y=c(1,1,1,2,2,2,1,1,0,1,2,0,0,0,0,0,0)) %>% mutate(p=Y/n)

#特徴量行列
X=as.matrix(cbind(rep(1,nrow(data)),data[,colnames(data) %in% c("WAIS")]))

#各変数
pis=data$p;n=data$n;y=data$Y

#回帰係数の初期値
alpha=0.05;beta=0.05

#反復計算回数
ite=100

#回帰係数の更新歴の保存データ
alpha_beta_data=data.frame(ite=1:ite,alpha=0,beta=0)

for(j in 1:ite){

#予測確率
p=1/(1+exp(-(X%*%c(alpha,beta))))

#Wの逆行列
V=array(0,dim=c(nrow(data),nrow(data)))

diag(V)=c(n*p*(1-p))

#回帰係数の更新
mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%(y-n*p)
alpha=alpha+mat[1];beta=beta+mat[2]

print(alpha);print(beta)

alpha_beta_data$alpha[j]=alpha

alpha_beta_data$beta[j]=beta

}

#予測確率の保存
data=data %>% dplyr::mutate(predict=1/(1+exp(-(X%*%c(alpha,beta)))))

#目的変量との比較(plot)
plot(data$WAIS,data$p,type="p",col=2,xlim=c(min(data$WAIS),max(data$WAIS)),ylim=c(min(data$p),max(data$p)),xlab="WAIS",ylab="p")

par(new=T)

plot(data$WAIS,data$predict,type="p",col=3,xlim=c(min(data$WAIS),max(data$WAIS)),ylim=c(min(data$p),max(data$p)),xlab="WAIS",ylab="p")

mu=data$n*data$predict

#残差にあたる部分（尤離度）

D_y_mu1=y*log(y/mu);D_y_mu2=(n-y)*log((n-y)/(n-mu))

D_y_mu1[is.nan(D_y_mu1)>0]=0;D_y_mu2[is.nan(D_y_mu2)>0]=0

D_y_mu=2*D_y_mu1+2*D_y_mu2

#回帰にあたる部分（尤離度）

D_y_y1=y*log(y/mean(y));D_y_y2=(n-y)*log((n-y)/(n-mean(y)))

D_y_y1[is.nan(D_y_y1)>0]=0;D_y_y2[is.nan(D_y_y2)>0]=0

D_y_y=2*D_y_y1+2*D_y_y2

sj=sqrt(diag(solve(t(X)%*%V%*%(X))))

V2=V;diag(V2)=sqrt(diag(V2))

Xs=V2%*%X

Hs=Xs%*%solve(t(Xs)%*%Xs)%*%t(Xs)

#尤離度決定係数

R_dev=1-D_y_mu/D_y_y

#H0:回帰係数=0の統計量

Z=c(alpha,beta)/sj

#尤離度

d=D_y_mu

#ピアソン残差

r_pi=(y-mu)/sqrt(n*p*(1-p))

#尤離度残差

r_di=sign(y-mu)*sqrt(d)

#標準ピアソン残差

r_pi/sqrt(1-diag(Hs))

#標準尤離度残差

r_di/sqrt(1-diag(Hs))






```








```{r}

#p.192 計数データのモデル

DEATH=c(32,104,206,186,102,2,12,28,28,31)

PERSON=c(52407,43248,28612,12663,5317,18790,10673,5710,2585,1462)

MR=100000*DEATH/PERSON

X3=c(1:5,1:5)

X4=X3^2

X5=c(1:5,rep(0,5))

X2=c(rep(1,5),rep(0,5))

X=cbind(X2,X3,X4,X5)



#ポアソン回帰

Y=DEATH

beta <- matrix(0, ncol = ncol(X), nrow = 100)

beta[1, ] <- rep(1,ncol(X))

for(m in 2:100){
  lambda <- exp(X %*% beta[m - 1, ])
  W <- diag(lambda[, 1])
  XtWX <- t(X) %*% W %*% X
  U <- t(Y - lambda) %*% X
  beta[m, ] <- beta[m - 1, ] + solve(XtWX) %*% t(U)
}

XtWX1 <- t(X) %*% W %*% X

params1=beta[nrow(beta),]

predict1=exp(X%*%params1)

plot(c(1:length(Y)),(Y),type="p",col=2,xlab="num",ylab="values",ylim=c(min((Y)),max((Y))))

par(new=T)

plot(c(1:length(Y)),predict1,type="p",col=3,xlab="num",ylab="values",ylim=c(min((Y)),max((Y))))


R2_DEV=1-sum(Y*log(Y/predict1)-(Y-predict1))/sum(Y*log(Y/mean(Y))-(Y-mean(Y)))

#E[Y]=var[Y]:ポアソン分布の平均と分散が等しいことを検定 p.212

#DEAN(1992) Q1~N(0,1) Q2~N(0,1)

Q1=sum(((Y-predict1)^2-Y))/sqrt(sum(2*predict1^2))

Q2=sum(((Y-predict1)^2-Y)/predict1)/sqrt(2*nrow(X))



```

```{r}

#p.137~p.161,p.216

#負の二項回帰 model p.105

#multi dimension

library(dplyr)

DEATH=c(32,104,206,186,102,2,12,28,28,31)

PERSON=c(52407,43248,28612,12663,5317,18790,10673,5710,2585,1462)

MR=DEATH/PERSON

X3=c(1:5,1:5)

X4=X3^2

X5=c(1:5,rep(0,5))

X2=c(rep(1,5),rep(0,5))

X=cbind(rep(1,length(X3)),X2,X3,X4,X5)

pis=DEATH

alpha=rep(0.01,(ncol(X)))

ite=1000;r=4.94

eta=0.1

alpha_beta_data=array(0,dim=c(ite,ncol(X)))

for(j in 1:ite){

V=diag(c(r*exp(X%*%alpha)/(exp(X%*%alpha)+r)))

p=(1-r/(exp(X%*%alpha)+r))

D=diag(c(1/exp(X%*%alpha)))

mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%V%*%D%*%(pis-exp(X%*%alpha))

alpha=alpha+eta*mat

alpha_beta_data[j,]=alpha

}

predict=exp(X%*%alpha)

data=data.frame(X,pis,predict)

plot(c(1:nrow(data)),data$pis,type="p",col=2,xlim=c(1,nrow(data)),ylim=c(min(c(pis,predict)),max(c(pis,predict))),xlab="samples",ylab="prob",main="負の二項回帰model")

par(new=T)

plot(c(1:nrow(data)),data$predict,type="p",col=3,xlim=c(1,nrow(data)),ylim=c(min(c(pis,predict)),max(c(pis,predict))),xlab="samples",ylab="prob")

```
```{r}
#p.247

#指数分布回帰 model 

#multi dimension

library(dplyr);library(statmod)

DEATH=c(32,104,206,186,102,2,12,28,28,31)

PERSON=c(52407,43248,28612,12663,5317,18790,10673,5710,2585,1462)

MR=DEATH

X3=c(1:5,1:5)

X4=X3^2

X5=c(1:5,rep(0,5))

X2=c(rep(1,5),rep(0,5))

X=cbind(rep(1,length(X3)),X2,X3,X4,X5)

#特徴量行列
X=cbind(rep(1,length(X3)),log(X2/1000),X3)

#目的変量
pis=DEATH

#初期値
alpha=rep(0.01,(ncol(X)))

#反復計算回数と学習率
ite=10000;eta=0.01

#パラメータを保存するデータ
alpha_beta_data=array(0,dim=c(ite,ncol(X)))

for(j in 1:ite){

V=array(0,dim=c(nrow(X),nrow(X)))

diag(V)=1

D=diag(c(1/exp(X%*%alpha)))

mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%V%*%D%*%(pis-exp(X%*%alpha))

alpha=alpha+eta*mat

alpha_beta_data[j,]=alpha

loglik=sum(-X%*%alpha-pis*exp(-X%*%alpha))

print(paste0("対数尤度:",loglik))

}

#予測値
predict=exp(X%*%alpha)

data=data.frame(X,pis,predict)

plot(c(1:nrow(data)),data$pis,type="p",col=2,xlim=c(1,nrow(data)),ylim=c(min(c(pis,predict)),max(c(pis,predict))),xlab="samples",ylab="predict",main="指数回帰model")

par(new=T)

plot(c(1:nrow(data)),data$predict,type="p",col=3,xlim=c(1,nrow(data)),ylim=c(min(c(pis,predict)),max(c(pis,predict))),xlab="samples",ylab="predict")



```

```{r}

#p.137~p.161,p.216

#指数分布回帰 model p.105

#multi dimension

library(dplyr)

DEATH=c(32,104,206,186,102,2,12,28,28,31)

PERSON=c(52407,43248,28612,12663,5317,18790,10673,5710,2585,1462)

MR=DEATH

X3=c(1:5,1:5)

X4=X3^2

X5=c(1:5,rep(0,5))

X2=c(rep(1,5),rep(0,5))

X=cbind(rep(1,length(X3)),X2,X3,X4,X5)

pis=MR

alpha=rep(0.01,(ncol(X)))

ite=1000;eta=0.01

alpha_beta_data=array(0,dim=c(ite,ncol(X)))

for(j in 1:ite){

V=array(0,dim=c(nrow(X),nrow(X)))

diag(V)=1

D=diag(c(1/exp(X%*%alpha)))

mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%V%*%D%*%(pis-exp(X%*%alpha))

alpha=alpha+eta*mat

alpha_beta_data[j,]=alpha

loglik=sum(-X%*%alpha-pis*exp(-X%*%alpha))

print(paste0("対数尤度:"loglik))

}

predict=exp(X%*%alpha)

data=data.frame(X,MR,predict)

plot(c(1:nrow(data)),data$MR,type="p",col=2,xlim=c(1,nrow(data)),ylim=c(min(c(MR,predict)),max(c(MR,predict))),xlab="samples",ylab="prob",main="指数回帰model")

par(new=T)

plot(c(1:nrow(data)),data$predict,type="p",col=3,xlim=c(1,nrow(data)),ylim=c(min(c(MR,predict)),max(c(MR,predict))),xlab="samples",ylab="prob")

R2_dev=1-sum(log(predict/pis)+(pis-predict)/predict)/sum(log(mean(pis)/pis)+(pis-mean(pis))/mean(pis))


```

```{r}

#ガンマ分布回帰 model 

library(dplyr)

#死亡者数
DEATH=c(32,104,206,186,102,2,12,28,28,31)

#観測者数
PERSON=c(52407,43248,28612,12663,5317,18790,10673,5710,2585,1462)

X3=c(1:5,1:5)

X4=X3^2

X5=c(1:5,rep(0,5))

X2=c(rep(1,5),rep(0,5))

#特徴量行列
X=cbind(rep(1,length(X3)),X2,X3,X4,X5)

#目的変量
pis=DEATH

#初期値
alpha=rep(0.1,(ncol(X)))

#反復計算回数、パラメータ、学習率
ite=2000;a=1;eta=0.01

#パラメータを保存するデータ
alpha_beta_data=array(0,dim=c(ite,ncol(X)))

for(j in 1:ite){

V=array(0,dim=c(nrow(X),nrow(X)))

#尤離度
D=2*sum(log(exp(X%*%alpha)/DEATH)+(DEATH-exp(X%*%alpha))/exp(X%*%alpha))

D=D/(nrow(X))

a=(1+sqrt(1+(2*D)/3))/(2*D)

diag(V)=c(a);d=diag(c(1/exp(X%*%alpha)))

mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%V%*%d%*%(pis-exp(X%*%alpha))

alpha=alpha+eta*mat

alpha_beta_data[j,]=alpha

loglik=sum(a*(-pis/exp(X%*%alpha)-X%*%alpha)+a*log(a)-log(gamma(a))+(a-1)*log(pis))

print(loglik)

}

predict=exp(X%*%alpha)

data=data.frame(X,pis,predict)

#尤離度決定係数
R2_dev=1-sum(log(exp(X%*%alpha)/DEATH)+(DEATH-exp(X%*%alpha))/exp(X%*%alpha))/sum(log(mean(DEATH)/DEATH)+(DEATH-mean(DEATH))/mean(DEATH))

plot(c(1:nrow(data)),data$pis,type="p",col=2,xlim=c(1,nrow(data)),ylim=c(min(c(pis,predict)),max(c(pis,predict))),xlab="samples",ylab="predict",main="gamma回帰model")

par(new=T)

plot(c(1:nrow(data)),data$predict,type="p",col=3,xlim=c(1,nrow(data)),ylim=c(min(c(pis,predict)),max(c(pis,predict))),xlab="samples",ylab="predict")


```

```{r}

#逆ガウス分布回帰 model 

library(dplyr);library(statmod)

DEATH=c(674,1414,3636,338,1022,1568,170,442,1140,370,1198,3184,266,620,1070,118,332,884,292,634,2000,210,438,566,90,220,360)

X2=rep(c(-1,0,1),9)

X3=rep(c(rep(-1,3),rep(0,3),rep(1,3)),3)

X4=c(rep(-1,9),rep(0,9),rep(1,9))

#特徴量行列
X=cbind(rep(1,length(X3)),X2,X3,X4)

#目的変量
pis=DEATH

#パラメータの初期値
alpha=rep(0.001,(ncol(X)))

#反復計算回数、パラメータ、学習率
ite=10000;lambda=1;eta=0.001

#パラメータを保存するデータ
alpha_beta_data=array(0,dim=c(ite,ncol(X)))

for(j in 1:ite){

V=array(0,dim=c(nrow(X),nrow(X)))

lambda=1/mean(((pis-exp(X%*%alpha))^2)/(pis*exp(X%*%alpha)^2))

diag(V)=1/(exp(X%*%alpha)/lambda)

D=diag(c(1/exp(X%*%alpha)))

mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%V%*%D%*%(pis-exp(X%*%alpha))

alpha=alpha+eta*mat

alpha_beta_data[j,]=alpha

#対数尤度
loglik=sum(lambda*(-pis/(2*exp(X%*%alpha)^2)+1/exp(X%*%alpha))-(log(2*pi)+3*log(pis))/2-lambda/(2*pis)+log(lambda)/2)

print(loglik)

}

predict=exp(X%*%alpha)

#尤離度
d=lambda*sum((pis-exp(X%*%alpha))^2/(pis*exp(X%*%alpha)^2))

#尤離度決定係数
R2_dev=1-d/(lambda*sum((pis-mean(pis))^2/(pis*mean(pis)^2)))

data=data.frame(X,pis,predict)

plot(c(1:nrow(data)),data$pis,type="p",col=2,xlim=c(1,nrow(data)),ylim=c(min(c(pis,predict)),max(c(pis,predict))),xlab="samples",ylab="prob",main="逆ガウス回帰model")

par(new=T)

plot(c(1:nrow(data)),data$predict,type="p",col=3,xlim=c(1,nrow(data)),ylim=c(min(c(pis,predict)),max(c(pis,predict))),xlab="samples",ylab="prob")



```

```{r}

#p.247

#レイリー分布回帰 model p.105

#multi dimension

library(dplyr);library(statmod)

DEATH=c(32,104,206,186,102,2,12,28,28,31)

PERSON=c(52407,43248,28612,12663,5317,18790,10673,5710,2585,1462)

MR=DEATH/PERSON

X3=c(1:5,1:5)

X4=X3^2

X5=c(1:5,rep(0,5))

X2=c(rep(1,5),rep(0,5))

X=cbind(rep(1,length(X3)),X2,X3,X4,X5)

pis=DEATH

alpha=rep(0.01,(ncol(X)))

ite=10000;lambda=1

alpha_beta_data=array(0,dim=c(ite,ncol(X)))

cost=c()

for(j in 1:ite){
  
sita=exp(X%*%alpha)/sqrt(pi/2)

vecs=c()

for(i in 1:length(DEATH)){
  
f=function(s){
  
return(exp(-s^2/(2*sita[i]^2)+log(s)-2*log(sita[i])))  
  
}  

vecs=c(vecs,integrate(f,0,DEATH[i])$value)
  
}  


V=array(0,dim=c(nrow(X),nrow(X)))

lambda=c((4-pi)/pi)

D=diag(c(1/exp(X%*%alpha)))

diag(V)=1/lambda

print(sum(abs(pis-(1-vecs))))

cost=c(cost,sum(abs(pis-(1-vecs))))

mat=solve(t(X)%*%V%*%(X))%*%t(X)%*%V%*%D%*%(pis-PERSON*(1-vecs))

alpha=alpha+mat

alpha_beta_data[j,]=alpha

}

predict=(1-vecs)

data=data.frame(X,MR,predict)

plot(c(1:nrow(data)),data$MR,type="p",col=2,xlim=c(1,nrow(data)),ylim=c(min(c(MR,predict)),max(c(MR,predict))),xlab="samples",ylab="prob",main="レイリー分布回帰model")

par(new=T)

plot(c(1:nrow(data)),data$predict,type="p",col=3,xlim=c(1,nrow(data)),ylim=c(min(c(MR,predict)),max(c(MR,predict))),xlab="samples",ylab="prob")




```

```{r}

#p.194

#指数分布;S(t)

data=data.frame(i=c(1:21),t=c(6,6,6,6,7,9,10,10,11,13,16,17,19,20,22,23,25,32,32,34,35),delta=c(1,1,1,0,1,0,1,0,0,1,1,0,0,0,1,1,0,0,0,0,0))

p_hat=sum(data$delta)/sum(data$t)

plot(seq(1,50,1),exp(-(seq(1,50,1)*p_hat)),main="S(t):exp")

plot(seq(1,50,1),-log(exp(-(seq(1,50,1)*p_hat))),main="H(t):exp")


#weibull分布:S(t)

d=sum(data$delta)

delta=data$delta

t=data$t

times=100
newton_raphson=data.frame(times=1:times,value=0,f=0,df=0)
value=10
for(j in 1:nrow(newton_raphson)){
f<-(d/value+sum(delta*log(t))-d*sum((t^value)*log(t))/sum(t^value))

df<-(-d/(value^2)-d*(sum(value*t^(value-1))/sum(t^value)+sum(log(t)*t^value)*sum(value*t^(value-1))/(sum(t^value)^2)))

  newton_raphson$f[j]=f
  newton_raphson$df[j]=df
  value=-f/df+value

  newton_raphson$value[j]=value
}

value=tail(newton_raphson$value,1)

k=value;p=(d/sum(t^k))^(1/k)

plot(seq(1,50,1),exp(-(seq(1,50,1)/k)^p),main="S(t):weibull")

plot(seq(1,50,1),-log(exp(-(seq(1,50,1)/k)^p)),main="H(t):weibull")


#kaplan-meier estimator

data=data[order(data$t),]

t_data=data %>% group_by(t) %>% dplyr::summarise(d=sum(delta))

t_data=t_data %>% dplyr::filter(d!=0)

t_data=rbind(data.frame(t=0,d=0),t_data)

t_data=t_data %>% dplyr::mutate(r=0,c=0)

t_data=t_data %>% dplyr::mutate(t1=c(t_data$t),t2=c(t_data$t[2:nrow(t_data)],NA))



for(j in 1:nrow(t_data)){



if(j==nrow(t_data)){

data_sub=data %>% dplyr::filter(t<=max(data$t) & t>=t_data$t1[j])  

}else{

data_sub=data %>% dplyr::filter(t<t_data$t2[j] & t>=t_data$t1[j])  

}  

t_data$c[j]=length(data_sub$delta[data_sub$delta==0])

}

for(j in 1:nrow(t_data)){

t_data$r[j]=nrow(data)-sum(c(t_data$d[1:j],t_data$c[1:j]))  

}

t_data$r=c(t_data$r[1],lag(t_data$r)[2:nrow(t_data)])

t_data=t_data %>% dplyr::mutate(ratio=1-d/r) %>% dplyr::mutate(S=cumprod(ratio))

plot(t_data$t,t_data$S,type="p")

t_data=t_data %>% mutate(var_S=0,var_H=cumsum(d/(r^2)),H=cumsum(d/r),sigma_hat=0)

for(j in 1:nrow(t_data)){
  
t_data$var_S[j]=sum(t_data$d[1:j]/(t_data$r[1:j]*(t_data$r[1:j]-t_data$d[1:j])))*t_data$S[j]^2  

t_data$sigma_hat[j]=sum(t_data$d[1:j]/(t_data$r[1:j]*(t_data$r[1:j]-t_data$d[1:j])))
  
}

alpha=0.01

t_data=t_data %>% mutate(KM_upper=S+qnorm(1-alpha/2)*sqrt(var_S),KM_under=S+qnorm(alpha/2)*sqrt(var_S))

#g(x)=log(-log(x))の信頼区間

sita=exp((qnorm(1-alpha/2)*t_data$sigma_hat)/log(t_data$S))

under=t_data$S^(1/sita)

upper=t_data$S^sita

#arcsin平方変換：g(x)=asin(sqrt(x))

upper=sqrt(asin(t_data$S))+qnorm(1-alpha/2)*sqrt(t_data$var_S)/(2*sqrt(t_data$S)*sqrt(1-t_data$S))

under=sqrt(asin(t_data$S))+qnorm(alpha/2)*sqrt(t_data$var_S)/(2*sqrt(t_data$S)*sqrt(1-t_data$S))

#Hの信頼区間

t_data=t_data %>% mutate(H_under=H+qnorm(alpha/2)*sqrt(var_H),H_upper=H+qnorm(1-alpha/2)*sqrt(var_H))

#平均生存時間の推定

t_data=t_data %>% mutate(mu=0,sigma_mu=0)

for(j in 2:nrow(t_data)){
  
t_data$mu[j]=t_data$S[j-1]*(t_data$t[j]-t_data$t[j-1])  
  
}

for(j in 2:nrow(t_data)){
  
mu_vec=t_data$mu[1:j]  
  
mu_vec=cumsum(mu_vec)^2

t_data$sigma_mu[j]=sqrt(sum(mu_vec*t_data$d[1:j]/(t_data$r[1:j]*(t_data$r[1:j]-t_data$d[1:j])))  )
  
}

mu_under=t_data$mu+qnorm(alpha/2)*t_data$sigma_mu

mu_upper=t_data$mu+qnorm(1-alpha/2)*t_data$sigma_mu

```

```{r}

#p.295

#MH検定

#H0:S1(t)=S2(t)

data=data.frame(i=c(1:21),t=c(6,6,6,6,7,9,10,10,11,13,16,17,19,20,22,23,25,32,32,34,35),delta=c(1,1,1,0,1,0,1,0,0,1,1,0,0,0,1,1,0,0,0,0,0))

data=data[order(data$t),]

t_data=data %>% group_by(t) %>% dplyr::summarise(d=sum(delta))

t_data=t_data %>% dplyr::filter(d!=0)

t_data=rbind(data.frame(t=0,d=0),t_data)

t_data=t_data %>% dplyr::mutate(r=0,c=0)

t_data=t_data %>% dplyr::mutate(t1=c(t_data$t),t2=c(t_data$t[2:nrow(t_data)],NA))



for(j in 1:nrow(t_data)){



if(j==nrow(t_data)){

data_sub=data %>% dplyr::filter(t<=max(data$t) & t>=t_data$t1[j])  

}else{

data_sub=data %>% dplyr::filter(t<t_data$t2[j] & t>=t_data$t1[j])  

}  

t_data$c[j]=length(data_sub$delta[data_sub$delta==0])

}

for(j in 1:nrow(t_data)){

t_data$r[j]=nrow(data)-sum(c(t_data$d[1:j],t_data$c[1:j]))  

}

t_data$r=c(t_data$r[1],lag(t_data$r)[2:nrow(t_data)])

t_data=t_data %>% dplyr::mutate(ratio=1-d/r) %>% dplyr::mutate(S=cumprod(ratio))

t_data1=t_data

data=data.frame(i=c(1:21),t=c(5,6,6,7,7,9,10,10,11,12,16,17,20,20,22,23,25,32,32,32,35),delta=c(1,1,1,0,1,0,1,0,0,1,1,0,1,0,1,1,0,0,1,0,1))

data=data[order(data$t),]

t_data=data %>% group_by(t) %>% dplyr::summarise(d=sum(delta))

t_data=t_data %>% dplyr::filter(d!=0)

t_data=rbind(data.frame(t=0,d=0),t_data)

t_data=t_data %>% dplyr::mutate(r=0,c=0)

t_data=t_data %>% dplyr::mutate(t1=c(t_data$t),t2=c(t_data$t[2:nrow(t_data)],NA))



for(j in 1:nrow(t_data)){



if(j==nrow(t_data)){

data_sub=data %>% dplyr::filter(t<=max(data$t) & t>=t_data$t1[j])  

}else{

data_sub=data %>% dplyr::filter(t<t_data$t2[j] & t>=t_data$t1[j])  

}  

t_data$c[j]=length(data_sub$delta[data_sub$delta==0])

}

for(j in 1:nrow(t_data)){

t_data$r[j]=nrow(data)-sum(c(t_data$d[1:j],t_data$c[1:j]))  

}

t_data$r=c(t_data$r[1],lag(t_data$r)[2:nrow(t_data)])

t_data=t_data %>% dplyr::mutate(ratio=1-d/r) %>% dplyr::mutate(S=cumprod(ratio))

t_data2=t_data

ts=intersect(t_data1$t,t_data2$t)

t_data1=t_data1 %>% filter(t %in% ts)

t_data2=t_data2 %>% filter(t %in% ts)

d=t_data1$d+t_data2$d

n1=t_data1$r;n2=t_data2$r

e=d*n1/(n1+n2)

v=(d*n1*n2*(n1+n2-d))/((n1+n2-1)*(n1+n2)^2)

#kai_MH~N(0,1)

kai_MH=(sum(d)-sum(e))/sqrt(sum(v))

#wilcoxson~chisq(1)

Q_W=sum((n1+n2)*(d-e))^2/(sum(v*(n1+n2)^2))

#Tarone-Ware~chisq(1)

Q_TW=sum(sqrt(n1+n2)*(d-e))^2/(sum(v*(n1+n2)))


#Kグループの生存関数の検定

K=5;ts=20

term1=array(0,dim=c(K-1,1))

term2=array(0,dim=c(K-1,K-1))

for(j in 1:ts){
  
d=sample(1:50,K)

n=sample(10:200,K)

e=sum(d)*n/sum(n)
  
v_mat=array(0,dim=c(K-1,K-1))

for(l in 1:(K-1)){
for(m in 1:(K-1)){  
 
if(l!=m){  
  
v_mat[l,m]=sum(d)*n[l]*n[m]*(1-sum(d)/sum(n))/((sum(n)-1)*sum(n))   
  
}else{
  
v_mat[l,m]=sum(d)*n[l]*(1-n[l]/sum(n))*(1-sum(d)/sum(n))/((sum(n)-1)) 
  
}  
  
}
}

W=diag(1,K-1)

term1=term1+W%*%(d[1:(K-1)]-e[1:(K-1)])

term2=term2+(W%*%v_mat%*%W)

}

Q=t(term1)%*%solve(term2)%*%term1

qchisq(1-0.01,K-1)


```
```{r}

#比例危険度モデル　p.304

library(dplyr);library(stringr)

#weibull hazard model

data=data.frame(i=c(1:65),t=c(1,1,2,2,2,3,5,5,6,6,6,6,7,7,7,9,11,11,11,11,11,13,14,15,16,16,17,17,18,19,19,24,25,26,32,35,37,41,42,51,52,54,58,66,67,88,89,92,4,4,7,7,8,12,11,12,13,16,19,19,28,41,53,57,77),X2=c(2.218,1.94,1.519,1.748,1.301,1.544,2.236,1.681,1.362,2.114,1.114,1.415,1.978,1.041,1.176,1.724,1.114,1.23,1.301,1.508,1.079,0.778,1.398,1.602,1.342,1.322,1.23,1.591,1.447,1.079,1.255,1.301,1,1.23,1.322,1.114,1.602,1,1.146,1.568,1,1.255,1.204,1.447,1.322,1.176,1.322,1.431,1.945,1.924,1.114,1.532,1.079,1.146,1.613,1.398,1.663,1.146,1.322,1.322,1.23,1.756,1.114,1.255,1.079),X3=c(9.4,12,9.8,11.3,5.1,6.7,10.1,6.5,9,10.2,9.7,10.4,9.5,5.1,11.4,8.2,14,12,13.2,7.5,9.6,5.5,14.6,10.6,9,8.8,10,11.2,7.5,14.4,7.5,14.6,12.4,11.2,10.6,7,11,10.2,5,7.7,10.1,9,12.1,6.6,12.8,10.6,14,11,10.2,10,12.4,10.2,9.9,11.6,14,8.8,4.9,13,13,10.8,7.3,12.8,12,12.5,14),X4=c(67,38,81,75,57,46,50,74,77,70,60,67,48,61,53,55,61,43,65,70,51,60,66,70,48,62,53,68,65,51,60,56,67,49,46,48,63,69,70,74,60,49,42,59,52,47,63,58,59,49,48,81,57,46,60,66,71,55,59,69,82,72,66,66,60),X5=c(0,0,0,0,0,1,1,0,0,1,0,1,0,1,1,0,0,0,0,0,1,1,0,0,0,1,0,0,1,0,1,1,0,1,0,0,0,0,1,0,1,0,1,0,0,1,0,1,0,1,1,0,1,1,0,1,1,0,1,1,1,0,0,0,0),X6=c(10,18,15,12,9,10,9,9,8,8,10,8,10,10,13,12,10,9,10,12,9,10,10,11,10,10,9,10,8,15,9,9,10,11,9,10,9,10,9,13,10,10,10,9,10,9,9,11,10,13,10,11,8,7,9,9,9,9,10,10,9,9,11,11,12),delta=c(rep(1,48),rep(0,17)))

t=data$t;delta=data$delta

X=as.matrix(data[,str_count(colnames(data),"X")>0])

for(j in 1:ncol(X)){
  
X[,j]=(X[,j]-mean(X[,j]))/sd(X[,j])  
  
}

data=data.frame(i=data$i,t=data$t,X,delta=data$delta)

X=cbind(rep(1,nrow(X)),X)

a=1;b=rep(1,ncol(X))

ite=100

param_box=array(0,dim=c(ite,ncol(X)+1))

log_lik=c()

for(j in 1:ite){
  
mat=array(0,dim=c(ncol(X)+1,ncol(X)+1))

U=array(0,dim=c(ncol(X)+1,1))

mat[1,1]=mat[1,1]-sum(delta/(a^2))-sum((log(t)^2)*(t^a)*exp(X%*%b))

mat[2:nrow(mat),1]=mat[2:nrow(mat),1]-as.numeric(apply((X*c((t^a)*log(t)))*c(exp(X%*%b)),2,sum))

mat[1,2:ncol(mat)]=mat[1,2:ncol(mat)]-as.numeric(apply((X*c((t^a)*log(t)))*c(exp(X%*%b)),2,sum))

for(i in 1:nrow(X)){

mat[2:nrow(mat),2:ncol(mat)]=mat[2:nrow(mat),2:ncol(mat)]-t(t(X[i,]))%*%t(X[i,])*(t[i]^a)*exp(X%*%b)[i]

}

mat=-mat

U[1,1]=U[1,1]+sum(delta*((1/a)+log(t)))-sum(log(t)*(t^a)*exp(X%*%b))

U[2:nrow(U),1]=U[2:nrow(U),1]+as.numeric(apply(X*c(delta)-X*c((t^a)*exp(X%*%b)),2,sum))

params=c(a,b)

params=params+as.numeric(solve(mat)%*%U)

a=params[1];b=params[-1]

param_box[j,1]=a

param_box[j,2:ncol(param_box)]=b

print(sum(delta*(log(a)+(a-1)*log(t)+X%*%b)-(t^a)*exp(X%*%b)))

log_lik=c(log_lik,sum(delta*(log(a)+(a-1)*log(t)+X%*%b)-(t^a)*exp(X%*%b)))

}

E_T=gamma(1+1/a)*b

S=exp(-exp(X%*%b)*t^a)

S_data=data.frame(t=t,S=S,X[,2:ncol(X)]) %>% mutate(H=-log(S))

S_data=S_data[order(S_data$t),]

plot(S_data$t,S_data$S,type="o",col=2)

RC=-log(S)

RM=delta-RC

RD=sign(RM)*sqrt(-2*(RM+delta*log(RC)))

plot(RD)

X2=X*c(delta)

#加法的ハザード回帰

times=sort(unique(t))

B_box=array(0,dim=c(length(times),ncol(X)))

for(j in 1:length(times)){
  
data_sub=data %>% filter(t<=times[j])  

ts=sort(unique(data_sub$t))

B=array(0,dim=c(1,ncol(X)))

data_sub_sub=data_sub 

delta=c(data_sub_sub$delta)

X=as.matrix(data_sub_sub[,str_count(colnames(data_sub_sub),"X")>0])

X=cbind(rep(1,nrow(data_sub_sub)),X)*c(delta)

if(abs(det(t(X)%*%X))>10^(-4)){

B=solve(t(X)%*%X)%*%t(X)%*%rep(1,nrow(X))

#var_B=solve(t(X)%*%X)%*%t(X)%*%rep(1,nrow(X))%*%X%*%t(solve(t(X)%*%X))

B_box[j,]=B_box[j,]+c(B)

}
  
}  

B_box=apply(B_box,2,cumsum)





```

```{r}

survreg(Surv(time, status) ~age+sex+frail, kidney,dist="weibull")

#生存時間モデル：ワイブル分布の場合(準ニュートン法)

library(survival)

data(kidney)

delta=kidney$status

t=kidney$time

XX=as.matrix(cbind(rep(1,length(t)),kidney[,colnames(kidney) %in% c("age","sex","frail")]))

#fiacco,mccormick(1968),murtagh,sargent(1969)

#関数定義
f=function(x){

a=x[1];b=x[-1]

f1=sum(delta*(1/a+log(t))-exp(XX%*%c(b))*log(t)*t^a)

f2=apply(XX*c(delta),2,sum)-apply(XX*c(exp(XX%*%b)*t^a),2,sum)

return(c(f1,f2))

}

#初期値(1,1)
X=c(0.1,rep(0.01,ncol(XX)))

#ヘッシアンの近似行列の初期値
H=diag(1,length(X))

#反復計算回数
ite=500000

#最小点から非常にかけ離れた点になるのを防ぐためのスカラ
eta=10^(-5)

for(l in 1:ite){

#各代数方程式の値が入ったベクトル
f_val=f(X)

#座標点に対する更新量ベクトル
sigma=t(t(-H%*%f_val))

#etaによって学習幅を設定し座標点を更新
X=X+eta*sigma

#以前の座標点と更新された座標点の関数値の差のベクトル(y)
y=f(X)-f_val

#ヘッシアンの逆行列の近似行列の更新
H=H+(eta*sigma-H%*%y)%*%t(eta*sigma-H%*%y)/as.numeric(t(y)%*%(eta*sigma-H%*%y))

alpha=X[1];beta=X[-1]

#対数尤度
loglik=sum(delta*(XX%*%c(beta)+log(rep(alpha,length(t)))+(rep(alpha,length(t))-1)*log(t)))-sum(exp(XX%*%c(beta))*t^alpha)

print(loglik)

}



```

```{r}

#生存時間モデル：ワイブル分布の場合（ニュートンラフソン法）

library(survival)

data(kidney)

delta=kidney$status

t=kidney$time

XX=as.matrix(cbind(rep(1,length(t)),kidney[,colnames(kidney) %in% c("age","sex","frail")]))

ite=10^5

eta=0.01

X=c(0.1,rep(0.01,ncol(XX)))

X=rep(1,ncol(XX)+1)

for(l in 1:ite){

a=X[1];b=X[-1]

f1=sum(delta*(1/a+log(t))-exp(XX%*%c(b))*log(t)*t^a)

f2=apply(XX*c(delta),2,sum)-apply(XX*c(exp(XX%*%b)*t^a),2,sum)

df=c(f1,f2)

ddf=array(0,dim=c(length(X),length(X)))

ddf[1,1]=sum(delta*(-1/(a^2)))-sum(exp(XX%*%c(b))*(t^a)*(log(t)^2))

ddf[1,2:ncol(ddf)]=c(-apply(XX*c(exp(XX%*%c(b))*(log(t))*t^a),2,sum))

ddf[2:ncol(ddf),1]=c(-apply(XX*c(exp(XX%*%c(b))*(log(t))*t^a),2,sum))


for(i in 2:ncol(ddf)){
for(j in 2:nrow(ddf)){

ddf[i,j]=-sum(XX[,(i-1)]*XX[,(j-1)]*exp(XX%*%c(b))*t^a)

}}

X=X-eta*solve(ddf)%*%df

alpha=X[1];beta=X[-1]

#対数尤度
loglik=sum(delta*(XX%*%c(beta)+log(rep(alpha,length(t)))+(rep(alpha,length(t))-1)*log(t)))-sum(exp(XX%*%c(beta))*t^alpha)

print(loglik)

}

```

```{r}


#例9.5　骨髄腫患者の生存期間　p.320

#生存時間モデル：ワイブル分布の場合（ニュートンラフソン法）

data=data.frame(i=c(1:65),t=c(1,1,2,2,2,3,5,5,6,6,6,6,7,7,7,9,11,11,11,11,11,13,14,15,16,16,17,17,18,19,19,24,25,26,32,35,37,41,42,51,52,54,58,66,67,88,89,92,4,4,7,7,8,12,11,12,13,16,19,19,28,41,53,57,77),X2=c(2.218,1.94,1.519,1.748,1.301,1.544,2.236,1.681,1.362,2.114,1.114,1.415,1.978,1.041,1.176,1.724,1.114,1.23,1.301,1.508,1.079,0.778,1.398,1.602,1.342,1.322,1.23,1.591,1.447,1.079,1.255,1.301,1,1.23,1.322,1.114,1.602,1,1.146,1.568,1,1.255,1.204,1.447,1.322,1.176,1.322,1.431,1.945,1.924,1.114,1.532,1.079,1.146,1.613,1.398,1.663,1.146,1.322,1.322,1.23,1.756,1.114,1.255,1.079),X3=c(9.4,12,9.8,11.3,5.1,6.7,10.1,6.5,9,10.2,9.7,10.4,9.5,5.1,11.4,8.2,14,12,13.2,7.5,9.6,5.5,14.6,10.6,9,8.8,10,11.2,7.5,14.4,7.5,14.6,12.4,11.2,10.6,7,11,10.2,5,7.7,10.1,9,12.1,6.6,12.8,10.6,14,11,10.2,10,12.4,10.2,9.9,11.6,14,8.8,4.9,13,13,10.8,7.3,12.8,12,12.5,14),X4=c(67,38,81,75,57,46,50,74,77,70,60,67,48,61,53,55,61,43,65,70,51,60,66,70,48,62,53,68,65,51,60,56,67,49,46,48,63,69,70,74,60,49,42,59,52,47,63,58,59,49,48,81,57,46,60,66,71,55,59,69,82,72,66,66,60),X5=c(0,0,0,0,0,1,1,0,0,1,0,1,0,1,1,0,0,0,0,0,1,1,0,0,0,1,0,0,1,0,1,1,0,1,0,0,0,0,1,0,1,0,1,0,0,1,0,1,0,1,1,0,1,1,0,1,1,0,1,1,1,0,0,0,0),X6=c(10,18,15,12,9,10,9,9,8,8,10,8,10,10,13,12,10,9,10,12,9,10,10,11,10,10,9,10,8,15,9,9,10,11,9,10,9,10,9,13,10,10,10,9,10,9,9,11,10,13,10,11,8,7,9,9,9,9,10,10,9,9,11,11,12),delta=c(rep(1,48),rep(0,17)))

t=data$t;delta=data$delta

XX=as.matrix(cbind(rep(1,length(t)),data[,colnames(data) %in% c("X2","X3","X4","X5","X6")]))

ite=10^5

eta=0.1

X=c(0.1,rep(0.01,ncol(XX)))

X=rep(1,ncol(XX)+1)

for(l in 1:ite){

a=X[1];b=X[-1]

f1=sum(delta*(1/a+log(t))-exp(XX%*%c(b))*log(t)*t^a)

f2=apply(XX*c(delta),2,sum)-apply(XX*c(exp(XX%*%b)*t^a),2,sum)

df=c(f1,f2)

ddf=array(0,dim=c(length(X),length(X)))

ddf[1,1]=sum(delta*(-1/(a^2)))-sum(exp(XX%*%c(b))*(t^a)*(log(t)^2))

ddf[1,2:ncol(ddf)]=c(-apply(XX*c(exp(XX%*%c(b))*(log(t))*t^a),2,sum))

ddf[2:ncol(ddf),1]=c(-apply(XX*c(exp(XX%*%c(b))*(log(t))*t^a),2,sum))


for(i in 2:ncol(ddf)){
for(j in 2:nrow(ddf)){

ddf[i,j]=-sum(XX[,(i-1)]*XX[,(j-1)]*exp(XX%*%c(b))*t^a)

}}

X=X-eta*solve(ddf)%*%df

alpha=X[1];beta=X[-1]

#対数尤度
loglik=sum(delta*(XX%*%c(beta)+log(rep(alpha,length(t)))+(rep(alpha,length(t))-1)*log(t)))-sum(exp(XX%*%c(beta))*t^alpha)

print(loglik)

}

```

```{r}


#例9.5　骨髄腫患者の生存期間　

#生存時間モデル：ワイブル分布の場合(準ニュートン法)

data=data.frame(i=c(1:65),t=c(1,1,2,2,2,3,5,5,6,6,6,6,7,7,7,9,11,11,11,11,11,13,14,15,16,16,17,17,18,19,19,24,25,26,32,35,37,41,42,51,52,54,58,66,67,88,89,92,4,4,7,7,8,12,11,12,13,16,19,19,28,41,53,57,77),X2=c(2.218,1.94,1.519,1.748,1.301,1.544,2.236,1.681,1.362,2.114,1.114,1.415,1.978,1.041,1.176,1.724,1.114,1.23,1.301,1.508,1.079,0.778,1.398,1.602,1.342,1.322,1.23,1.591,1.447,1.079,1.255,1.301,1,1.23,1.322,1.114,1.602,1,1.146,1.568,1,1.255,1.204,1.447,1.322,1.176,1.322,1.431,1.945,1.924,1.114,1.532,1.079,1.146,1.613,1.398,1.663,1.146,1.322,1.322,1.23,1.756,1.114,1.255,1.079),X3=c(9.4,12,9.8,11.3,5.1,6.7,10.1,6.5,9,10.2,9.7,10.4,9.5,5.1,11.4,8.2,14,12,13.2,7.5,9.6,5.5,14.6,10.6,9,8.8,10,11.2,7.5,14.4,7.5,14.6,12.4,11.2,10.6,7,11,10.2,5,7.7,10.1,9,12.1,6.6,12.8,10.6,14,11,10.2,10,12.4,10.2,9.9,11.6,14,8.8,4.9,13,13,10.8,7.3,12.8,12,12.5,14),X4=c(67,38,81,75,57,46,50,74,77,70,60,67,48,61,53,55,61,43,65,70,51,60,66,70,48,62,53,68,65,51,60,56,67,49,46,48,63,69,70,74,60,49,42,59,52,47,63,58,59,49,48,81,57,46,60,66,71,55,59,69,82,72,66,66,60),X5=c(0,0,0,0,0,1,1,0,0,1,0,1,0,1,1,0,0,0,0,0,1,1,0,0,0,1,0,0,1,0,1,1,0,1,0,0,0,0,1,0,1,0,1,0,0,1,0,1,0,1,1,0,1,1,0,1,1,0,1,1,1,0,0,0,0),X6=c(10,18,15,12,9,10,9,9,8,8,10,8,10,10,13,12,10,9,10,12,9,10,10,11,10,10,9,10,8,15,9,9,10,11,9,10,9,10,9,13,10,10,10,9,10,9,9,11,10,13,10,11,8,7,9,9,9,9,10,10,9,9,11,11,12),delta=c(rep(1,48),rep(0,17)))

t=data$t;delta=data$delta

XX=as.matrix(cbind(rep(1,length(t)),data[,colnames(data) %in% c("X2","X3","X4","X5","X6")]))

#fiacco,mccormick(1968),murtagh,sargent(1969)

#関数定義
f=function(x){

a=x[1];b=x[-1]

f1=sum(delta*(1/a+log(t))-exp(XX%*%c(b))*log(t)*t^a)

f2=apply(XX*c(delta),2,sum)-apply(XX*c(exp(XX%*%b)*t^a),2,sum)

return(c(f1,f2))

}

#初期値(1,1)
X=c(1,rep(10^(-4),ncol(XX)))

#ヘッシアンの近似行列の初期値
H=diag(1,length(X))

#反復計算回数
ite=150000

#最小点から非常にかけ離れた点になるのを防ぐためのスカラ
eta=10^(-8)

for(l in 1:ite){

#学習率を徐々に上昇
eta=eta+10^(-9)

#各代数方程式の値が入ったベクトル
f_val=f(X)

#座標点に対する更新量ベクトル
sigma=t(t(-H%*%f_val))

#etaによって学習幅を設定し座標点を更新
X=X+eta*sigma

#以前の座標点と更新された座標点の関数値の差のベクトル(y)
y=f(X)-f_val

#ヘッシアンの逆行列の近似行列の更新
H=H+(eta*sigma-H%*%y)%*%t(eta*sigma-H%*%y)/as.numeric(t(y)%*%(eta*sigma-H%*%y))

alpha=X[1];beta=X[-1]

#対数尤度
loglik=sum(delta*(XX%*%c(beta)+log(rep(alpha,length(t)))+(rep(alpha,length(t))-1)*log(t)))-sum(exp(XX%*%c(beta))*t^alpha)

print(loglik)

}

```

```{r}

#例9.5　骨髄腫患者の生存期間　

#生存時間モデル：ワイブル分布の場合(準ニュートン法)

data=data.frame(i=c(1:65),t=c(1,1,2,2,2,3,5,5,6,6,6,6,7,7,7,9,11,11,11,11,11,13,14,15,16,16,17,17,18,19,19,24,25,26,32,35,37,41,42,51,52,54,58,66,67,88,89,92,4,4,7,7,8,12,11,12,13,16,19,19,28,41,53,57,77),X2=c(2.218,1.94,1.519,1.748,1.301,1.544,2.236,1.681,1.362,2.114,1.114,1.415,1.978,1.041,1.176,1.724,1.114,1.23,1.301,1.508,1.079,0.778,1.398,1.602,1.342,1.322,1.23,1.591,1.447,1.079,1.255,1.301,1,1.23,1.322,1.114,1.602,1,1.146,1.568,1,1.255,1.204,1.447,1.322,1.176,1.322,1.431,1.945,1.924,1.114,1.532,1.079,1.146,1.613,1.398,1.663,1.146,1.322,1.322,1.23,1.756,1.114,1.255,1.079),X3=c(9.4,12,9.8,11.3,5.1,6.7,10.1,6.5,9,10.2,9.7,10.4,9.5,5.1,11.4,8.2,14,12,13.2,7.5,9.6,5.5,14.6,10.6,9,8.8,10,11.2,7.5,14.4,7.5,14.6,12.4,11.2,10.6,7,11,10.2,5,7.7,10.1,9,12.1,6.6,12.8,10.6,14,11,10.2,10,12.4,10.2,9.9,11.6,14,8.8,4.9,13,13,10.8,7.3,12.8,12,12.5,14),X4=c(67,38,81,75,57,46,50,74,77,70,60,67,48,61,53,55,61,43,65,70,51,60,66,70,48,62,53,68,65,51,60,56,67,49,46,48,63,69,70,74,60,49,42,59,52,47,63,58,59,49,48,81,57,46,60,66,71,55,59,69,82,72,66,66,60),X5=c(0,0,0,0,0,1,1,0,0,1,0,1,0,1,1,0,0,0,0,0,1,1,0,0,0,1,0,0,1,0,1,1,0,1,0,0,0,0,1,0,1,0,1,0,0,1,0,1,0,1,1,0,1,1,0,1,1,0,1,1,1,0,0,0,0),X6=c(10,18,15,12,9,10,9,9,8,8,10,8,10,10,13,12,10,9,10,12,9,10,10,11,10,10,9,10,8,15,9,9,10,11,9,10,9,10,9,13,10,10,10,9,10,9,9,11,10,13,10,11,8,7,9,9,9,9,10,10,9,9,11,11,12),delta=c(rep(1,48),rep(0,17)))

t=data$t;delta=data$delta

XX=as.matrix(cbind(rep(1,length(t)),data[,colnames(data) %in% c("X2","X3","X4","X5","X6")]))

#Broyden Quasi-Newton Algorithm 

#関数定義
f=function(x){

a=x[1];b=x[-1]
f1=sum(delta*(1/a+log(t))-exp(XX%*%c(b))*log(t)*t^a)

f2=apply(XX*c(delta),2,sum)-apply(XX*c(exp(XX%*%b)*t^a),2,sum)

return(c(f1,f2))

}

#パラメータの初期値
X=c(1,rep(10^(-7),ncol(XX)))

#反復計算回数
ite=10^6

#ヤコビアンの初期値
H=diag(f(X))

#学習率
eta=10^(-11)

for(l in 1:ite){

eta=eta+10^(-11)

#以前の座標点を保存
X_pre=X  

#座標点を更新
X=X-eta*H%*%f(X)  

#以前の座標点と更新された座標点の差のベクトル(s)
s=X-X_pre

#以前の座標点と更新された座標点の関数値の差のベクトル(y)
y=f(X)-f(X_pre)

#ヤコビアンの近似行列を更新する
H=H+((s-H%*%y)/as.numeric(t(s)%*%H%*%y))%*%t(s)%*%H

alpha=X[1];beta=X[-1]

#対数尤度
loglik=sum(delta*(XX%*%c(beta)+log(rep(alpha,length(t)))+(rep(alpha,length(t))-1)*log(t)))-sum(exp(XX%*%c(beta))*t^alpha)

print(loglik)

}
```
```{r}

#加速故障時間モデル

#指数分布 p.336

t=c(19,25,30,34,37,46,47,51,56,57,61,66,67,74,78,86,122,123,130,130,133,134,136,141,143,148,151,152,153,154,156,162,164,165,182,189,22,23,38,42,73,77,89,115,144)

IMM=c(rep(1,36),rep(0,9))

X=cbind(rep(1,length(IMM)),IMM)

delta=c(rep(1,16),rep(0,20),rep(1,8),0)

ite=10000

eta=0.001

beta=rep(0.01,ncol(X))

for(l in 1:ite){
  
dbeta=apply(t(X*c(delta)),1,sum)-apply(X*c(exp(X%*%beta)*t),2,sum)

beta=beta+eta*dbeta

loglik=sum(delta*X%*%beta)-sum(exp(X%*%beta)*t)

print(loglik)
 
}

```


```{r}

#対数正規分布 p.341

t=c(19,25,30,34,37,46,47,51,56,57,61,66,67,74,78,86,122,123,130,130,133,134,136,141,143,148,151,152,153,154,156,162,164,165,182,189,22,23,38,42,73,77,89,115,144)

IMM=c(rep(1,36),rep(0,9))

X=cbind(rep(1,length(IMM)),IMM)

delta=c(rep(1,16),rep(0,20),rep(1,8),0)

ite=10000

eta=0.001

beta=rep(0.01,ncol(X))

sig=1

h=0.001

for(l in 1:ite){
  
vec=beta

db1=c();db2=c()

for(j in 1:ncol(X)){
  
vec_sub=vec;vec_sub[j]=vec_sub[j]+h  
  
db1=c(db1,sum(delta*(log(dnorm((log(t)-X%*%vec_sub)/sig))-log(dnorm((log(t)-X%*%vec)/sig)))/h))

db2=c(db2,sum((delta-1)*(log(1-pnorm((log(t)-X%*%vec_sub)/sig))-log(1-pnorm((log(t)-X%*%vec)/sig)))/h))

}

dbeta=db1-db2

dsig1=sum(delta*((log(dnorm((log(t)-X%*%vec)/(sig+h)))-log(dnorm((log(t)-X%*%vec)/sig)))/h-1/sig))

dsig2=sum((delta-1)*(log(1-pnorm((log(t)-X%*%vec)/(sig+h)))-log(1-pnorm((log(t)-X%*%vec)/sig)))/h)

dsig=dsig1-dsig2

beta=beta+eta*dbeta

sig=sig+eta*dsig
  
z=(log(t)-X%*%beta)/sig

loglik=sum(delta*(log(dnorm(z))-log(sig)-log(t))-(delta-1)*log(1-pnorm(z)))

print(loglik)
 
}

```


```{r}

#対数ロジスティック分布 p.343 確認中

t=c(19,25,30,34,37,46,47,51,56,57,61,66,67,74,78,86,122,123,130,130,133,134,136,141,143,148,151,152,153,154,156,162,164,165,182,189,22,23,38,42,73,77,89,115,144)

IMM=c(rep(1,36),rep(0,9))

X=cbind(rep(1,length(IMM)),IMM)

delta=c(rep(1,16),rep(0,20),rep(1,8),0)

ite=10000

eta=0.001

del=rep(-10,ncol(X))

sig=10

for(l in 1:ite){
  
dsig=sum(delta*(-1/sig-(log(t)-X%*%del)/(sig^2)))+sum((delta+1)*(log(t)-X%*%del)*exp((log(t)-X%*%del)/sig)/((1+exp((log(t)-X%*%del)/sig))*sig^2))  
  
ddel=apply((-t(X*c(delta))*c(del))/sig,1,sum)+apply(X*c((delta+1)*exp((log(t)-X%*%del)/sig)/(sig*(1+exp((log(t)-X%*%del)/sig)))),2,sum)

sig=sig+eta*dsig

del=del+eta*ddel

loglik=sum(delta*(-log(sig)+(log(t)-X%*%del)/sig))-sum((delta+1)*log(1+exp((log(t)-X%*%del)/sig)))

print(loglik)
 
}

```

```{r}

#サービスサイエンスの事訳　ーデータサイエンスと数理科学の融合ー　（筑波大学出版）p.175

#Page Rank

#ぶら下がりページが存在しない
A=matrix(c(0,1,rep(0,4),1,0,1,1,0,0,1,1,0,0,0,1,1,rep(0,5),1,0,0,1,0,0,0,0,0,1,1,0),nrow=6)

#ぶら下がりページが存在
A=matrix(c(rep(0,6),1,0,1,1,0,0,1,1,0,0,0,1,rep(0,6),rep(c(1,0,0),2),rep(0,3),1,1,0),nrow=6)

H=t(t(A)/apply(A,2,sum))
H[is.nan(H)>0]=0

a=ifelse(apply(A,2,sum)>0,0,1)
S=H+t(t(rep(1,nrow(H))))%*%a/nrow(H)

n=100
P=rep(1/ncol(S),ncol(S))
for(i in 1:n){
P=S%*%P  
}


eig_values=eigen(S)$values
eigen_vector=Re(eigen(S)$vectors[,1])
eigen_vector/sum(eigen_vector)



```

```{r}

#サービスサイエンスの事訳　ーデータサイエンスと数理科学の融合ー　（筑波大学出版）p.175

#HITS

#ぶら下がりページが存在しない
A=matrix(c(0,1,rep(0,4),1,0,1,1,0,0,1,1,0,0,0,1,1,rep(0,5),1,0,0,1,0,0,0,0,0,1,1,0),nrow=6)

#ぶら下がりページが存在
A=matrix(c(rep(0,6),1,0,1,1,0,0,1,1,0,0,0,1,rep(0,6),rep(c(1,0,0),2),rep(0,3),1,1,0),nrow=6)

H=t(t(A)/apply(A,2,sum))
H[is.nan(H)>0]=0

a=ifelse(apply(A,2,sum)>0,0,1)
S=H+t(t(rep(1,nrow(H))))%*%a/nrow(H)

n=100
x=rep(1/ncol(H),ncol(H))
for(i in 1:n){
x=A%*%t(A)%*%x
x=x/sum(x)
}
y=t(A)%*%x
y=y/sum(y)


```

```{r}

#サービスサイエンスの事訳　ーデータサイエンスと数理科学の融合ー　（筑波大学出版）p.175

#SALSA

A=matrix(c(rep(0,6),1,0,1,1,0,0,1,1,0,0,0,1,rep(0,6),rep(c(1,0,0),2),rep(0,3),1,1,0),nrow=6)

Ar=t(t(A)/apply(A,2,sum))
Ar[is.nan(Ar)>0]=0
Ac=((A)/apply(A,1,sum))

#権威行列
A_childer=Ar%*%t(Ac)
A_childer_vector=Re(eigen(A_childer)$vector[,1])
A_childer_vector=A_childer_vector/sum(A_childer_vector)

#ハブ行列
H_childer=t(Ac)%*%Ar
H_childer_vector=Re(eigen(H_childer)$vector[,1])
H_childer_vector=H_childer_vector/sum(H_childer_vector)


```