In [1]:
#==============================================================================
# Gerando uma amostra aleatoria de uma populacao exponencial
# com parametro lambda (media 1/lambda). 
#==============================================================================
set.seed(12345)  # fixa a sequencia de numeros (pseudo-aleatorios)
n      <- 100    # tamanho da amostra aleatoria
lambda <- 0.0003 # valor verdadeiro do parametro (desconhecido na pratica)
x      <- rexp(n,lambda) # gerando a amostra
print(x,digits=3)

  [1]  1472.7  1758.2  2694.9    61.5  1518.0    79.8 21340.6  4193.2  3207.3
 [10]  6060.0   755.0  1489.2  4174.9  1323.6   293.7  5740.5 17949.5  6293.3
 [19]  1212.2  3911.9  3798.8  1426.6  4861.8  1881.3  4097.1  2439.6  3260.2
 [28]  9342.4   837.0   727.9  5774.5  1531.8  7882.7  2857.1  2426.3   505.9
 [37]   345.3  1168.0  3949.6  3454.6 11055.4   793.0  1179.0  5290.3    20.0
 [46] 12314.0  5141.3  3044.2  2171.0    15.6  2023.8  9424.1   672.9  1431.9
 [55]    92.3  1467.4  1666.3  8699.0  4281.5  2903.4   781.7  1204.6   916.4
 [64]   496.9   729.8  4089.8  1701.8  4040.3  1966.5  1358.4 14567.6  3254.4
 [73]   396.4    99.5  8382.8  1221.7  1399.8  2001.5   431.2  1828.4  9476.4
 [82]  1146.1  2085.1  1087.7  1611.3  6222.1  8800.5  2051.7   145.2  1290.3
 [91]  4960.6  1147.0   859.8   393.6  1971.1  3670.8  6848.4  1898.5  1220.3
[100]  3701.3
In [2]:
#==============================================================================
# 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
A matrix: 1 × 5 of type dbl
MédiaDev.Pad.MedianaAssimetriaExec. curtose
3332375019332.3086.417
In [3]:
#==============================================================================
# Histograma da amostra (observada) dos tempos de vida
#==============================================================================
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)
hist(x,lwd=2,col="darkorange",main="",xlab="Tempos de vida",ylab="Densidade",
     prob=T,xlim=c(0,25000))
In [4]:
#==============================================================================
# Box-plot da amostra (observada) dos tempos de vida
#==============================================================================
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)
boxplot(x,lwd=2,col="darkorange",main="",ylab="Tempos de vida",pch=16,
    ylim=c(0,25000))
In [5]:
#==============================================================================
# Funcao de verossimilhanca do modelo exponencial
#==============================================================================
fvero.exp <- function(lambda,x,logT=T){
    if(!is.vector(x)) return(NA)
    n    <- length(x)
    out  <- n*log(lambda)-lambda*sum(x) 
    if(!logT) out <- exp(out-max(out))
    return(out)
}
In [6]:
#==============================================================================
# Grafico da funcao de verossimilhanca do modelo exponencial
#==============================================================================
xvero <- seq(0.000001,0.001,length=1000) # sequencia de valores de lambda
yvero   <-  fvero.exp(xvero,x,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(lambda),
    ylab=expression(L(lambda)))
In [7]:
#==============================================================================
# Distribuicao a priori, funcao de verossimilhanca e distribuicao a 
# posteriori do modelo exponencial
#==============================================================================
a      <- 4      # priori: parametro de forma
b      <- 20000  # priori: parametro de escala
yprior <- dgamma(xvero,a,b)
ypost  <- dgamma(xvero,a+n,b+sum(x))

# 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,4,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(lambda),
    ylab="",ylim=c(0,15000))
lines(xvero,yprior,lwd=2,col="red")
lines(xvero,ypost,lwd=3,lty=2,col="black")
legend(6e-4,13000,legend=c("Priori","Verossimilhança","Posteriori"),
    lwd=c(2,2,3),lty=c(1,1,2),col=c("red","blue","black"),bty="n")