rm(list=ls())
set.seed(12345)

# simulando uma regressao linear simples
n     = 200
sig2  = 1.0
alpha = 0.0
beta  = 2.0
sig   = sqrt(sig2)
x     = rnorm(n,0,1)
erro  = rnorm(n,0,sig)
y     = alpha + beta*x + erro

plot(x,y)

# Combinando as colunas x e y
matrix.xy = cbind(x,y)

# salvando a matriz de dados matri.xy
write(t(matrix.xy),"myfirstdataset.txt",ncolumns=2)


# vendo os dados

x[1:10]

y[190:n]

matrix.xy[1:3,]


# Regressao OLS
X = cbind(1,x)

b.hat = solve(t(X)%*%X)%*%t(X)%*%y

yfit = X%*%b.hat

ehat = (y-yfit)

S2e = sum((y-yfit)^2)/n


# Grafico dos observados e ajustados
plot(x,y)
points(x,yfit,col=2,pch=16)


# Nossa primeira inferencia Bayesiana
b0 = matrix(0,2,1)
B0 = diag(100,2)

B1 = solve(solve(B0) + t(X)%*%X)

b1 = B1%*%(solve(B0)%*%b0 + t(X)%*%y)

cbind(b.hat,b1)

# Nossa segunda inferencia Bayesiana (tight priors)
B00 = diag(0.0001,2)
B2  = solve(solve(B00) + t(X)%*%X)
b2  = B2%*%(solve(B00)%*%b0 + t(X)%*%y)

# Comparando as estimativas
cbind(b.hat,b1,b2)



# Grafico da posterior para a inclinacao
L = b1[2] - 4*sqrt(sig2*B1[2,2])
U = b1[2] + 4*sqrt(sig2*B1[2,2])
b1s = seq(L,U,length=500)

plot(b1s,dnorm(b1s,b1[2],sqrt(sig2*B1[2,2])),type="l",
     xlab=expression(beta[2]),ylab="Densidade",main="")
title("Posteriori para a inclinacao")



# Minha primeira funcao
# Avaliar a normal bivariada no ponto x com media b e 
# matrix de precisao P
dnorm2 = function(x,b,P){
  1/(2*pi)*sqrt(det(P))*exp(-0.5*t(x-b)%*%P%*%(x-b))
}


x1 = c(0,0)
b = c(0,0)
P = diag(1,2)
dnorm2(x1,b,P)


# Grafico da superficie p(alpha,beta|X,y)
L1  = b1[1] - 4*sqrt(sig2*B1[1,1])
U1  = b1[1] + 4*sqrt(sig2*B1[1,1])
L2  = b1[2] - 4*sqrt(sig2*B1[2,2])
U2  = b1[2] + 4*sqrt(sig2*B1[2,2])
b1s = seq(L1,U1,length=100)
b2s = seq(L2,U2,length=100)

P = solve(sig2*B1)
post = matrix(0,100,100)
for (i in 1:100)
  for (j in 1:100)
    post[i,j] = dnorm2(c(b1s[i],b2s[j]),b1,P)


contour(b1s,b2s,post,xlab=expression(alpha),ylab=expression(beta))
points(b.hat[1],b.hat[2],col=2,pch=16,cex=3)


image(b1s,b2s,post,xlab=expression(alpha),ylab=expression(beta))
contour(b1s,b2s,post,add=TRUE,nlevels=20,drawlabels=FALSE)


persp(b1s,b2s,post,xlab=expression(alpha),ylab=expression(beta))



pdf(file="myfirstpdfgraph.pdf",height=10,width=10)
persp(b1s,b2s,post, theta = 30, phi = 30, expand = 0.5, 
      col = "lightblue",
      ltheta = 120, shade = 0.75, ticktype = "detailed",
      xlab = "alpha", ylab = "beta", zlab = "posteriori")
dev.off()




lm.bayes = function(y,x,sig2,b0,B0){
  X  = cbind(1,x)
  B1 = solve(solve(B0) + t(X)%*%X)
  b1 = B1%*%(solve(B0)%*%b0 + t(X)%*%y)
  B11 = sig2*B1
  return(list(mean=b1,var=B11))
}


regbayes = lm.bayes(y,x,sig2,b0,B0)



names(regbayes)

regbayes$mean

regbayes$var