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

Gere as seguintes variáveis aleatórias, comparando-as com as respectivas funções de densidade: (a) Normal usando amostragem por rejeição com candidatos de uma Cauchy, (b) Gama(4,3, 6,2) usando Gama(4,7), (c) Normal Padrão Truncada em \((2,\infty)\).

0 votos
32 visitas
perguntada Set 29 em Programação Computacional por Thiago Trafane (21 pontos)  
editado Set 29 por Thiago Trafane

Exercício 2.3 do capítulo 2 do livro "Monte Carlo Statistical Methods" de Christian P. Robert e George Casella.

Compartilhe

1 Resposta

0 votos
respondida Set 29 por Thiago Trafane (21 pontos)  
editado Set 29 por Thiago Trafane

1. Discussão itens (a) e (b)

Para resolver os dois primeiros itens, eu usei o método de amostragem por rejeição. Seja \(f\) a função densidade que queremos amostrar, \(g\) a função de densidade que vamos usar para gerar os possíveis candidatos e \(M\) uma constante tal que \(f(x) \leq Mg(x)\) \(\forall x\). Então, como descrito no algoritmo A.4 do capítulo 2 do livro "Monte Carlo Statistical Methods" de Christian P. Robert e George Casella, devemos seguir os seguintes passos para obter uma observação: (1) gere \(X \sim g\) e \(U \sim \mathcal{U}_{[0,1]}\) e (2) aceite \(X\) se \(U \leq f(x)/Mg(x)\) ou retorne para o passo (1) caso contrário.

Assim, o grande desafio de implementar esse método é obter o valor de \(M\). Ou, equivalentemente, sendo \(h(x) \equiv f(x)/g(x)\), o desafio é encontrar o máximo da função \(h\).

