母集団薬物データの解析

```{r}

#p.43

data("Indometh")


ite=10000;eta=0.001

Subjects=sort(unique(Indometh$Subject))

K_el=rep(0.01,length(Subjects))

D=rep(1,length(Subjects))

V=rep(1,length(Subjects))


for(l in 1:ite){
  
cost=c()  
  
for(j in 1:length(Subjects)){  
  
dat=Indometh %>% filter(Subject==j)  
 
dK=2*sum((D[j]*dat$time/V[j])*exp(-K_el[j]*dat$time)*(dat$conc-D[j]*exp(-K_el[j]*dat$time)/V[j]))

dV=2*sum((D[j]/(V[j]^2))*exp(-K_el[j]*dat$time)*(dat$conc-D[j]*exp(-K_el[j]*dat$time)/V[j]))

K_el[j]=K_el[j]-eta*dK 

V[j]=V[j]-eta*dV

val=sum((dat$conc-(D[j]/V[j])*exp(-K_el[j]*dat$time))^2)

cost=c(cost,val)
  
}
  
print(sum(cost))  

}

```
```{r}

#p.59

library(cubature)

data("Indometh")

sig=1;sig_v=1;sig_k=1;mu_k_val=1;mu_v_val=1

Subjects=sort(unique(Indometh$Subject))

cost=function(sigma,sigma_v,sigma_k,mu_k,mu_v){

val=c()

for(j in 1:length(Subjects)){
  
dat=Indometh %>% filter(Subject==j)  

log_f=function(x){
  
value=prod(dnorm(dat$conc,D[j]*exp(-x[1]*dat$time)/x[2],sigma))*dnorm(x[1],mu_k,sigma_k)*dnorm(x[2],mu_v,sigma_v)  

return(value)
  
}    

monte=c()

h=0.1

k_interval=seq(0,20,h);v_interval=seq(h,5,h)

for(i in 1:length(k_interval)){
for(j in 1:length(v_interval)){  
  
vec=c(k_interval[i],v_interval[j])  
  
monte=c(monte,log_f(vec)) 
  
}}
  
integral=log(sum(monte[is.na(monte)==0]))

val=c(val,integral)
  
}
  
return((sum(val)))  
  
}



ite=100;eta=0.0001

h=0.01

for(l in 1:ite){
  
values=cost(sig,sig_v,sig_k,mu_k_val,mu_v_val)  
  
sig=sig+eta*(cost(sig+h,sig_v,sig_k,mu_k_val,mu_v_val)-cost(sig,sig_v,sig_k,mu_k_val,mu_v_val))/h  

sig_v=sig_v+eta*(cost(sig,sig_v+h,sig_k,mu_k_val,mu_v_val)-cost(sig,sig_v,sig_k,mu_k_val,mu_v_val))/h  

sig_k=sig_k+eta*(cost(sig,sig_v,sig_k+h,mu_k_val,mu_v_val)-cost(sig,sig_v,sig_k,mu_k_val,mu_v_val))/h  

mu_k_val=mu_k_val+eta*(cost(sig,sig_v,sig_k,mu_k_val+h,mu_v_val)-cost(sig,sig_v,sig_k,mu_k_val,mu_v_val))/h  

mu_v_val=mu_v_val+eta*(cost(sig,sig_v,sig_k,mu_k_val,mu_v_val+h)-cost(sig,sig_v,sig_k,mu_k_val,mu_v_val))/h  

print(cost(sig,sig_v,sig_k,mu_k_val,mu_v_val))
  
}

```
```{r}

#p.71

data("Indometh")

Indometh$conc=Indometh$conc

Subjects=sort(unique(Indometh$Subject))

K_el=rep(0.01,length(Subjects))

D=rep(1,length(Subjects))

V=rep(1,length(Subjects))

cost=function(k_el,v){
  
val=c()

for(j in 1:length(Subjects)){
  
dat=Indometh %>% filter(Subject==j)  

var=exp(-(dat$conc-D[j]*exp(-k_el[j]*dat$time)/v[j])^2)
  
val=c(val,sum(-log(var)-(dat$conc-D[j]*exp(-k_el[j]*dat$time)/v[j])^2/(var)))  
  
}

return(sum(val))
  
}




ite=100000

eta=0.001

h=0.01

for(l in 1:ite){
  
value=1
  
vec=K_el

for(j in 1:length(K_el)){
  
vec_sub=vec;vec_sub[j]=vec_sub[j]+h

K_el[j]=K_el[j]+eta*value*(cost(vec_sub,V)-cost(vec,V))/h
  
}

vec=V

for(j in 1:length(V)){
  
vec_sub=vec;vec_sub[j]=vec_sub[j]+h  
  
V[j]=V[j]+eta*value*(cost(K_el,vec_sub)-cost(K_el,vec))/h

}
  

print(cost(K_el,V))  

}



```

