# schiz.f.r
#
# Will try marginal posterior ....
# log marginal posterior distribution for p(alpha,phi|y) for the
# schizophrenic example, p 431
#
# THETA is a 23xK matrix with, each column representing a single
# occurence of:
#     (alpha,sigmaa,beta,lambda,tao,mu,sigmay), with
#     length(alpha)=17
# In all, there are K occurences of these tuples.
#f<-function(theta)
f<-function(theta,dolog=FALSE)

{

#for testing
#theta<-seq(0,0,length=23)
#theta[1:17]<- c(310.6914,359.7207,306.3999,304.1072,270.2453,313.9323,352.1958,273.1259,263.6610,323.5660,309.1630,321.1054,422.9851,480.1710,538.6872,488.4738,461.9173)
#theta[18:23]<-c(46.81092,144.3316,0.07207053,758.19298,307.8917,93.86176)

# Data
m<-30
n<-17
nnon<-11

if (is.null(dim(theta))) {
  K<-1
  alpha<-theta[1:17]
  sigmaa<-theta[18]
  beta<-theta[19]
  lambda<-theta[20]
  tao<-theta[21]
  mu<-theta[22]
  sigmay<-theta[23]
}
else {
  K<-dim(theta)[2]
  alpha<-theta[1:17,1:K]
  sigmaa<-theta[18,1:K]
  beta<-theta[19,1:K]
  lambda<-theta[20,1:K]
  tao<-theta[21,1:K]
  mu<-theta[22,1:K]
  sigmay<-theta[23,1:K]
}

#x<-matrix(scan("schiz.dat"),m,n)    # columns are the measurements for
                                    # the individuals
xmat<-matrix(x,m,n*K)     # block matrix with data X: XMAT=[X|X|X|...|X]

alpham<-t(matrix(alpha,n*K,m))
sigmaymat1<-t(matrix(sigmay,K,n))
sigmaymat<-t(matrix(sigmaymat1,K*n,m))
mumat1<-t(matrix(mu,K,n))
mumat<- t(matrix(mumat1,K*n,m))
betamat1<-t(matrix(beta,K,n))
betamat<-t(matrix(betamat1,K*n,m))
sigmaamat1<-t(matrix(sigmaa,K,n))
sigmaamat<-t(matrix(sigmaamat1,K*n,m))
lambdamat1<-t(matrix(lambda,K,n))
lambdamat<-t(matrix(lambdamat1,K*n,m))
taomat1<-t(matrix(tao,K,n))
taomat<-t(matrix(taomat1,K*n,m))
Schiz<-seq(1,1,length=n)
Schiz[1:nnon]<-0
Schizmat1<-matrix(Schiz,n,K)
Schizmat<-t(matrix(Schiz,K*n,m))

term1<-dnorm(alpha,(mumat1+betamat1*Schizmat1),sigmaamat1,log=dolog)
if (dolog) {
   term2<-log(apply(((1-lambdamat*Schizmat)*dnorm(xmat,alpham,sigmaymat)+lambdamat*Schizmat*dnorm(xmat,(alpham+taomat),sigmaymat)),2,prod))
   term2<-matrix(term2,n,K)
   return(apply((term1+term2),2,sum))
 }
else {
  term2<-apply(((1-lambdamat*Schizmat)*dnorm(xmat,alpham,sigmaymat)+lambdamat*Schizmat*dnorm(xmat,(alpham+taomat),sigmaymat)),2,prod)
   term2<-matrix(term2,n,K)
   return(apply((term1*term2),2,prod))
}

}