No caso do item (a), \(f\) é uma Normal e \(g\) é uma Cauchy. Como o enunciado não especifica os parâmetros dessas distribuições, eu considerei uma Normal Padrão e uma Cauchy Padrão. Nesse caso,
\( h(x) = \frac{\left(\frac{1}{\sqrt{2\pi}}\right)\exp(-x^2/2)}{\left(\frac{1}{\pi}\right)\left(\frac{1}{x^2+1}\right)} = \left(\frac{\pi}{\sqrt{2\pi}}\right)\left(\frac{x^2+1}{\exp(x^2/2)}\right)\)
\( h' (x) = \left(\frac{\pi}{\sqrt{2\pi}}\right)\left(\frac{2x-(x^2+1)x}{\exp(x^2/2)}\right) = \left(\frac{\pi}{\sqrt{2\pi}}\right)\left(\frac{x-x^3}{\exp(x^2/2)}\right)\)
\( h'' (x) = \left(\frac{\pi}{\sqrt{2\pi}}\right)\left(\frac{1-3x^2-(x-x^3)x}{\exp(x^2/2)}\right) = \left(\frac{\pi}{\sqrt{2\pi}}\right)\left(\frac{1-4x^2+x^4}{\exp(x^2/2)}\right)\)

\( h' (x) = 0 \rightarrow x(x^2-1) = 0 \)

Então, os pontos críticos são 1, -1 e 0. É fácil ver que \(h'' (0) = \frac{\pi}{\sqrt{2\pi}} \gt 0\) e que \(h'' (1)= h''(-1)=-\left(\frac{\pi}{\sqrt{2\pi}}\right)\left(\frac{2}{\exp(1/2)}\right) \lt 0 \). Logo, apenas 1 e -1 são máximos locais. De fato, 1 e -1 são máximos globais, já que a função é derivável em \(\Re\) e \(\lim_{x \rightarrow \infty}h (x) = \lim_{x \rightarrow -\infty} h(x) = 0\). Assim, \(M=h(1)=h(-1)\).

No caso do item (b), tanto \(f\) como \(g\) são distribuições Gama, mas com parâmetros distintos. Assim, considere que \(f\) é uma \(Gama(\alpha,\beta)\) e \(g\) é uma \(Gama(\tilde{\alpha},\tilde{\beta})\). Então,
\( h(x) = \frac{\frac{\beta^{-\alpha} x^{\alpha-1}\exp(-x/\beta)}{\Gamma(\alpha)}}{\frac{\tilde{\beta}^{-\tilde{\alpha}} x^{\tilde{\alpha}-1}\exp(-x/\tilde{\beta})}{\Gamma(\tilde{\alpha})}} = \left(\frac{\Gamma(\tilde{\alpha})\tilde{\beta}^{\tilde{\alpha}}}{\Gamma(\alpha)\beta^{\alpha}} \right) \left[ \frac{x^{\alpha-\tilde{\alpha}}}{\exp\left((\beta^{-1}-\tilde{\beta}^{-1})x\right)} \right]\)
\( h' (x) = \left(\frac{\Gamma(\tilde{\alpha})\tilde{\beta}^{\tilde{\alpha}}}{\Gamma(\alpha)\beta^{\alpha}} \right) x^{\alpha-\tilde{\alpha}-1} \left[ \frac{(\alpha-\tilde{\alpha})-(\beta^{-1}-\tilde{\beta}^{-1})x}{\exp\left((\beta^{-1}-\tilde{\beta}^{-1})x\right)} \right] \)
\( h'' (x) = \left(\frac{\Gamma(\tilde{\alpha})\tilde{\beta}^{\tilde{\alpha}}}{\Gamma(\alpha)\beta^{\alpha}} \right) x^{\alpha-\tilde{\alpha}-1} A \)
em que \( A = \frac{(\alpha-\tilde{\alpha})\left((\alpha-\tilde{\alpha}-1)x^{-1}-2(\beta^{-1}-\tilde{\beta}^{-1})\right) + (\beta^{-1}-\tilde{\beta}^{-1})^2 x}{\exp\left((\beta^{-1}-\tilde{\beta}^{-1})x\right)} \).

Note que se \(\alpha \leq \tilde{\alpha}\) e \(\tilde{\beta} \gt \beta\), \(h' (x) \lt 0\) \( \forall x \in \Re^{++}\). Se \(\alpha \leq \tilde{\alpha}\) e \(\tilde{\beta} \lt \beta\), \(\lim_{x\rightarrow \infty} h(x) = \infty\). Logo, se \(\alpha \leq \tilde{\alpha}\), a função \(h\) não tem ponto de máximo em \(\Re^{++}\) para \(\beta \neq \tilde{\beta}\). Para \( \beta = \tilde{\beta}\), \(h' (x) \lt 0\) \( \forall x \in \Re^{++}\) se \(\alpha \lt \tilde{\alpha}\), implicando novamente na ausência de ponto de máximo. Ou seja, se \(\alpha \leq \tilde{\alpha}\), só podemos considerar o caso trivial \(\alpha = \tilde{\alpha}\) e \(\beta = \tilde{\beta}\).

Nesse contexto, considere que \(\alpha \gt \tilde{\alpha}\). Então, é fácil verificar que não existe ponto de máximo se \( \tilde{\beta} \leq \beta\), já que \(\lim_{x\rightarrow \infty} h(x) = \infty\) sob tais condições. Por isso, considere também que \( \tilde{\beta} \gt \beta\).

Então,
\( h' (x) = 0 \rightarrow x = \frac{\alpha-\tilde{\alpha}}{\beta^{-1}-\tilde{\beta}^{-1}} \gt 0 \)

Para verificar a condição de segunda ordem, vamos avaliar \(A\) no ponto crítico:
\( A = \frac{(\alpha-\tilde{\alpha}-1)(\beta^{-1}-\tilde{\beta}^{-1}) - (\beta^{-1}-\tilde{\beta}^{-1})(\alpha-\tilde{\alpha}) }{\exp\left(\alpha-\tilde{\alpha}\right)} = -\frac{\beta^{-1}-\tilde{\beta}^{-1}}{\exp\left(\alpha-\tilde{\alpha}\right)} \)

Logo, como \(\tilde{\beta} \gt \beta\), \(A \lt 0\) e, consequentemente, a derivada segunda no ponto crítico é negativa. Assim, o ponto crítico é ponto de máximo local. É, de fato, ponto de máximo global, dado que \(h\) é derivável em \(\Re^{++}\), \(h' (x) \lt 0\) para \(x \gt \frac{\alpha-\tilde{\alpha}}{\beta^{-1}-\tilde{\beta}^{-1}}\) e \(h' (x) \gt 0\) para \(0 \lt x \lt \frac{\alpha-\tilde{\alpha}}{\beta^{-1}-\tilde{\beta}^{-1}}\). Assim, \(M = h \left( \frac{\alpha-\tilde{\alpha}}{\beta^{-1}-\tilde{\beta}^{-1}} \right)\).

Por fim, vale ressaltar que o caso proposto no item (b), em que \(\alpha=4,3\), \(\beta = 6,2\), \(\tilde{\alpha} = 4\) e \(\tilde{\beta}=7\), satisfaz as condições \(\alpha \gt \tilde{\alpha}\) e \( \tilde{\beta} \gt \beta\) e, logo, o exercício pode ser resolvido usando a estratégia aqui proposta.

2. Discussão item (c)

A solução do item (c) é mais direta. Eu simplesmente extraí observações de uma Normal Padrão, descartando as observações fora do suporte.

3. Solução dos problemas

O código que resolve os três itens dessa questão é apresentado abaixo. Como se vê, eu criei cinco funções. As duas primeiras são apenas funções auxiliares que geram gráficos e as três últimas são as funções que efetivamente geram as amostras das distribuições. O output principal dessas três últimas funções é um gráfico comparando o histograma das observações geradas com a função densidade desejada (a da Normal Truncada é a densidade da Normal no suporte divida pela probabilidade da variável aleatória estar no suporte, isto é, a densidade da Normal no suporte ajustada para integrar 1). No caso das duas primeiras funções, que são usadas para resolver os itens (a) e (b), o programa também gera um gráfico para revelar possíveis inconsistências no cálculo de \(M\), que foi discutido na seção 1.

# 0) Importando pacotes
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt


# 1) Definindo as funções
# 1.1) Função auxiliar que constrói os gráficos mostrando a adequação do M obtido
def plot_M(x_grid,M_grid,x_opt,M,title,file_name):       
    plt.plot(x_grid,M_grid,label='Razão')

    tam_x_opt = len(x_opt)        
    for i in range(tam_x_opt):
        if i == 0:        
            plt.plot(x_opt[i],M,'ro',label='M')
        else:
            plt.plot(x_opt[i],M,'ro')

    plt.title(title)
    plt.legend()
    plt.savefig(file_name)
    plt.show()   

# 1.2) Função auxiliar que constrói os gráficos dos histogramas
def plot_hist(x_pdf,pdf,draws,title,file_name):
    plt.plot(x_pdf,pdf,'r-', lw=2, label='PDF')    
    plt.hist(draws,bins='auto',density=True,label='Histograma')
    plt.title(title)
    plt.legend()
    plt.savefig(file_name)
    plt.show()    

# 1.3) Função que gera draws de uma distribuição Normal Padrão usando uma Cauchy Padrão
def st_normal_draws(n):      
    # a) Obtendo M
    x_opt = [-1,1]
    M = stats.norm.pdf(x_opt[0],loc=0,scale=1)/stats.cauchy.pdf(x_opt[0],loc=0,scale=1)

    # b) Construindo gráfico mostrando que estamos usando o M adequado
    x_grid = np.linspace(0-2,0+2,100)
    M_grid = stats.norm.pdf(x_grid,loc=0,scale=1)/stats.cauchy.pdf(x_grid,loc=0,scale=1)
    titulo_gM = 'Razão N(0,1)/Cauchy(0,1) e M'
    plot_M(x_grid,M_grid,x_opt,M,titulo_gM,'gM_st_norm')  

    # c) Obtendo draws
    n_obs = 0
    valid_draws = []            
    while n_obs < n:              
        draw = np.random.standard_cauchy() #Draw de teste vem de uma distribuição Cauchy Padrão             
        aceit_crit = np.random.uniform(low=0,high=1) #Critério de aceitação vem de uma distribuição Uniforme entre 0 e 1

        if aceit_crit <= stats.norm.pdf(draw,loc=0,scale=1)/(M*stats.cauchy.pdf(draw,loc=0,scale=1)): #Método de amostragem por rejeição
            valid_draws.append(draw)
            n_obs = n_obs+1

    valid_draws = np.array(valid_draws) #Ajustando formato

    # d) Construindo gráfico draws    
    x_pdf = np.linspace(min(valid_draws),max(valid_draws), 100)
    pdf = stats.norm.pdf(x_pdf,loc=0,scale=1)  
    titulo_gH = 'Histograma Normal(0,1)'
    plot_hist(x_pdf,pdf,valid_draws,titulo_gH,'gH_st_norm')

    # e) Definindo output    
    return valid_draws

# 1.4) Função que gera draws de uma distribuição Gama usando outra Gama
def gamma_draws(alpha,beta,alpha_draw,beta_draw,n):     
    # 1.4.1) Caso 1: caso trivial
    if alpha==alpha_draw or beta==beta_draw:
        valid_draws = np.random.gamma(shape=alpha_draw,scale=beta_draw,size=n)
        return valid_draws    

    # 1.4.2) Caso 2: caso não válido
    elif (alpha<=alpha_draw or beta>=beta_draw):
        print('Erro: condições do alfa e beta não são satisfeitas')
        return

    # 1.4.3) Caso 3: caso principal
    else:
        # a) Obtendo M        
        x_opt = [(alpha-alpha_draw)/((1/beta)-(1/beta_draw))]
        M = stats.gamma.pdf(x_opt[0],a=alpha,loc=0,scale=beta)/stats.gamma.pdf(x_opt[0],a=alpha_draw,loc=0,scale=beta_draw)

        # b) Construindo gráfico mostrando que estamos usando o M adequado  
        x_grid = np.linspace(0.000001,alpha*beta+2*alpha*beta**2,100)
        M_grid = stats.gamma.pdf(x_grid,a=alpha,loc=0,scale=beta)/stats.gamma.pdf(x_grid,a=alpha_draw,loc=0,scale=beta_draw)
        titulo_gM = 'Razão Gama('+str(alpha)+','+str(beta)+')/Gama('+str(alpha_draw)+','+str(beta_draw)+') e M'   
        plot_M(x_grid,M_grid,x_opt,M,titulo_gM,'gM_gama')  

        # c) Obtendo draws
        n_obs = 0
        valid_draws = []     
        while n_obs < n:              
            draw = np.random.gamma(shape=alpha_draw,scale=beta_draw) #Draw de teste vem de uma distribuição Gama com parâmetros alpha_draw e beta_draw        
            aceit_crit = np.random.uniform(low=0,high=1) #Critério de aceitação vem de uma distribuição Uniforme entre 0 e 1

            if aceit_crit <= stats.gamma.pdf(draw,a=alpha,loc=0,scale=beta)/(M*stats.gamma.pdf(draw,a=alpha_draw,loc=0,scale=beta_draw)): #Método de amostragem por rejeição
                valid_draws.append(draw)
                n_obs = n_obs+1

        valid_draws = np.array(valid_draws) #Ajustando formato     

        # d) Construindo gráfico draws      
        x_pdf = np.linspace(min(valid_draws),max(valid_draws), 100)
        pdf = stats.gamma.pdf(x_pdf,a=alpha,loc=0,scale=beta)
        titulo = "Histograma Gama("+str(alpha)+','+str(beta)+')'
        plot_hist(x_pdf,pdf,valid_draws,titulo,'gH_gama')

        # e) Definindo output
        return valid_draws

# 1.5) Função que gera draws de uma distribuição Normal Padrão Truncada
def st_normal_trunc_draws(lim_inf,lim_sup,n):
    # a) Obtendo draws
    n_obs = 0
    valid_draws = []

    while n_obs < n:
        draw = np.random.normal(loc=0,scale=1) #Draw de teste vem de uma distribuição Normal Padrão

        if draw>lim_inf and draw<lim_sup: #Draw é descartado se estiver fora do suporte da Normal Padrão Truncada
            valid_draws.append(draw)
            n_obs = n_obs+1

    valid_draws = np.array(valid_draws) #Ajustando formato  

    # b) Construindo gráfico draws    
    x_pdf = np.linspace(min(valid_draws),max(valid_draws), 100)  
    pdf = stats.norm.pdf(x_pdf,loc=0,scale=1)/(stats.norm.cdf(lim_sup,loc=0,scale=1)-stats.norm.cdf(lim_inf,loc=0,scale=1))
    titulo_gH = 'Histograma Normal Padrão Truncada em ('+str(lim_inf)+','+str(lim_sup)+')'
    plot_hist(x_pdf,pdf,valid_draws,titulo_gH,'gH_st_norm_trunc')      

    # c) Definindo output    
    return valid_draws


# 2) Gerando os resultados solicitados no exercício 2.3
if __name__ == '__main__':        
    n = 100000 #Tamanho da amostra

    draws_st_normal = st_normal_draws(n) #Item (a)
    draws_st_gamma = gamma_draws(4.3,6.2,4,7,n) #Item (b)
    draws_st_normal_trunc = st_normal_trunc_draws(2,float('inf'),n) #Item (c)    

Os gráficos gerados considerando 100.000 observações (válidas) para cada uma das três distribuições são mostrados abaixo. Como se vê, em cada um dos três casos, os histogramas se mostram consistentes com as funções de densidade. Ademais, os valores de \(M\) considerados nos itens (a) e (b) parecem ser de fato os pontos de máximo das respectivas funções \(h\), como seria esperado.

A imagem será apresentada aqui.
A imagem será apresentada aqui.
A imagem será apresentada aqui.

A imagem será apresentada aqui.
A imagem será apresentada aqui.

...