#
# Simulación de una distribución normal-gamma
#
# Parámetros
#
m <- 0
c <- 1
a <- 5
b <- 2
#
# Simulaciones
#
simuls <- 10000
nc <- 50
#
phi <- rgamma(simuls,a/2,b/2)
mu <- rnorm(simuls,m,1/sqrt(phi*c))
maxh <- max(hist(mu,nclass=nc,plot=FALSE)$density)
#
# Verdadera densidad de mu
#
rejillamu <- min(mu)+c(0:1000)*(max(mu)-min(mu))/1000
fmu <- dt((rejillamu-m)/sqrt(b/(a*c)),a)*sqrt(a*c/b)
maxf <- max(fmu)
#
# Histograma de los valores generados.
#
hist(mu,freq=FALSE,nclass=nc,xlab=expression(mu),ylab='f',ylim=c(0,max(maxh,maxf)),main='')
#
# Incluir la densidad verdadera de mu.
#
lines(rejillamu,fmu,lwd=2,col='green')
#
# Simular de la densidad predictiva.
#
y <- rnorm(simuls,mu,1/sqrt(phi))
maxhy <- max(hist(y,nclass=nc,plot=FALSE)$density)
rejillay <- min(y)+c(0:1000)*(max(y)-min(y))/1000
cy <- c/(c+1)
fy <- dt((rejillay-m)/sqrt(b/(a*cy)),a)*sqrt(a*cy/b)
maxfy <- max(fy)
hist(y,freq=FALSE,nclass=nc,xlab='y',ylab='f',ylim=c(0,max(maxhy,maxfy)),main='')
lines(rejillay,fy,type='l',lwd=2,col='blue')