#
# Regresión lineal simple
#
#
# Datos sobre el tiempo necesario para parar coches moviendo a distintas velocidades.
#
rm(list=ls())
help(cars)
data(cars)
plot(cars)
#
# Mirar la correlación
#
cor(cars$speed, cars$dist)
#
# REGRESIÓN CLÁSICA
#
modelolineal <- lm(dist ~ speed, data=cars)  
print(modelolineal)
abline(modelolineal,col='red',lwd=2)
#
# Tabla de ANOVA
#
anova(modelolineal)
#
# Predicción
#
xnuevo <- data.frame(20)
names(xnuevo) <- "speed"
predict(modelolineal,xnuevo,interval="confidence")
predict(modelolineal,xnuevo,interval="prediction")
#
# BAYESIANA CON A PRIOR CONJUGADA
#
m <- c(0,1)
V <- 1000*diag(2)
Vinv <- solve(V)
a <- 1
b <- 1
#
# Valores a posteriori
#
X <- cars[,1]
n <- length(X)
X <- matrix(c(rep(1,n),X),nrow=n)
XTX <- t(X)%*%X
Vpost <- solve(XTX+Vinv)
y <- cars$dist
mpost <- Vpost%*%(Vinv%*%m+t(X)%*%y)
apost <- a+n
bpost <- b+t(m)%*%Vinv%*%m+t(y)%*%y-t(mpost)%*%solve(Vpost)%*%mpost
#
# Simulación Monte Carlo.
#
simuls <- 10000
library(MASS)
phi <- rgamma(simuls,apost/2,bpost/2)
beta <- matrix(rep(NA,simuls*2),nrow=simuls)
for (i in 1:simuls){
  beta[i,] <- mvrnorm(n = 1, mpost, Vpost/phi[i])
}
#
# Media y intervalos de credibilidad para los coeficientes de la regresión
#
apply(beta,2,mean)
apply(beta,2,quantile,probs=c(0.025,0.975))
#
# Predicción de la media
#
xnuevo <- as.numeric(xnuevo)
mediap <- beta[,1]+beta[,2]*xnuevo
mean(mediap)
quantile(mediap,p=c(0.025,0.975))
#
# Predicción de una nueva observación
#
yp <- rnorm(simuls,mediap,1/sqrt(phi))
hist(yp,freq=FALSE,xlab='y',ylab='f',main="")
mean(yp)
quantile(yp,c(0.025,0.975))
#
# BAYESIANA CON A PRIORI NO INFORMATIVA
#
mpost <- modelolineal$coefficients
Vpost <- solve(XTX)
k <- 2
apost <- n-k
bpost <- t(y-X%*%mpost)%*%(y-X%*%mpost)
phi <- rgamma(simuls,apost/2,bpost/2)
beta <- matrix(rep(NA,simuls*2),nrow=simuls)
for (i in 1:simuls){
  beta[i,] <- mvrnorm(n = 1, mpost, Vpost/phi[i])
}
#
# Media y intervalos de credibilidad para los coeficientes de la regresión
#
apply(beta,2,mean)
apply(beta,2,quantile,probs=c(0.025,0.975))
mediap <- beta[,1]+beta[,2]*xnuevo
mean(mediap)
quantile(mediap,p=c(0.025,0.975))
#
# Predicción de una nueva observación
#
yp <- rnorm(simuls,mediap,1/sqrt(phi))
hist(yp,freq=FALSE,xlab='y',ylab='f',main="")
mean(yp)
quantile(yp,c(0.025,0.975))
#
# BAYESIANA CON A PRIORI SEMI CONJUGADA
#
m <- c(0,1)
V <- 1000*diag(2)
a <- 1
b <- 1
#
# SIMULACIONES PARA GIBBS
#
if (!require("MCMCpack")) install.packages("MCMCpack")
library(MCMCpack)
cc <- MCMCregress(dist ~ speed, data=cars, burnin = 1000, 
            mcmc = 10000, thin = 1,
            beta.start = modelolineal$coefficients, 
            marginal.likelihood = c("none"))
summary(cc)
beta <- cc[,1:2]
apply(beta,2,mean)
apply(beta,2,quantile,probs=c(0.025,0.975))
mediap <- beta[,1]+beta[,2]*xnuevo
mean(mediap)
quantile(mediap,p=c(0.025,0.975))
phi <- cc[,3]
yp <- rnorm(simuls,mediap,1/sqrt(phi))
hist(yp,freq=FALSE,xlab='y',ylab='f',main="")
mean(yp)
quantile(yp,c(0.025,0.975))