softdrinks<-read.table('ex1_softdrin.txt',header=TRUE) softdrinks y<-softdrinks$Time x<-matrix(ncol=3,nrow=25) x[,1]<-1 x[,2]<-softdrinks[,2] x[,3]<-softdrinks[,3] x1<-x[,1] x2<-x[,1:2] x3<-x[,c(1,3)] x4<-x[,1:3] n<-length(y) c<-sqrt(10) c<-exp(350) #test constant <- a*log(b)-lgamma(a)-0.5*n*log(2*3.14)+lgamma(0.5*n+a) # ------------------model 1 ----------------------------- x0<-x1 # design matrix Pm<-1 # number of parameters in the linear predictor # ------------------------------------------------------- xtx <- t(x0)%*%x0 # xtx matrix a<-0.001 b<-0.001 mu<-rep(0,Pm) # prior mean PS <- diag(c^2,nrow=Pm) # prior variance IPS <- solve(PS) # inverse prior variance/prior precision matrix IS <- xtx + IPS # inverse S tilde S <- solve(IS) # S tilde # mle estimate beta.mle <- solve(xtx) %*% t(x0) %*% y # posterior mode beta.post <- S %*%( xtx %*% beta.mle + IPS%*%mu) # ss ss <- t(y)%*%y - t(beta.post)%*%IS%*%beta.post + t(mu)%*%IPS%*%mu logmar <- -Pm*log(c) + 0.5*log( det(S) ) + ( -0.5*n-a)*log(0.5*ss+b)+constant logmar1<-logmar # ------------------------------------------------------- # ------------------model 2 ----------------------------- x0<-x2 # design matrix Pm<-2 # number of parameters in the linear predictor # ------------------------------------------------------- xtx <- t(x0)%*%x0 # xtx matrix a<-0.001 b<-0.001 mu<-rep(0,Pm) # prior mean PS <- diag(c^2,nrow=Pm) # prior variance IPS <- solve(PS) # inverse prior variance/prior precision matrix IS <- xtx + IPS # inverse S tilde S <- solve(IS) # S tilde # mle estimate beta.mle <- solve(xtx) %*% t(x0) %*% y # posterior mode beta.post <- S %*%( xtx %*% beta.mle + IPS%*%mu) # ss ss <- t(y)%*%y - t(beta.post)%*%IS%*%beta.post + t(mu)%*%IPS%*%mu logmar <- -Pm*log(c) + 0.5*log( det(S) ) + ( -0.5*n-a)*log(0.5*ss+b)+constant logmar2<-logmar # ------------------------------------------------------- # ------------------model 3 ----------------------------- x0<-x3 # design matrix Pm<-2 # number of parameters in the linear predictor # ------------------------------------------------------- xtx <- t(x0)%*%x0 # xtx matrix a<-0.001 b<-0.001 mu<-rep(0,Pm) # prior mean PS <- diag(c^2,nrow=Pm) # prior variance IPS <- solve(PS) # inverse prior variance/prior precision matrix IS <- xtx + IPS # inverse S tilde S <- solve(IS) # S tilde # mle estimate beta.mle <- solve(xtx) %*% t(x0) %*% y # posterior mode beta.post <- S %*%( xtx %*% beta.mle + IPS%*%mu) # ss ss <- t(y)%*%y - t(beta.post)%*%IS%*%beta.post + t(mu)%*%IPS%*%mu logmar <- -Pm*log(c) + 0.5*log( det(S) ) + ( -0.5*n-a)*log(0.5*ss+b)+constant logmar3<-logmar # ------------------------------------------------------- # ------------------model 4 ----------------------------- x0<-x4 # design matrix Pm<-3 # number of parameters in the linear predictor # ------------------------------------------------------- xtx <- t(x0)%*%x0 # xtx matrix a<-0.001 b<-0.001 mu<-rep(0,Pm) # prior mean PS <- diag(c^2,nrow=Pm) # prior variance IPS <- solve(PS) # inverse prior variance/prior precision matrix IS <- xtx + IPS # inverse S tilde S <- solve(IS) # S tilde # mle estimate beta.mle <- solve(xtx) %*% t(x0) %*% y # posterior mode beta.post <- S %*%( xtx %*% beta.mle + IPS%*%mu) # ss ss <- t(y)%*%y - t(beta.post)%*%IS%*%beta.post + t(mu)%*%IPS%*%mu logmar <- -Pm*log(c) + 0.5*log( det(S) ) + ( -0.5*n-a)*log(0.5*ss+b)+constant logmar4<-logmar # ------------------------------------------------------- lmarginal<-c(logmar1,logmar2,logmar3,logmar4) mar<-exp(lmarginal) lmarginal mar/sum(mar)