In [1]:
#==============================================================================
# Distribuicao amostral assintotica do estimador de maxima 
# verossimilhanca do desvio padrao de uma normal com media 0 e 
# variancia sigma^2 (um pequeno estudo de Monte Carlo).
#
#==============================================================================
# Declarando valores para o estudo
#==============================================================================
set.seed(5438)
M      <- 100000         # numero de replicas
mu     <- 0              # valor verdadeiro (neste estudo)
sigma2 <- 9              # valor verdadeiro (neste estudo)
sigma  <- sqrt(sigma2)   # valor verdadeiro (neste estudo)
In [2]:
#==============================================================================
# Raiz quadrada da media dos quadrados ordinarios
#==============================================================================
sqrtmeansq <- function(x){
    sqrt(mean(x^2))
}
In [3]:
#==============================================================================
#==============================================================================
# Supondo tamanho de amostra igual a 5
#==============================================================================
#==============================================================================
n      <- 5              # tamanho da amostra por replica 5
x      <- matrix(NA,M,n) # para salvar as amostras
In [4]:
#==============================================================================
# Gerando M amostras de tamanho n
#==============================================================================
for(j in 1:M){
    #==========================================================================
    # salvando as amostras de tamanho n por linha
    #==========================================================================
    x[j,] <- rnorm(n,mu,sigma)
}
In [5]:
#==============================================================================
# Calculando a raiz quadrada da media dos quadrados ordinarios das
# M amostras
#==============================================================================
sigmahat <- apply(x,1,"sqrtmeansq") # vetor de tamanho M
print(length(sigmahat))

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

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(sigmahat,prob=T,lwd=2,nclass=30,main="",
  xlim=c(max(0,sigma-1.1*aux),sigma+1.1*aux),
  xlab=expression(hat(sigma)),ylab=expression(f(hat(sigma))))
lines(xseq,yseq,lwd=4,col="blue")
lines(rep(sigma,2),c(0,max(yseq)),lwd=4,col="red")
In [7]:
#==============================================================================
#==============================================================================
# Supondo tamanho de amostra igual a 10
#==============================================================================
#==============================================================================
n      <- 10             # tamanho da amostra por replica 10
x      <- matrix(NA,M,n) # para salvar as amostras
In [8]:
#==============================================================================
# Gerando M amostras de tamanho n
#==============================================================================
for(j in 1:M){
    #==========================================================================
    # salvando as amostras de tamanho n por linha
    #==========================================================================
    x[j,] <- rnorm(n,mu,sigma)
}
In [9]:
#==============================================================================
# Calculando a raiz quadrada da media dos quadrados ordinarios das
# M amostras
#==============================================================================
sigmahat <- apply(x,1,"sqrtmeansq") # vetor de tamanho M
print(length(sigmahat))

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

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(sigmahat,prob=T,lwd=2,nclass=30,main="",
  xlim=c(max(0,sigma-1.1*aux),sigma+1.1*aux),
  xlab=expression(hat(sigma)),ylab=expression(f(hat(sigma))))
lines(xseq,yseq,lwd=4,col="blue")
lines(rep(sigma,2),c(0,max(yseq)),lwd=4,col="red")
In [11]:
#==============================================================================
#==============================================================================
# Supondo tamanho de amostra igual a 100
#==============================================================================
#==============================================================================
n      <- 100            # tamanho da amostra por replica 100
x      <- matrix(NA,M,n) # para salvar as amostras
In [12]:
#==============================================================================
# Gerando M amostras de tamanho n
#==============================================================================
for(j in 1:M){
    #==========================================================================
    # salvando as amostras de tamanho n por linha
    #==========================================================================
    x[j,] <- rnorm(n,mu,sigma)
}
In [13]:
#==============================================================================
# Calculando a raiz quadrada da media dos quadrados ordinarios das
# M amostras
#==============================================================================
sigmahat <- apply(x,1,"sqrtmeansq") # vetor de tamanho M
print(length(sigmahat))

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

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(sigmahat,prob=T,lwd=2,nclass=30,main="",
  xlim=c(max(0,sigma-1.1*aux),sigma+1.1*aux),
  xlab=expression(hat(sigma)),ylab=expression(f(hat(sigma))))
lines(xseq,yseq,lwd=4,col="blue")
lines(rep(sigma,2),c(0,max(yseq)),lwd=4,col="red")