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)