graphical2

3.c

theta4 = -1
theta1 = theta2 = theta3 = theta4 

theta34 = -0.5
theta12 = theta23 = theta14 = theta34

x1=rpois(n=10000,lambda = exp(theta1))
x2=sapply(x1,FUN = function(x1) rpois(n=1,lambda = exp(theta2 + theta12*x1)))
x3=sapply(x2,FUN = function(x2) rpois(n=1,lambda = exp(theta3 + theta23*x2)))
x = data.frame(x1=x1,x2=x2,x3=x3)
x['x4'] = apply(x,1,FUN = function(x) rpois(n=1,lambda = exp(theta4 + theta14*x[1]+theta34*x[3])))

obj1<-function(t1) {
  sum(exp(t1) - t1*x1)
}

optimize(obj1,c(-2,2))$minimum

## [1] -1.0105864535973974

obj2 <- function(t){
  t2<-t[1]
  t12 <- t[2]
  sum(exp(t2+t12*x1)-t2*x2-t12*x1*x2)
}

obj2.g <- function(t){
  t2<-t[1]
  t12 <- t[2]
  c(sum(exp(t2+t12*x1)-x2)
    ,sum(exp(t2+t12*x1)*x1-x1*x2))
}
optim(c(-1,-0.5),obj2,obj2.g)$par

## [1] -1.01770500019192678 -0.48976898379623895

obj3 <- function(t){
  t3<-t[1]
  t23 <- t[2]
  sum(exp(t3+t23*x2)-t3*x3-t23*x2*x3)
}
obj3.g <- function(t){
  t3<-t[1]
  t23 <- t[2]
  c(sum(exp(t3+t23*x2)-x3)
    ,sum(exp(t3+t23*x2)*x2-x2*x3))
}
optim(c(-1,-0.5),obj3,obj3.g)$par

## [1] -1.01279158927500257 -0.42055039014667273

obj4 <- function(t){
  t4=t[1]
  t14=t[2]
  t34=t[3]
  x4=x$x4
  sum(exp(t4+t14*x1+t34*x3)-t4*x4-t14*x1*x4-t34*x3*x4)
}

obj4.g <- function(t){
  t4=t[1]
  t14=t[2]
  t34=t[3]
  x4=x$x4
  c(sum(exp(t4+t14*x1+t34*x3)-x4),
    sum(exp(t4+t14*x1+t34*x3)*x1-x1*x4),
    sum(exp(t4+t14*x1+t34*x3)*x3-x3*x4))
}

optim(c(-1,-0.5,-0.5),obj4,obj4.g)$par

## [1] -0.98873528049536996 -0.49892500751476943 -0.48072774122346995

5.d

count= c(293,4,176,3,23,2,197,17)
c1=c(0,0,0,0,1,1,1,1)
x1<-unlist(sapply(1:8,FUN = function(x) rep(c1[x],count[x])))
c2=c(1,1,0,0,1,1,0,0)
x2<-unlist(sapply(1:8,FUN = function(x) rep(c2[x],count[x])))
c3=c(1,0,1,0,1,0,1,0)
x3<-unlist(sapply(1:8,FUN = function(x) rep(c3[x],count[x])))

a1 <- function(t){
  t1=t[1]
  t2=t[2]
  t3=t[3]
  t12=t[4]
  t13=t[5]
  (exp(t1)+exp(t1+t2+t12)+exp(t1+t3+t13)+exp(t1+t2+t3+t12+t13))*
    (exp(0)+exp(t1)+exp(t2)+exp(t3)+exp(t1+t2+t12)+exp(t1+t3+t13)
     +exp(t2+t3)+exp(t1+t2+t3+t12+t13))^(-1)
  
}
a2 <- function(t){
  t1=t[1]
  t2=t[2]
  t3=t[3]
  t12=t[4]
  t13=t[5]
  (exp(t2)+exp(t1+t2+t12)+exp(t2+t3)+exp(t1+t2+t3+t12+t13))*
    (exp(0)+exp(t1)+exp(t2)+exp(t3)+exp(t1+t2+t12)+exp(t1+t3+t13)
     +exp(t2+t3)+exp(t1+t2+t3+t12+t13))^(-1)
}

a3 <- function(t){
  t1=t[1]
  t2=t[2]
  t3=t[3]
  t12=t[4]
  t13=t[5]
  (exp(t3)+exp(t1+t3+t13)+exp(t2+t3)+exp(t1+t2+t3+t12+t13))*(exp(0)+exp(t1)+exp(t2)+exp(t3)+
                                                                        exp(t1+t2+t12)+exp(t1+t3+t13)+exp(t2+t3)+
                                                                        exp(t1+t2+t3+t12+t13))^(-1)
}

a12 <- function(t){
  t1=t[1]
  t2=t[2]
  t3=t[3]
  t12=t[4]
  t13=t[5]
  (exp(t1+t2+t12)+exp(t1+t2+t3+t12+t13))*
    (exp(0)+exp(t1)+exp(t2)+exp(t3)+exp(t1+t2+t12)+exp(t1+t3+t13)+exp(t2+t3)+exp(t1+t2+t3+t12+t13))^(-1)
}

