Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Comparação do histograma da função densidade após a geração de variáveis de forma randômicas.

0 votos
21 visitas
perguntada Out 5 em Economia por CICERO FILHO (26 pontos)  

Gere as variáveis aleatórias para os casos abaixo, e compare o histograma com as respectivas funções de densidade:

(a) Normal usando amostragem por rejeição com candidatos de uma Cauchy.

(b) Variáveis aleatórias Gama (4,3; 6,2) usando Gama(4,7).

(c) Normal Padrão Truncada para (2,∞).

Referência: Exercício 2.20 - Cap 2 – do livro “Monte Carlo Statistical Methods"- 2ed - pag. 66 - Christian P. Robert e George Casella.

Compartilhe

1 Resposta

0 votos
respondida Out 5 por CICERO FILHO (26 pontos)  
editado Out 6 por CICERO FILHO

a) Para resolver a questão, foram utilizados métodos de amostragem que servem exatamente para gerar amostras de uma determinada variável aleatória. A Amostragem por Rejeição é um dos métodos que fazem isso, pois é bem intuitivo e consiste basicamente em dois passos. O primeiro consiste em gerar uma variável aleatória usando uma distribuição auxiliar. E segundo realizar o teste para verificar se essa variável aleatória poderia ter vindo da distribuição de interesse. Caso positivo, é aceita. Caso contrário, é rejeitada.

Considerando uma normal padrão f e uma Cauchy padrão g.

A densidade da distribuição normal padrão é:

\[f(x)=1\sqrt2\pi\ exp \frac{-x^2}{2}\]

A densidade da distribuição padrão de Cauchy é:

\[f(x) = \frac{1}{\pi\left(1+x^2\right)}\]

O método de aceitar-rejeitar requer:

\[M\geq\frac{fx\left(x\right)}{gx\left(x\right)}=\frac{\pi\left(1+x^2\right)}{\sqrt{2\pi}}exp\left\{\frac{-x^2}{2}\right\}=h(x)\]

O menor limite superior de h (x) é:

\[\frac{d\log{h\left(x\right)}}{dx}=-x+\frac{2x}{1+x^2}=0\]

Assim,

\[x=1\]
\[M\geq h\left(1\right)=\sqrt{\frac{2\pi}{e}}\]

O uso de uma distribuição de proposta normal padrão para simular a distribuição de Cauchy padrão requer:

\[h\left(x\right)=\frac{\sqrt{2\pi}}{x(1+x^2)}exp\left\{\frac{x^2}{2}\right\}\le M\]

No entanto, uma função exponencial cresce mais rápido do que uma função polinomial quando:

\[x \rightarrow \infty \]

de modo que não há ponto de máximo.

Resolução computacional letra (a):

n=10000
k=1.6
m=1
x <- seq(m-5,m+5,.001)
f <- 1.6/(pi*(1+(x-m)^2))
f <- k*dcauchy(x,location = m)
plot(x,dnorm(x,m,1),type="l",ylim=c(0,.5),xlab="x",ylab="Densidade", main =" ")
linhas(x,f,lty=4)
xc <- rep(0,n)
yc <- rep(0,n)
para(i in 1:n){
       rejeitar<- VERDADEIRO
       Enquanto(rejeitar){
           p <- runif(1)
           if (p < 0.5){
           xc[i] <- qcauchy(p/k,location=m)+m
       }
      else{
         xc[i] <- -(qcauchy((1-p)/k,location=m))+m
       }
      p <- runif(1)
      yc[i] <- p*k*dcauchy(xc[i],location = m)
      Se (yc[i] < dnorm(xc[i])){
          rejeitar<- FALSO
      }
    }
  }
  points(xc,yc)
  hist(xc,probability=T,xlim=c(m-4,m+4),ylim=c(0,1),xlab="x", breaks = 32, main =" ")
  lines(x,dnorm(x,m,1))

A imagem será apresentada aqui.

b) Tanto f como g são distribuições Gama, mas com parâmetros distintos.

Assim, sendo f uma Gama 1 (a, b) e g uma Gama 2 \((a^\ast,b^\ast)\) então, ao se considerar \(a>a^\ast\), é possível verificar que não há ponto de máximo caso \(b\geq b^\ast\), tendo em vista que o limite de x, para a função h(x), tende ao infinito. Com isso pode-se considerar \(b^\ast>b\).

Portanto, dado que \(a>a^\ast\) e \(b^\ast>b\), os dados apresentados Gamma (4,3; 6,2) e Gamma (4; 7), são satisfeitos, pois 4,3>4 e 6,2<7.

Resolução computacional letra (b):

n=1000
alpha = 4.3 #shape
beta = 6.2 #rate
k = 2.8 #constant
x <- seq(0,10,length.out=n)
f <- dgamma(x,shape=alpha, rate = beta)
plot(x,k*dgamma(x,shape=4, rate=7),type="l",ylim=c(0,5),xlab="x",ylab="Densidade", main = "") # Plots the gamma(4,7)
lines(x,f,lty=4) # Representa a função alvo
xc <- rep(0,n)
yc <- rep(0,n)
para (i in 1:n){
       rejeitar <- VERDADEIRA
       Enquanto (rejeitar){
             p <- runif(1)
             if (p < 0.5){
                 xc[i]<- qgamma(p,shape = 4, rate=7)
             }
             else{
                xc[i] <- qgamma(p,shape = 4,rate = 7)
             }
            p <- runif(1)
            yc[i] <- p*k*dgamma(xc[i], shape= 4,rate = 7)
            if (yc[i] < dgamma(xc[i],shape=alpha,rate=beta)){ # Então o ponto sob o ponto gamma(alpha,beta)
                rejeitar <- FALSO
           }
      }
 }
pontos (xc,yc,type='p')   # Os pontos não rejeitados
hist(xc,probabilidade=T,xlim=c(0,10),ylim=c(0,2),xlab="x",col=16, breaks= (n/16), main = "")  # plota o histograma
linhas (x,dgamma(x,shape=alpha,rate=beta)) # Traça a função de destino e é verificado se elas se parecem.

A imagem será apresentada aqui.

c) Resolução computacional.

library(truncnorm)
x=seq(0,10,length.out=n)
random<- rtruncnorm(x, a=2)
prandom<-ptruncnorm(random, a=2)
qrandom<-qtruncnorm(prandom, a=2)
#Comparação do histograma com a função densidade
plot(x,dtruncnorm(x, a = 2),type="l",xlab="x",ylab="Densidade", main = "")
hist(qrandom,probability=T,xlab="x",col=16, breaks = (n/16), main="") 
# Plota histograma
lines(x, dtruncnorm(x, a=2))

A imagem será apresentada aqui.

Conforme pode ser observado nos gráficos apresentados nas letras "a", "b" e "c", existe aderência da função densidade com os seus respectivos histogramas.

...