rm(list = ls())
setwd("c:\\E-drive\\Books\\LOGLIN3\\BAYES\\")
load("post-samp.Rda")
# enter data
n=23
y=c(1,1,1,1,0,0,0,0,0,0,0,0,1,1,0,0,0,1,0,0,0,0,0)
N=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
tau=c(53,57,58,63,66,67,67,67,68,69,70,70,70,70,72,
73,75,75,76,76,78,79,81)
iterates=length(beta[1,])

# define likelihood contribution function
LC <- function(beta,x,N,y)
{
LP = sum(beta*x)
lc=((exp(LP)/(1+exp(LP)))**y)*((1-(exp(LP)/(1+exp(LP))))**(N-y))
return(lc)}

# Find s dim. vector Pr(y=1|Y,x_i^f)
# Xf is predictive model matrix
Xf=matrix(c(rep(1,11),seq(31,51,2)),11)
JJf=exp(Xf %*% beta)
JJf=JJf/(1+JJf)
PPf=rowMeans(JJf)

# Find n vector Pr(y=1|Y,x_i)
# X is model matrix
X=matrix(c(rep(1,n),tau),nrow=n)
JJ=exp(X %*% beta)
JJ=JJ/(1+JJ)
PP=rowMeans(JJ)

# Define dummy vectors for diagnostics
# to be redefined below
D1=rep(1,n)
D2=rep(1,n)
D=rep(1,n)
Dp=rep(1,n)
Df=rep(1,n)
lc=rep(1,iterates)
qti=rep(1,iterates)

# Diagnostics
# outer "for" loop finds diagnostics for each i
for(i in 1:n)
{
# inner "for" loop creates "iterates" length vectors of
# L(\beta^r|y_i) and unstandardized \qtilde_{r(i)} values,
# with i fixed.  Note \qtilde_{r} =  1/t
for(r in 1:iterates)
{
lc[r]=LC(beta[,r],X[i,],N[i],y[i])
qti[r]=1/lc[r]
}
# standardize \qtilde_{r(i)}
qti=qti/sum(qti)
D1[i]=log(sum(lc*qti))-sum(log(lc)*qti)
D2[i]=mean(log(lc))-log(sum(lc*qti))
D[i]=D1[i]+D2[i]
PPi=JJ%*% qti
PPfi=JJf%*% qti
Dp[i]=sum( (PP-PPi)*(log(PP/(1-PP)) - log(PPi/(1-PPi))) )
Df[i]=sum( (PPf-PPfi)*(log(PPf/(1-PPf)) - log(PPfi/(1-PPfi))) )
}

# Figure 13.11 in text
par(mfrow=c(2,1),mar=c(5, 5, 4, 2) + 0.1)
plot(seq(1:23),Dp,xlab="index",
      ylab=expression(italic(D[i]^p)),cex.lab=1.5)
plot(seq(1:23),Df,xlab="index",
      ylab=expression(italic(D[i]^f)),cex.lab=1.5)
#mtext(expression(D[1]^beta),side=2,cex=1.5,at=1.1)
par(mfrow=c(1,1),mar=c(5, 4, 4, 2) + 0.1)

# Figure not included in text
par(mfrow=c(2,2),mar=c(5, 5, 4, 2) + 0.1)
plot(seq(1:23),D1,ylab=expression(D["1i"]^beta),
                                   xlab="index")
plot(seq(1:23),D2,ylab=expression(D["2i"]^beta),
                                   xlab="index")
plot(seq(1:23),D,ylab=expression(D[i]^beta),
                                   xlab="index")
par(mfrow=c(1,1))