```{r}


#p.102

N=10000;n=1000;B=200

x=rnorm(N)

pis=c()

for(j in 1:B){
  
x_sub=sample(x,n,replace=T)

pis=c(pis,mean(x_sub))
  
}

se=sum((pis-mean(pis))^2)/(B-1)

p=mean(x)

p-se*qnorm(1-0.05/2)

p+se*qnorm(1-0.05/2)

```


```{r}

#p.95

data("Indometh")

t=mean(Indometh$time)

time=Indometh$time;x=Indometh$conc

seg=sort(as.numeric(unique(Indometh$Subject)))

sig=1;k0=10

log_lik=function(X,t_vec,k_vec,sigma){
  
k10=k_vec[1];k12=k_vec[2];k21=k_vec[3];V1=k_vec[4]  
  
a=(k10+k12+k21+sqrt((k10+k12+k21)^2-4*k21*k10))/2

b=(k10+k12+k21-sqrt((k10+k12+k21)^2-4*k21*k10))/2

val=c()

for(s in 1:length(t_vec)){
  
if(t_vec[s]<=t){
  
value=as.numeric((k0*(k21-a)*(exp(-a*t_vec[s])-1)/(V1*a*(a-b))+k0*(b-k21)*(exp(-b*t_vec[s])-1)/(V1*b*(a-b))))  
  
val=c(val,value) 
  
}else{
  
value=as.numeric((k0*(k21-a)*(exp(-a*t_vec[s])-1)*exp(-a*(t_vec[s]-t)))/(V1*b*(a-b))+k0*(b-k21)*(exp(-b*t_vec[s])-1)*exp(-b*(t_vec[s]-t))/(V1*(a-b)*b))  
  
val=c(val,value)
  
}  
  
}

return(sum(-(log(X)-log(val))^2/(2*sigma^2)-log(2*pi*(X*sigma)^2)))

}




ite=10000

eta=0.0001

h=0.01;sig=1

ks=array(sample(1:10,4*length(seg),replace=T),dim=c(4,length(seg)))

for(l in 1:ite){
  
cost=c()  

dsig=c()
  
for(j in seg){

dat_sub=Indometh %>% filter(Subject==j)

x_sub=dat_sub$conc;t_sub=dat_sub$time

vec=ks[,j]

for(i in 1:length(vec)){
    
vec_sub=vec;vec_sub[i]=vec_sub[i]+h  
  
ks[i,j]=ks[i,j]+eta*(log_lik(x_sub,t_sub,vec_sub,sig)-log_lik(x_sub,t_sub,vec,sig))/h  
  
}

dsig=c(dsig,(log_lik(x_sub,t_sub,vec,sig+h)-log_lik(x_sub,t_sub,vec,sig))/h)
  
cost=c(cost,log_lik(x_sub,t_sub,vec,sig))

}

sig=sig+eta*sum(dsig)

print(sum(cost))

}

```