In [1]:
#==============================================================================
# Distribuicao amostral assintotica do estimador de maxima verossimilhanca 
# da probabilidade (parametro) de uma Bernoulli 
# (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
    #==========================================================================
    x[j,] <- rbinom(n,1,prob)
}
In [4]:
#==============================================================================
# Calculando a media amostral das M amostras
#==============================================================================
pihat <- apply(x,1,"mean") # vetor de tamanho M
print(length(pihat))

[1] 100000
In [5]:
#==============================================================================
# Analise da distribuicao de pi chapeu
#==============================================================================
desvpad <- sqrt( prob*(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))))
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
    #==========================================================================
    x[j,] <- rbinom(n,1,prob)
}
In [8]:
#==============================================================================
# Calculando a media amostral das M amostras
#==============================================================================
pihat <- apply(x,1,"mean") # vetor de tamanho M
print(length(pihat))

[1] 100000
In [9]:
#==============================================================================
# Analise da distribuicao de pi chapeu
#==============================================================================
desvpad <- sqrt( prob*(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))))
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
    #==========================================================================
    x[j,] <- rbinom(n,1,prob)
}
In [12]:
#==============================================================================
# Calculando a media amostral das M amostras
#==============================================================================
pihat <- apply(x,1,"mean") # vetor de tamanho M
print(length(pihat))

[1] 100000
In [13]:
#==============================================================================
# Analise da distribuicao de pi chapeu
#==============================================================================
desvpad <- sqrt( prob*(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))))
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
    #==========================================================================
    x[j,] <- rbinom(n,1,prob)
}
In [16]:
#==============================================================================
# Calculando a media amostral das M amostras
#==============================================================================
pihat <- apply(x,1,"mean") # vetor de tamanho M
print(length(pihat))

[1] 100000
In [17]:
#==============================================================================
# Analise da distribuicao de pi chapeu
#==============================================================================
desvpad <- sqrt( prob*(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))))
lines(xseq,yseq,lwd=4,col="blue")
lines(rep(prob,2),c(0,max(yseq)),lwd=4,col="red")