a13 <- function(t){
  t1=t[1]
  t2=t[2]
  t3=t[3]
  t12=t[4]
  t13=t[5]
  (exp(t1+t3+t13)+exp(t1+t2+t3+t12+t13))*
    (exp(0)+exp(t1)+exp(t2)+exp(t3)+exp(t1+t2+t12)+exp(t1+t3+t13)+exp(t2+t3)+exp(t1+t2+t3+t12+t13))^(-1)
}

aprime<-function(t){
  return(c(a1(t),a2(t),a3(t),a12(t),a13(t)))
}


initial=c(-1,1,-1,1,1)
step=0.01
container<-seq(5)
for(i in seq(10000)){
  container[1]<- initial[1] - step*(aprime(initial)[1]-mean(x1))
  container[2]<- initial[2] - step*(aprime(initial)[2]-mean(x2))
  container[3]<- initial[3] - step*(aprime(initial)[3]-mean(x3))
  container[4]<- initial[4] - step*(aprime(initial)[4]-mean(x1*x2))
  container[5]<- initial[5] - step*(aprime(initial)[5]-mean(x1*x3))
  initial <- container
  }

f123<-function(x,t){
  t1=t[1]
  t2=t[2]
  t3=t[3]
  t12=t[4]
  t13=t[5]
  x1=x[1]
  x2=x[2]
  x3=x[3]

  temp=t1*x1+t2*x2+t3*x3+t12*x1*x2+t13*x1*x3
  exp(temp-log(exp(0)+exp(t1)+exp(t2)+exp(t3)+
                  exp(t1+t2+t12)+exp(t1+t3+t13)+exp(t2+t3)+
                  exp(t1+t2+t3+t12+t13)))
}

X<-matrix(c(0,0,0,1,0,0,0,1,0,0,0,1,0,1,1,1,0,1,1,1,0,1,1,1)
          ,nrow = 3)
nor<-sum(apply(X,2,f123,t=initial))#nor=1

xxx<-data.frame(x1=c1,x2=c2,x3=c3)
apply(xxx,1,f123,t=initial)*715

## [1] 280.52803012227940371  13.11316722511849875 180.46872511348232138
## [4]   8.43593622457440162  26.96391610360172209   0.77058370993059899
## [7] 199.03164922049126062   5.68799228052209038

6.d

f1 <- function(x) {
  exp(sum(x) + 0.5*(x[1]*x[2]+x[1]*x[3]+x[1]*x[4]))
  }

f2 <- function(x) {
  exp(sum(x) + 0.5*(x[1]*x[2]+x[2]*x[3]+x[3]*x[4]+x[1]*x[4]))
}

l<-sapply(1:4,function(x) combn(1:4,x))#0
ll<-sapply(l,function(x) {
  apply(x,2,function(i) {
    t=rep(0,4)
    t[i]<-1
    t})
})
X<-matrix(c(rep(0,4),unlist(ll)),nrow = 4)
X

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,]    0    1    0    0    0    1    1    1    0     0     0     1     1
## [2,]    0    0    1    0    0    1    0    0    1     1     0     1     1
## [3,]    0    0    0    1    0    0    1    0    1     0     1     1     0
## [4,]    0    0    0    0    1    0    0    1    0     1     1     0     1
##      [,14] [,15] [,16]
## [1,]     1     0     1
## [2,]     0     1     1
## [3,]     1     1     1
## [4,]     1     1     1

G1

iteration=100000
u<-runif(n=iteration)
burn<-100

candi.s<-sample(1:16,size = iteration,replace = TRUE)
candidate<-X[,candi.s]
out2=out<-candidate#out sample
out2.s=out.s<-candi.s

for(i in seq(iteration-1)){
  if(u[i]<min(1,f1(candidate[,i])/f1(out[,i]))){
    out[,i+1]<-candidate[,i]
    out.s[i+1]<-candi.s[i]
  }else{
    out[,i+1]<-out[,i]
    out.s[i+1]<-out.s[i]
  }
}

count<-table(out.s[burn:iteration])
sum(count[which(X[1,]==1)])/sum(count)

## [1] 0.89983083252419893

sum(count[which(X[2,]==1)])/sum(count)

## [1] 0.81095284331488171

sum(count[which(X[3,]==1)])/sum(count)

## [1] 0.81068257574999247

sum(count[which(X[4,]==1)])/sum(count)

## [1] 0.81113302169147461

G2

for(i in seq(iteration-1)){
  if(u[i]<min(1,f2(candidate[,i])/f2(out[,i]))){
    out2[,i+1]<-candidate[,i]
    out2.s[i+1]<-candi.s[i]
  }else{
    out2[,i+1]<-out2[,i]
    out2.s[i+1]<-out2.s[i]
  }
}

count<-table(out2.s[burn:iteration])
sum(count[which(X[1,]==1)])/sum(count)

## [1] 0.83466631965645988

sum(count[which(X[2,]==1)])/sum(count)

## [1] 0.83866027367093421

sum(count[which(X[3,]==1)])/sum(count)

## [1] 0.83967127456181623

sum(count[which(X[4,]==1)])/sum(count)

## [1] 0.83801963944304858