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

Uma integral é aproximada através do Método de Monte Carlo. Como encontrar um valor \(\delta\) tal que a probabilidade de o erro da estimativa ser menor ou igual a \(\delta\) é de 5%? Na solução, eu devo utilizar o Teorema do Limite Central.

0 votos
74 visitas
perguntada Abr 12 em Estatística por Ricardo Saldanha (1 ponto)  
editado Abr 19 por Ricardo Saldanha

[Adaptado de RICE, J, A. Mathematical Statistics and Data Analysis (2006) - cap. 5, exc. 22]
Seja \(\hat{I}(f)\) a aproximação de \(I(f) = \int_{0}^{1} cos(2\pi x) dx \) obtida através do Método de Monte Carlo, com 1000 termos. Use o Teorema do Limite Central para encontrar um \(\delta\) tal que \(P(|\hat{I}(f) - I(f)| \leq \delta)=0.05\).

Compartilhe

1 Resposta

0 votos
respondida Abr 12 por Ricardo Saldanha (1 ponto)  
editado Mai 22 por Ricardo Saldanha

Caso você não esteja familiarizado com o Método de Monte Carlo, veja Rice - Exemplo A, capítulo 5, o mesmo do exercício. Em suma, \({\displaystyle \hat{I}(f) = \frac{1}{N} \sum_{n=1}^{N} f(x_{n}) }\) , em que \(\{x_{n}\}\) são \(N\) variáveis aleatórias uniformemente distribuídas no intervalo em que a integral é avaliada. Decorrerá da Lei dos Grandes Números e da definição de \(E(f(x_{n}))\) (i.e.. a definição de esperança de uma função sobre uma variável aleatória) que \(\hat{I}(f)\) aproxima \(I(f)\), sendo mais precisa a estimativa a medida que \(N\) aumenta. Neste caso, \({\displaystyle \hat{I}(f) = \frac{1}{1000} \sum_{n=1}^{1000} cos(2\pi x_{n}) }\).

Antes de passarmos para a solução, apresentemos uma versão do Teorema do Limite Central e dois fatos:

[Adaptado de STACHURSKI, J. (2016). A Primer in Econometric Theory - Theorem 6.1.2]:
Teorema do Limite Central (CLT) Seja um \(x\) com segundo momento finito e \(\{{x_n}\}\) \(N\) cópias i.i.d. de \(x\). Então
\( \kern{5em} \sqrt{N}\frac{(\bar{x}-\mu)}{\sigma} \xrightarrow{d} N(0,1)\) quando \(N \rightarrow \infty \)
, em que \(\bar{x}\) é a média dos \(\{x_n\}\), \(\mu=E(x)\) e \(\sigma^{2}=var(x)\).

[IDEM - Fact 5.1.10]:
Fato 1 Se \(x_{1},...,x_{N}\) são variáveis aleatórias independentes e \(f_{1},...,f_{N}\) são funções Borel-mensuráveis, então \(f_{1}(x_{1}),...,f_{N}(x_{N})\) também são independentes.

Embora a demonstração desse fato seja bem curta, evidentemente exige o conhecimento do que é uma função Borel-Mensurável (vd. Stachurski, cap. 4). Para nosso propósito, a intuição é suficiente: se, por exemplo, \(x_{1}\) não contém informação sobre \(x_{2}\), então \(f_{1}(x_{1})\) não contém informação sobre \(f_{2}(x_{2})\). Pense como isso faz sentido no nosso caso. Sortear um \(x_{1}\) em \([0,1]\) e avaliar \(cos(2\pi x_{1})\) não nos dá informação sobre qual seria \(cos(2\pi x_{2})\), já que \(x_{2}\) também será sorteado.

[Da resolução da questão 20, do mesmo capítulo de Rice]
Fato 2 \(var(\hat{I}(f)) = \frac{1}{N}( I(f^{2})-(I(f))^2) = \frac{1}{2} \)

A expressão da variância deve-se à definição de \(E(f(x_{n}))\) e ao fato de \(\{x_{n}\}\) serem i.i.d. (o que, pelo Fato 1 acima e pelas distribuições idênticas, leva a que as \(var(f(x_{n}))\) são iguais para qualquer \(n\)). O resultado \(\frac{1}{2}\) é a avaliação da expressão dada a mesma integral \(I(f)\) deste exercício.

