#==============================================================================
# Distribuicao preditiva (x_n | x_1,x_2,...,x_{n-1})
#==============================================================================
x      <- c(2911,3403,3237,3509,3118) 
xbarra <- mean(x)

#==============================================================================
# Funcao de densidade de probabilidade preditiva (x_n | x_1,x_2,...,x_{n-1})
#==============================================================================
fpred <- function(xn,xbarra,a,b,logT=T){
    if(!is.vector(xbarra) || !is.vector(a) || !is.vector(b)) return(NA)
    if(length(xbarra)!=1 || length(a)!=1 || length(b)!=1) return(NA)
    n    <- length(xn)
    out  <- ( (a+n-1)*log(b+(n-1)*xbarra)-(a+n)*log(b+(n-1)*xbarra+xn)
            + lgamma(a+n)-lgamma(a+n-1) )
    if(!logT) out <- exp(out)
    return(out)
}

#==============================================================================
# Grafico da f.d.p. preditiva (x_n | x_1,x_2,...,x_{n-1})
#==============================================================================
xpred <- seq(0,30000,length=10000) # sequencia de valores de lambda
ypred <-  fpred(xpred,xbarra=xbarra,a=4,b=20000,logT=F)

par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(4,4,5),
    mar=c(5,5,2,2.5),xpd=T,cex.main=2.0)
plot(xpred,ypred,type="l",lwd=2,col="blue",xlab=expression(x[n]),
    ylab=expression(paste("f(",x[n],"|",x[1],",...,",x[n-1],")",sep="")))