In [1]:
#==============================================================================
# Distribuicao amostral assintotica do estimador de maxima 
# verossimilhanca do probabilidade (parametro) de uma Geometrica 
# (um pequeno estudo de Monte Carlo).
#
#==============================================================================
# Declarando valores para o estudo
#==============================================================================
set.seed(5438)
M      <- 100000         # numero de replicas
In [2]:
#==============================================================================
#==============================================================================
# Supondo probabilidade igual a 0,1 e tamanho de amostra igual a 50
#==============================================================================
#==============================================================================
prob   <- 0.1            # valor verdadeiro (neste estudo) 0.1
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
    #==========================================================================
    # Esta versao da geometrica conta o numero de ensaios  
    # independentes de Bernoulli ate a ocorrencia do primeiro
    # sucesso ( por isto o fator +1 abaixo)
    #==========================================================================
    x[j,] <- rgeom(n,prob) + 1 
}
In [4]:
#==============================================================================
# Calculando a media amostral das M amostras
#==============================================================================
pihat <- 1/apply(x,1,"mean") # vetor de tamanho M
print(length(pihat))

[1] 100000
In [5]:
#==============================================================================
# Analise da distribuicao de pi chapeu
#==============================================================================
desvpad <- sqrt( (prob^2)*(1-prob)/n )
xseq    <- seq(max(0,prob-4*desvpad),min(prob+4*desvpad,1),length=1001)
yseq    <- dnorm(xseq,prob,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,2,2.5),xpd=T,cex.main=2.0)
hist(pihat,prob=T,lwd=2,nclass=30,main="",
  xlim=c(max(0,prob-4*desvpad),min(prob+4*desvpad,1)),
  xlab=expression(hat(pi)),ylab=expression(f(hat(pi))),
  ylim=c(0,1.1*max(yseq)))
lines(xseq,yseq,lwd=4,col="blue")
lines(rep(prob,2),c(0,max(yseq)),lwd=4,col="red")
In [6]:
#==============================================================================
#==============================================================================
# Supondo probabilidade igual a 0,5 e tamanho de amostra igual a 50
#==============================================================================
#==============================================================================
prob   <- 0.5            # valor verdadeiro (neste estudo) 0,5
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
    #==========================================================================
    # Esta versao da geometrica conta o numero de ensaios  
    # independentes de Bernoulli ate a ocorrencia do primeiro
    # sucesso ( por isto o fator +1 abaixo)
    #==========================================================================
    x[j,] <- rgeom(n,prob) + 1 
}
In [8]:
#==============================================================================
# Calculando a media amostral das M amostras
#==============================================================================
pihat <- 1/apply(x,1,"mean") # vetor de tamanho M
print(length(pihat))

[1] 100000
In [9]:
#==============================================================================
# Analise da distribuicao de pi chapeu
#==============================================================================
desvpad <- sqrt( (prob^2)*(1-prob)/n )
xseq    <- seq(max(0,prob-4*desvpad),min(prob+4*desvpad,1),length=1001)
yseq    <- dnorm(xseq,prob,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,2,2.5),xpd=T,cex.main=2.0)
hist(pihat,prob=T,lwd=2,nclass=30,main="",
  xlim=c(max(0,prob-4*desvpad),min(prob+4*desvpad,1)),
  xlab=expression(hat(pi)),ylab=expression(f(hat(pi))),
  ylim=c(0,1.1*max(yseq)))
lines(xseq,yseq,lwd=4,col="blue")
lines(rep(prob,2),c(0,max(yseq)),lwd=4,col="red")
In [10]:
#==============================================================================
#==============================================================================
# Supondo probabilidade igual a 0,1 e tamanho de amostra igual a 1000
#==============================================================================
#==============================================================================
prob   <- 0.1            # valor verdadeiro (neste estudo) 0,1
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
    #==========================================================================
    # Esta versao da geometrica conta o numero de ensaios  
    # independentes de Bernoulli ate a ocorrencia do primeiro
    # sucesso ( por isto o fator +1 abaixo)
    #==========================================================================
    x[j,] <- rgeom(n,prob) + 1 
}
In [12]:
#==============================================================================
# Calculando a media amostral das M amostras
#==============================================================================
pihat <- 1/apply(x,1,"mean") # vetor de tamanho M
print(length(pihat))

[1] 100000
In [13]:
#==============================================================================
# Analise da distribuicao de pi chapeu
#==============================================================================
desvpad <- sqrt( (prob^2)*(1-prob)/n )
xseq    <- seq(max(0,prob-4*desvpad),min(prob+4*desvpad,1),length=1001)
yseq    <- dnorm(xseq,prob,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,2,2.5),xpd=T,cex.main=2.0)
hist(pihat,prob=T,lwd=2,nclass=30,main="",
  xlim=c(max(0,prob-4*desvpad),min(prob+4*desvpad,1)),
  xlab=expression(hat(pi)),ylab=expression(f(hat(pi))),
  ylim=c(0,1.1*max(yseq)))
lines(xseq,yseq,lwd=4,col="blue")
lines(rep(prob,2),c(0,max(yseq)),lwd=4,col="red")
In [14]:
#==============================================================================
#==============================================================================
# Supondo probabilidade igual a 0,5 e tamanho de amostra igual a 1000
#==============================================================================
#==============================================================================
prob   <- 0.5            # valor verdadeiro (neste estudo) 0,5
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
    #==========================================================================
    # Esta versao da geometrica conta o numero de ensaios  
    # independentes de Bernoulli ate a ocorrencia do primeiro
    # sucesso ( por isto o fator +1 abaixo)
    #==========================================================================
    x[j,] <- rgeom(n,prob) + 1 
}
In [16]:
#==============================================================================
# Calculando a media amostral das M amostras
#==============================================================================
pihat <- 1/apply(x,1,"mean") # vetor de tamanho M
print(length(pihat))

[1] 100000
In [17]:
#==============================================================================
# Analise da distribuicao de pi chapeu
#==============================================================================
desvpad <- sqrt( (prob^2)*(1-prob)/n )
xseq    <- seq(max(0,prob-4*desvpad),min(prob+4*desvpad,1),length=1001)
yseq    <- dnorm(xseq,prob,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,2,2.5),xpd=T,cex.main=2.0)
hist(pihat,prob=T,lwd=2,nclass=30,main="",
  xlim=c(max(0,prob-4*desvpad),min(prob+4*desvpad,1)),
  xlab=expression(hat(pi)),ylab=expression(f(hat(pi))),
  ylim=c(0,1.1*max(yseq)))
lines(xseq,yseq,lwd=4,col="blue")
lines(rep(prob,2),c(0,max(yseq)),lwd=4,col="red")