Estamos prontos para a solução propriamente dita:

Pelo Fato 1, \(f(x_{1}),...,f(x_{N})\) são independentes. E, como a função \(f\) é sempre a mesma, claramente também serão identicamente distribuídos. Além disso, \(\hat{I}(f)\) é a média dos \(\{f(x_{n})\}\). Dessa maneira, no lugar de \(\bar{x}\) e o parâmetro \(\mu\) do qual o mesmo se aproxima, podemos aplicar o CLT a \(\hat{I}(f)\) e \(I(f)\). Uma aproximação de \(1000\) termos é grande o bastante para considerarmos que:
\( \hat{Z}(f) \mathrel{\mathop:}= \sqrt{N}\frac{\hat{I}(f)-I(f)}{dp(\hat{I}(f))} \approx N(0,1) \)

Valendo-se do método de substituição, é fácil encontrar a solução exata da integral: \(I(f)=0\). Dessa forma, temos:
\(P(|\hat{I}(f) - I(f)| \leq \delta)=0.05\)
\(P(-\delta \leq \hat{I}(f) \leq \delta)=0.05\)
\( {\displaystyle P\bigg(\sqrt{N}\frac{-\delta}{dp(\hat{I}(f))} \leq \hat{Z}(f) \leq \sqrt{N}\frac{\delta}{dp(\hat{I}(f))}\bigg)=0.05 }\)

Os valores críticos da distribuição amostral de \( \hat{Z}(f)\) serão associados a uma região parecida com a faixa roxa ilustrada abaixo. A faixa será bem estreita porque a probabilidade é de apenas \(5\%\), assim, esperamos encontrar valores críticos bem próximos de \(0\).
A imagem será apresentada aqui.

Pela tabela de distribuição normal acumulada no apêndice B de Rice, o valor crítico \( Z^{*} = \sqrt{N}\frac{\delta}{dp(\hat{I}(f))} \) mais próximo do que seria correspondente à probabilidade \(0.525\) (\(0.5+0.0250\)) é \(Z^{*} = 0.06\) (associado na tab. à prob. \(0.5239\)).

E, conforme Fato 2, o \(dp(\hat{I}(f)) = \sqrt{0.5} \). Logo, concluímos que:
\( \delta = 0.06 \frac{ \sqrt{0.5}}{ \sqrt{1000}} \approx 0.001341641\)

Em complemento à resposta da questão. O código abaixo testa para várias simulações a proporção que fica no intervalo (-0.001341641,0.001341641). Se os cálculos estiverem corretos, devemos encontrar proporção próxima de 0.05.

import numpy as np
import random as rd
import scipy.stats as st

def simulacao(N):
    soma = 0
    for i in range(N):
        x = rd.uniform(0,1)
        y = np.cos(2*x*np.pi)
        soma += y
    media = soma/N
    return media

def within_delta(num_reps=100,N=1000,delta=0.00134):
    print('%s simulações, cada uma com N = %s.' %(num_reps,N))
    ''' Queremos contar quantas das médias simuladas tem módulo menor que 
    delta; i.e.,  distanciam-se menos do que delta da solução exata.'''
    count_cond = 0
    for i in range(1,num_reps+1):
        media = simulacao(N)
        if abs(media)<delta:
            count_cond += 1
    # Proporção de simulações que atenderam à condição
    prop_cond = count_cond/num_reps
    print('Propração de simulações no intervalo +-%10.9f: %10.9f' %(delta,prop_cond))

if __name__ == '__main__':

    num_reps = 10000
    N = 1000

    # Testando o delta encontrado a partir tabela de distb..
    within_delta(num_reps,N,0.001341641)
    # Usando delta mais preciso.
    delta = st.norm.ppf(.525)*0.5**0.5/N**0.5
    within_delta(num_reps,N,delta) 

A minha execução do código retornou:
"10000 simulações, cada uma com N = 1000.
Propração de simulações no intervalo +-0.001341641: 0.048400000
10000 simulações, cada uma com N = 1000.
Propração de simulações no intervalo +-0.001402166: 0.049800000"
Conforme pode ser visto na rotina, esse segundo intervalo é calculado pela biblioteca spipy.stats, no lugar da busca manual numa tabela.

Os resultados de fato são próximos a 0.05.

...