
#p.87 多項反応モデル

require(MASS)
data(crabs)

df=data.frame(crabs)

sp=unique(df$sp)
sex=unique(df$sex)
df$BD=round(df$BD)

X=as.matrix(cbind(1,df[,colnames(df) %in% c("FL","RW","CL","CW")]))
Yij=t(matrix(c(df$BD),nrow=50))

ite=10
beta=rep(0.05,ncol(X))
eta=10^(-0)

for(l in 1:ite){
dbeta=apply(t(X*c(Yij)),1,sum)-apply(X*c(exp(X%*%beta)),2,sum)
ddbeta=-t(X)%*%diag(c(exp(X%*%beta)))%*%X
beta=beta-eta*solve(ddbeta)%*%dbeta
loglik=sum(c(Yij)*c(X%*%beta))-sum(exp(X%*%beta))
#print(loglik)
dev=2*sum(c(Yij)*log(c(Yij)/exp(X%*%beta))-(c(Yij)-exp(X%*%beta)))
print(dev)
}

exp_val=matrix(exp(X%*%beta),nrow=50)
pi_ij=t(exp_val)/apply(exp_val,2,sum)

#p.115 ゼロ過剰ポアソン分布

n=100
var=sample(1:5,5)
Y=sample(0:5,n,replace=T,prob=c(0.4,0.6*var/sum(var)))

x=rep(0,length(Y))
x[Y>0]=sample(3:5,length(x[Y>0]),replace=T)
x[Y<2]=sample(1:2,length(x[Y<2]),replace=T)
z=ifelse((Y>0)&(x>0),1,0)
x=cbind(1,x);z=cbind(1,z)

loglik=function(b,r){
pis=1/(1+exp(-z%*%r));lam=exp(x%*%b)
return(sum(log(pis+(1-pis)*exp(-lam))[Y==0])+sum((log(1-pis)+Y*log(lam)-lam-log(factorial(Y)))[Y>0]))
}

beta=sample(1:5,ncol(x))/10
rs=sample(1:5,ncol(z))/10

ite=10^4
h=0.01
eta=10^(-5)
for(l in 1:ite){
for(j in 1:length(beta)){
betasub=beta;betasub[j]=betasub[j]+h  
beta=beta+eta*(loglik(betasub,rs)-loglik(beta,rs))/h
}
for(j in 1:length(rs)){
rssub=rs;rssub[j]=rssub[j]+h
rs=rs+eta*(loglik(beta,rssub)-loglik(beta,rs))/h
}  
print(loglik(beta,rs))    
}

#p.147

#メトロポリスヘイスティング法
#ベイズ線形回帰
#事前分布：beta~N(0,tau^2)
#提案分布：beta~N(beta,tau^2) 

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))


Y=data$y
X=cbind(1,data$x1,data$x2)

posterior=function(b,sig,Tau){
return(sum(-b^2/(2*Tau))-sum((Y-X%*%b)^2)/(2*sig^2)-log(sqrt(2*pi*sig^2))-log(sqrt(2*pi*Tau)))  
}


beta=solve(t(X)%*%X)%*%t(X)%*%Y
tau=1
sigma=sqrt(sum((Y-X%*%beta)^2)/length(Y))


ite=10^4
eta=0.1
for(l in 1:ite){

for(j in 1:length(beta)){  
betasub=beta;betasub[j]=rnorm(1,beta[j],sqrt(tau))
r=min(c(1,exp(posterior(betasub,sigma,tau)-posterior(beta,sigma,tau))))
p=sample(0:1,1,prob=c(1-r,r))
if(p==1){
beta=beta*(1-eta)+betasub*eta  
}
}
  
sigma=sqrt(sum((Y-X%*%beta)^2)/length(Y))
print(posterior(beta,sigma,tau))  
    
}


#p.164 ベイズポアソン回帰
#メトロポリスヘイスティング法
#事前分布：beta[j]~N(0,1)
#提案分布：beta[j]~N(beta[j],1)

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))


Y=round(data$y)
X=cbind(1,data$x1,data$x2)
fact_Y=0
for(i in 1:length(Y)){
fact_Y=fact_Y+sum(log(c(1:Y[i])))  
}

posterior=function(b){
mu=exp(X%*%b)
return(sum(Y*log(mu)-mu)-fact_Y+sum(log(dnorm(b,0,1))))  
}

beta=rep(0.01,ncol(X))
ite=10^4
eta=0.01
for(l in 1:ite){

for(j in 1:length(beta)){  
betasub=beta;betasub[j]=rnorm(1,beta[j],1)
r=min(c(1,exp(posterior(betasub)-posterior(beta))))
p=sample(0:1,1,prob=c(1-r,r))
if(p==1){
beta=beta*(1-eta)+betasub*eta  
}
}
  
print(posterior(beta))
}