# =====================================
# Regression model:
# =====================================

hyades = read.table("hyades.dat", header=T) 
x  = hyades$Vmag
y  = hyades$B.V

plot(x,y,xlab="Vmag",ylab="B.V",lwd=2)

# linear model: y = b0 + b1*x
reg1 <- lm(y~x)
summary(reg1)
AIC(reg1)
logLik(reg1)
plot(reg1)

# linear model: y = b0 + b1*x + b2*x^2
x2   = x^2
reg2 = lm(y~x+x2)
summary(reg2)
AIC(reg2)
logLik(reg2)
plot(reg2)


# =====================================
# Normal model:
# =====================================
n   = 30
mu  = 3.5
sig = 1.5
y   = rnorm(n,mu,sig)
plot(y,type="l")

L      = 50
seqmu  = seq(2,5,length=L)
seqsig = seq(1,2,length=L)
conj   = matrix(0,L,L)
for(i in 1:L)
  for(j in 1:L)
    conj[i,j] = prod(dnorm(y,seqmu[i],seqsig[j]))
    
image(seqmu,seqsig,conj,xlab="mu",ylab="sigma")
contour(seqmu,seqsig,conj,add=T)
points(mu,sig,pch=19,cex=2,lwd=2,col=4)

# Regression model:
reg1 = lm(y~1)
summary(reg1)

# ========================================
# Plots: Likelihood, prior and posterior distribution
# p(mu|y) = p(y|mu)p(mu)/p(y)
# ========================================

# Likelihood (fixed sig): p(y|mu)
likemu = rep(0,L)
for(i in 1:L)
  likemu[i] = prod(dnorm(y,seqmu[i],sig))

# ========================================                                       

# prior: p(mu)
priormu = dnorm(seqmu,4,0.5)

# ========================================
# likelihood x prior
likeprior = likemu*priormu

# ~p(y)
delta = seqmu[2]-seqmu[1]
py    = delta*(sum(likeprior) -0.5*(likeprior[1]+likeprior[L]))

# ========================================
# Posterior distribution: 
posterior = likemu*priormu/py

# plot
aux = max(posterior,priormu,likemu/py)
plot(seqmu,priormu,col=2,lwd=3,type="l",ylim=range(0,aux))
lines(seqmu,likemu/py,col=1,lwd=3)
lines(seqmu,posterior,col=4,lwd=3)
legend(2.5,max(likemu/py),c("likelihood","prior","posterior"),col=c(1,2,4),lty=1,lwd=3,bty="n")

# =========================================
# Bayesian inference using MCMC
# Aim: estimate parameters mu and sig
# =========================================


L  = 500                 # number of iterations
y1 = matrix(y,n,1)       # dependent variable
x1 = matrix(1,n,1)       # independent variable

# prior parameters
mp = 5
sp = 0.5^2
a  = b = 0.02
post.mu  = rep(0,L)
post.sig = rep(0,L)

mu.i  = 5
sig.i = 5
post.mu[1]  = mu.i
post.sig[1] = sig.i
for(i in 2:L){
    aux1 = 1/(sig.i^(-2)*t(x1)%*%x1 + sp^(-1)) 
    aux2 = aux1*(sig.i^(-2)*t(x1)%*%y1 + sp^(-1)*mp)
    mu.i = rnorm(1,aux2,sqrt(aux1))
    post.mu[i] = mu.i
    
    aux1  = 0.5*(n+a)
    aux2  = 0.5*(b+t(y1-x1%*%mu.i)%*%(y1-x1%*%mu.i))
    sig.i = sqrt(1/rgamma(1,aux1,rate=aux2))
    post.sig[i] = sig.i
    print(i)
}

par(mfrow=c(2,2))

plot(post.mu,type="l")
abline(h=mu,col=2)
plot(density(post.mu))
abline(v=mu,col=2)

plot(post.sig,type="l")
abline(h=sig,col=2)
plot(density(post.sig))
abline(v=sig,col=2)


summary(post.mu)
summary(post.sig)

