
#==============================================================================
# Gerando uma amostra aleatoria de uma populacao normal com media MU e 
# variancia SIGMA2 (desvio padrao SIGMA)
#==============================================================================
set.seed(12345)   # fixa a sequencia de numeros (pseudo-aleatorios)
n      <- 10      # tamanho da amostra aleatoria
MU     <- 5       # valor verdadeiro da media (desconhecido na pratica)
SIGMA  <- 3       # valor verdadeiro do desvio padrao (supor conhecido)
SIGMA2 <- SIGMA^2 # valor verdadeiro da variancai (supor conhecido)
x      <- rnorm(n,MU,SIGMA) # gerando a amostra
print(x,digits=3)

 [1]  6.757  7.128  4.672  3.640  6.818 -0.454  6.890  4.171  4.148  2.242


#==============================================================================
# Algumas estatisticas descritivas
#==============================================================================
library(psych)
options(digits=4)
x.describe      <- describe(x)
x.media         <- x.describe[,3]
x.desvpad       <- x.describe[,4]
x.mediana       <- x.describe[,5]
x.assimetria    <- x.describe[,11]
x.exc.curtose   <- x.describe[,12]
STATS           <- cbind(x.media,x.desvpad,x.mediana,x.assimetria,
                         x.exc.curtose)
colnames(STATS) <- c("Média","Dev.Pad.","Mediana","Assimetria","Exec. curtose")
STATS


#==============================================================================
# Funcao de verossimilhanca do modelo exponencial
#==============================================================================
fvero.norm <- function(mu,x,sigma2,logT=T){
    if(!is.vector(x) || !is.vector(mu)) return(NA)
    m    <- length(mu)
    out  <- rep(NA,m)
    for(i in 1:m)
        out[i] <- sum(dnorm(x,mu[i],SIGMA,log=T))        
    if(!logT) out <- exp(out)
    return(out)
}


#==============================================================================
# Grafico da funcao de verossimilhanca do modelo normal com media 
# desconhecida e variancia conhecida
#==============================================================================
xvero <- seq(-12,12,length=5000) # sequencia de valores de lambda
yvero   <-  fvero.norm(xvero,x,SIGMA2,logT=F)

par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(6,4,5),
    mar=c(5,5,2,2.5),xpd=T,cex.main=2.0)
plot(xvero,yvero,type="l",lwd=2,col="blue",xlab=expression(mu),
    ylab=expression(L(mu)))


#==============================================================================
# Distribuicao a priori, funcao de verossimilhanca e distribuicao a 
# posteriori do modelo normal
# Consideramos uma priori normal com media 0 (zero) e variancia 100.
#==============================================================================
a      <- 0   # media a priori
b      <- 3  # desvio padrao a priori (variancia a priori igual a 100)
b2     <- b^2
tau2   <- 1/( n/SIGMA2 + 1/b2 )        # variancia a posteriori
tau    <- sqrt(tau2)                   # desvio padrao a posteriori
xi     <- tau2*(sum(x)/SIGMA2 + a/b2 ) # media a posteriori
yprior <- dnorm(xvero,a,b)
ypost  <- dnorm(xvero,xi,tau)

# Constante auxiliar para o grafico
cst    <- sum(yvero)*(xvero[2]-xvero[1])  # integracao numerica de Riemann

par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(6,6,5),
    mar=c(5,5,2,2.5),xpd=T,cex.main=2.0)
plot(xvero,yvero/cst,type="l",lwd=2,col="blue",xlab=expression(mu),
    ylab="",ylim=c(0,0.5))
lines(xvero,yprior,lwd=2,col="red")
lines(xvero,ypost,lwd=3,lty=2,col="black")
legend(-12,0.4,legend=c("Priori","Verossimilhança","Posteriori"),
    lwd=c(2,2,3),lty=c(1,1,2),col=c("red","blue","black"),bty="n")