In [1]:
#==============================================================================
# Distribuicao amostral assintotica do estimador de maxima 
# verossimilhanca do parametro (lambda) de uma Poisson 
# (um pequeno estudo de Monte Carlo).
#
#==============================================================================
# Declarando valores para o estudo
#==============================================================================
set.seed(5438)
M      <- 100000         # numero de replicas
In [2]:
#==============================================================================
#==============================================================================
# Supondo lambda igual a 2 e tamanho de amostra igual a 50
#==============================================================================
#==============================================================================
lambda <- 2              # valor verdadeiro (neste estudo) 2
n      <- 50             # tamanho da amostra por replica 50
x      <- matrix(NA,M,n) # para salvar as amostras
In [3]:
#==============================================================================
# Gerando M amostras de tamanho n
#==============================================================================
for(j in 1:M){
    #==========================================================================
    # salvando as amostras de tamanho n por linha
    #==========================================================================
    x[j,] <- rpois(n,lambda)
}
In [4]:
#==============================================================================
# Calculando a media amostral das M amostras
#==============================================================================
lambdahat <- apply(x,1,"mean") # vetor de tamanho M
print(length(lambdahat))

[1] 100000
In [5]:
#==============================================================================
# Analise da distribuicao de lambda chapeu
#==============================================================================
desvpad <- sqrt( lambda/n )
xseq    <- seq(max(0,lambda-4*desvpad),lambda+4*desvpad,length=1001)
yseq    <- dnorm(xseq,lambda,desvpad)

par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,5,5),
    mar=c(5,5.5,2,2.5),xpd=T,cex.main=2.0)
hist(lambdahat,prob=T,lwd=2,nclass=30,main="",
  xlim=c(max(0,lambda-4*desvpad),lambda+4*desvpad),
  xlab=expression(hat(lambda)),ylab=expression(f(hat(lambda))),
  ylim=c(0,1.1*max(yseq)))
lines(xseq,yseq,lwd=4,col="blue")
lines(rep(lambda,2),c(0,max(yseq)),lwd=4,col="red")
In [6]:
#==============================================================================
#==============================================================================
# Supondo lambda igual a 20 e tamanho de amostra igual a 50
#==============================================================================
#==============================================================================
lambda <- 20             # valor verdadeiro (neste estudo) 20
n      <- 50             # tamanho da amostra por replica 50
x      <- matrix(NA,M,n) # para salvar as amostras
In [7]:
#==============================================================================
# Gerando M amostras de tamanho n
#==============================================================================
for(j in 1:M){
    #==========================================================================
    # salvando as amostras de tamanho n por linha
    #==========================================================================
    x[j,] <- rpois(n,lambda)
}
In [8]:
#==============================================================================
# Calculando a media amostral das M amostras
#==============================================================================
lambdahat <- apply(x,1,"mean") # vetor de tamanho M
print(length(lambdahat))

[1] 100000
In [9]:
#==============================================================================
# Analise da distribuicao de lambda chapeu
#==============================================================================
desvpad <- sqrt( lambda/n )
xseq    <- seq(max(0,lambda-4*desvpad),lambda+4*desvpad,length=1001)
yseq    <- dnorm(xseq,lambda,desvpad)

par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,5,5),
    mar=c(5,5.5,2,2.5),xpd=T,cex.main=2.0)
hist(lambdahat,prob=T,lwd=2,nclass=30,main="",
  xlim=c(max(0,lambda-4*desvpad),lambda+4*desvpad),
  xlab=expression(hat(lambda)),ylab=expression(f(hat(lambda))),
  ylim=c(0,1.1*max(yseq)))
lines(xseq,yseq,lwd=4,col="blue")
lines(rep(lambda,2),c(0,max(yseq)),lwd=4,col="red")
In [10]:
#==============================================================================
#==============================================================================
# Supondo lambda igual a 2 e tamanho de amostra igual a 1000
#==============================================================================
#==============================================================================
lambda <- 2              # valor verdadeiro (neste estudo) 2
n      <- 1000           # tamanho da amostra por replica 1000
x      <- matrix(NA,M,n) # para salvar as amostras
In [11]:
#==============================================================================
# Gerando M amostras de tamanho n
#==============================================================================
for(j in 1:M){
    #==========================================================================
    # salvando as amostras de tamanho n por linha
    #==========================================================================
    x[j,] <- rpois(n,lambda)
}
In [12]:
#==============================================================================
# Calculando a media amostral das M amostras
#==============================================================================
lambdahat <- apply(x,1,"mean") # vetor de tamanho M
print(length(lambdahat))

[1] 100000
In [13]:
#==============================================================================
# Analise da distribuicao de lambda chapeu
#==============================================================================
desvpad <- sqrt( lambda/n )
xseq    <- seq(max(0,lambda-4*desvpad),lambda+4*desvpad,length=1001)
yseq    <- dnorm(xseq,lambda,desvpad)

par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,5,5),
    mar=c(5,5.5,2,2.5),xpd=T,cex.main=2.0)
hist(lambdahat,prob=T,lwd=2,nclass=30,main="",
  xlim=c(max(0,lambda-4*desvpad),lambda+4*desvpad),
  xlab=expression(hat(lambda)),ylab=expression(f(hat(lambda))),
  ylim=c(0,1.1*max(yseq)))
lines(xseq,yseq,lwd=4,col="blue")
lines(rep(lambda,2),c(0,max(yseq)),lwd=4,col="red")
In [14]:
#==============================================================================
#==============================================================================
# Supondo lambda igual a 20 e tamanho de amostra igual a 1000
#==============================================================================
#==============================================================================
lambda <- 20             # valor verdadeiro (neste estudo) 20
n      <- 1000           # tamanho da amostra por replica 1000
x      <- matrix(NA,M,n) # para salvar as amostras
In [15]:
#==============================================================================
# Gerando M amostras de tamanho n
#==============================================================================
for(j in 1:M){
    #==========================================================================
    # salvando as amostras de tamanho n por linha
    #==========================================================================
    x[j,] <- rpois(n,lambda)
}
In [16]:
#==============================================================================
# Calculando a media amostral das M amostras
#==============================================================================
lambdahat <- apply(x,1,"mean") # vetor de tamanho M
print(length(lambdahat))

[1] 100000
In [17]:
#==============================================================================
# Analise da distribuicao de lambda chapeu
#==============================================================================
desvpad <- sqrt( lambda/n )
xseq    <- seq(max(0,lambda-4*desvpad),lambda+4*desvpad,length=1001)
yseq    <- dnorm(xseq,lambda,desvpad)

par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(10,5,5),
    mar=c(5,5.5,2,2.5),xpd=T,cex.main=2.0)
hist(lambdahat,prob=T,lwd=2,nclass=30,main="",
  xlim=c(max(0,lambda-4*desvpad),lambda+4*desvpad),
  xlab=expression(hat(lambda)),ylab=expression(f(hat(lambda))),
  ylim=c(0,1.1*max(yseq)))
lines(xseq,yseq,lwd=4,col="blue")
lines(rep(lambda,2),c(0,max(yseq)),lwd=4,col="red")