Uma comparação entre técnicas de propagação de erros em Astrofísica: Monte Carlo × Bootstrap Alexandre Zabot, Raymundo Baptista Departamento de Física, Universidade Federal de Santa Catarina, Campus Trindade, 88040-900, Florianópolis-SC, Brazil
Resumo Neste trabalho é feito um estudo comparativo entre dois algoritmos numéricos usados para propagação de erros em dados experimentais. Um deles é conhecido por Método de Monte Carlo e o outro por Método de Bootstrap. Recentemente, Dhillon & Watson argüiram que a aplicação do método de Monte Carlo introduz ruído nos dados, e propuseram então a utilização do Bootstrap como alternativa capaz de produzir resultados superiores. O objetivo deste trabalho é testar a validade dessa afirmação. As duas técnicas foram aplicadas a três problemas diferentes: o ajuste de modelos de emissao LTE simples e atmosfera estelar a espectros estelares observados e o ajuste de curvas de luz de eclipses de Variáveis Cataclísmicas para a determinação da distribuição radial de brilho dos seus discos de acréscimo. Os métodos foram testados quanto à sua robusteza, ou seja, a capacidade de prover resultados coerentes entre si. Além disso, as soluções dos métodos foram comparadas. Os resultados indicam que não existe evidência de superioridade de um método em relação ao outro. Palavras Chave: Propagação de erros, Monte Carlo, Bootstrap
1
Introdução
Quando uma grandeza qualquer é obtida a partir de um conjunto de dados experimentais, inevitavelmente surge a questão da propagação de erros. Uma vez que os dados são intrinsecamente afetados por erros, o que se deseja saber Email:
[email protected] (Alexandre Zabot),
[email protected] (Raymundo Baptista).
é qual a incerteza associada à determinação da grandeza em questão. A propagação de erros em dados experimentais é um problema constante em qualquer área da física experimental e também em Astrofísica [9]. Saber como os erros nos dados se propagam para a grandeza calculada requer uma análise estatística. Quando a determinação da grandeza em questão a partir dos dados experimentais é feita através de uma função, pode-se usar uma expressão analítica para encontrar o erro propagado. A fórmula é inspirada no truncamento em primeira ordem da expansão em série de Taylor para funções de várias variáveis, mais o pressuposto estatístico de que os erros se somam quadraticamente, 2
(∆y) =
" n → X ∂y(− x) i=1
∂xi
∆xi
#2
(1)
,
onde xi é o valor do dado experimental, ∆xi sua incerteza e ∆y é a incerteza na grandeza final obtida, y. Entretando, nem sempre a grandeza a ser calculada a partir dos dados pode ser expressa de forma analítica. Nestas situações, é preciso encontrar outra maneira de avaliar o erro propagado. Um procedimento numérico comumente usado consiste em realizar simulações utilizando o método de Monte Carlo (MC). Entretanto, como será visto na próxima seção, o método de MC requer que o valor dos pontos experimentais seja alterado. Isso levou Dhillon & Watson [3] a afirmar que a aplicação do MC superestima o erro propagado pois insere ruído adicional nos dados. Como alternativa ao problema, eles propuseram a utilização do método de Bootstrap (BS) [4]. Segundo eles, o BS seria capaz de produzir resultados superiores, uma vez que esta técnica não altera os valores dos dados.
2
Os métodos de análise
A implementação do método de MC é muito simples. Toma-se o arquivo de dados original e multiplica-se o valor da incerteza em cada dado por um número aleatório. Esse número é gerado por uma função com distribuição gaussiana de média zero e variância unitária. Na nossa implementação, usamos a função gasdev do Numerical Recipes in C [10]. A seguir, soma-se o resultado ao valor original do dado. Assim, o valor do dado experimental é alterado aleatoriamente de acordo com uma distribuição normal com desvio padrão igual à sua incerteza. É como se tivesse sido feita uma outra medida ( Fig. 1 ). Este procedimento é repetido várias vezes, para simular várias medidas. Depois, os dados experimentais sintéticos são processados pelos mesmos métodos numéricos usados para obter a grandeza final desejada. O desvio padrão da grandeza 2
final é tomado como o erro propagado pelo método numérico. Sendo uma alternativa ao MC, o método de BS é baseado no princípio de não alterar o valor do dado experimental. Como a obtenção da grandeza desejada leva em conta os erros nos dados, para calcular o erro propagado pode-se alterar o valor da barra de erro ao invés de alterar os dados experimentais em si. Ou seja, o BS também simula uma outra medida, sendo que o que muda não é o valor de cada dado, mas sim o seu peso ( i.e., a incerteza associada ). Assim como no MC, repete-se este procedimento várias vezes, para simular várias medidas. Depois, os dados experimentais sintéticos são processados pelos mesmos métodos numéricos para obter a grandeza final desejada. Novamente, o desvio padrão do valor final é considerado uma estimativa do erro propagado. Porém, como mexer nas barras de erros? Suponha que haja n dados experimentais. Sorteia-se n vezes um número entre 1 e n. Assim se houverem 5 dados experimentais, sorteia-se 5 vezes um número entre 1 e 5. A seguir, conta-se quantas vezes cada número apareceu. Se, por exemplo, num dado sorteio o número 3 apareceu √ 2 vezes, então divide-se a barra de erro do terceiro dado experimental por 2. Se, por outro lado, o número 1 não foi sorteado, multiplica-se a barra de erro do primeiro dado experimental por um número grande, de modo que este ponto torne-se desprezível para o ajuste ( por exemplo, multiplica-se o erro por 10 ). De maneira geral, podemos dizer que se o número i apareceu j vezes, divide-se a barra de erro do i-ésimo dado experi√ mental por j caso j seja maior que zero. Se j for zero, multiplica-se a barra de erro por 10. Na Fig. 1 é mostrado um exemplo da aplicação dos métodos de MC e do BS a uma situação que simula dados experimentais.
3
Procedimentos e Resultados
Para avaliar a eficiência e a robustez dos dois métodos, eles foram aplicados a três problemas diferentes. Primeiro, estudou-se o erro propagado na temperatura e no fator de escala (eq. 2) quando se ajusta modelos de atmosferas ( emissão simples LTE de hidrogênio na subseção 3.1 e modelos de atmosferas estelares na subseção 3.2 ) a espectros sintéticos com ruído ( visando simular dados reais ). Por fim, na subseção 3.3, avaliou-se o erro propagado para o valor da taxa de acréscimo para um disco em estado estacionário, a partir do ajuste de curvas de luz em eclipse usando técnicas de imageamento indireto. Definimos o fator de escala através da seguinte equação: A≡
y , m
(2)
onde A é o fator de escala, y é o fluxo médio sobre todos os comprimentos de 3
Figura 1. Em cima: Um conjunto de dados gerados a partir de uma distribuição senoidal cujas barras de erro representam 10% da ordenada do ponto. No meio: Dados alterados de acordo com o método de MC. Note que os valores dos pontos mudaram mas o tamanho das respectivas barras de erro se manteve constante. Em baixo: Dados alterados de acordo com o método de BS. Neste caso os valores dos pontos permanecem constantes. Entretanto, as barras de erros foram alteradas conforme o procedimento de sorteio discutido na seção 2.
onda do espectro sintético e m é o fluxo médio dos modelos que estão sendo ajustados ao espectro sintético. Para todos os casos foram feitos ajustes a conjuntos de dados de entrada originais para 8 diferentes razões sinal/ruído ( S/R = 10, 13, 20, 30, 40, 50, 100 e 200 ). O objetivo de usar várias razões S/R é verificar qual a dependência do erro propagado e da sua incerteza com a razão S/R. 3.1 Ajuste Espectral: Atmosfera de Hidrogênio Neste problema gerou-se um espectro sintético de hidrogênio com temperatura T= 104 K e Densidade de Coluna Σ = 1018 barions/cm2 usando a ferramenta 4
calcspec do pacote iraf 1 . Adotou-se esse espectro como sendo o espectro original sem erros. Em seguida, utilizando uma rotina geradora de número com distribuição normal [10], inseriu-se ruído artificialmente nos dados de modo a produzir espectros sintéticos para as razões S/R mencionadas anteriormente. Deste modo, foram criados 8 espectros sintéticos com ruído. Para cada um desses novos espectros aplicou-se os métodos de MC e de BS, gerando 100 espectros sintéticos para cada S/R e para cada método. Ao todo foram criados 1600 espectros. Para fazer o ajuste, foi criado um grid de modelos de espectros de hidrogênio. Os espectros deste grid têm a mesma densidade de coluna que os espectros anteriores. Entretanto, varrem uma faixa de temperaturas que vai de 5000 K até 20000 K, em intervalos de 500 K. Para ajustar cada um dos 1600 espectros criados pelos métodos, foi criado um programa de ajuste. Inicialmente varre-se todo o espaço de parâmetros ( temperatura e fator de escala ) em busca do par de valores de fator de escala e temperatura que produz um mínimo da função χ2 , 2
χ =
N X i=1
Ã
fimod − fiobs σiobs
!2
(3)
,
onde N é o número de dados, fimod o fluxo do espectro modelo, fiobs o fluxo do espectro observado e σiobs a incerteza no fluxo observado, que é igual ao fluxo dividido pelo S/R. Uma vez localizada a região de mínimo global, o programa chama a rotina amoeba do Numerical Recipes [10] para refinar a posição do mínimo. A Fig. 2 mostra o exemplo de um espectro sintético com S/R=10 e o espectro de melhor ajuste encontrado. Desta maneira, para cada método e para cada S/R, foram obtidos 100 valores de melhor ajuste para a temperatura (que coincidem com o valor inicial de 10.000 K ao nivel de 1-σ) e para o fator de escala. Calculou-se então o valor médio e o desvio padrão para essas grandezas. Adotou-se o desvio padrão como sendo o erro propagado. Para avaliar a robustez e a consistência dos resultados, repetiu-se todo o procedimento 10 vezes, gerando 10 valores para o erro propagado na temperatura e no fator de escala, para cada método e para cada S/R. Com o intuito de testar se os métodos fornecem sempre os mesmos valores para o erro propagado, calculou-se a média e o desvio padrão desses 10 valores de erro propagado. Os resultados podem ser vistos na Fig. 3. São mostrados o erro relativo percentual propagado contra o inverso da S/R ( vezes 100 ). Ou seja, mostramos o valor médio do erro propagado dividido pelo valor médio de melhor ajuste contra a incerteza dos dados originais. 1
iraf é distribuído pelos National Optical Astronomy Observatories, que é operado pela Association of Universities for Research in Astronomy, Inc., sob contrado com a National Science Foundation.
5
Figura 2. Espectro dos dados sintéticos ( pontos ) gerado com o método de BS para uma razão S/R=10 a partir do espectro de atmosfera de Hidrogênio. A linha cheia mostra o espectro de melhor ajuste.
Pode-se ver que, à medida que os dados de entrada tem razões S/R mais baixas, o valor do erro propagado aumenta, o que é algo natural de se esperar. Nota-se também que a variação no erro percentual dos pontos é pequena para altas S/R, mas aumenta progressivamente à medida que a qualidade dos dados de entrada é deteriorada. Isso mostra que os métodos são robustos apresentando comportamento similar em função da razão S/R. No entanto, no limite em que os dados de entrada são muito ruins, as respostas podem oscilar bastante. Nota-se ainda que os dois métodos produzem resultados semelhantes; os erros relativos em cada parâmetro são comparáveis para uma dada S/R. Ou seja, não existe evidência de que o MC superestima o erro propagado neste exemplo. Se ele faz isto, o BS o faz da mesma forma. Entretanto, o MC parece ser uma técnica mais robusta porque fornece barras de erro (que refletem a variação dos valores de erros propagados) menores que as do BS. Esse efeito será mais evidenciado nos resultados da subseção 3.2, onde as barras de erro associadas ao BS são sempre maiores que as do MC. Esse efeito está provavelmente relacionado com a estratégia do BS de usar um número reduzido de pontos, uma vez que o método induz o programa de ajuste a descartar alguns pontos que tem grandes barras de erro. O ajuste do modelo parece ser mais sensivel à escolha particular do subconjunto de dados usado no ajuste do que às oscilações nos dados em torno de seus valores médios. Avaliando por este aspecto, podemos afirmar que o método de MC é superior 6
ao método de BS. A similaridade entre os gráficos do erro relativo no fator de escala e na temperatura é digna de nota. Isso não é obra do acaso. Ao contrário, este é um resultado esperado pois existe uma correlação entre a temperatura e o fator de escala. Através da eq. 2, nota-se que maiores (menores) temperaturas indicam fluxos maiores (menores), o que implica em menores (maiores) fatores de escala. Logo, com uma razão S/R mais baixa (alta), tanto a temperatura quanto o fator de escala são menos (mais) bem determinados, uma vez que estão amarrados pela eq. 2.
3.2 Ajuste Espectral: Espectro Estelar
Um problema similiar ao anterior é o ajuste de um modelo de atmosfera estelar a um espectro estelar observado. O espectro utilizado foi escolhido dentre os disponíveis no atlas espectral de Bruzual-Persson-Gunn-Stryker [5] e corresponde à estrela BD 010935. Esta é uma estrela da seqüência principal de tipo espectral B1V. O grid de modelos foi composto por espectros estelares do atlas de Kurucz [8] com metalicidade solar e log g = 4.0 − 4.5. Eles cobrem uma faixa (não igualmente espaçada) de temperaturas que vai de 5500 K a 50000 K. Ajustamos cada modelo do grid ao espectro sintético e, assim como no caso da atmosfera de hidrogênio, foram recuperados valores para a temperatura ( cujo valor médio é algo em torno de 26.000 K ) e o fator de escala de melhor ajuste. Os procedimentos são os mesmos aplicados ao problema da seção anterior. Na Fig. 4 está mostrado um espectro sintético criado pelo método de MC a partir do espectro estelar observado para uma razão S/R=200, mais o modelo de atmosfera de melhor ajuste. A figura 5 mostra o erro propagado para o fator de escala e para a temperatura. Assim como no caso anterior, à medida que a qualidade inicial dos dados é deteriorada, os erros propagados são maiores e as respostas dos métodos oscilam mais. Neste caso também não é possível afirmar que o método de MC superestima o erro propagado, pois os resultados dos dois métodos se superpõem. Isto é, para uma dada S/R, os métodos de MC e BS resultam em erros propagados comparáveis entre si. Nota-se que para este caso, assim como na subseção anterior (3.1), os BS tem barras de erro maiores que o MC. Além disso, para altas razões S/R, o BS produz erros propagados sistematicamente maiores que o MC, o que está em contradição com a sugestão de Dhillon & Watson. 7
Figura 3. Acima: a razão entre o valor médio do erro propagado no fator de escala e o valor médio do fator de escala. Os círculos preenchidos representam o método de MC e os quadrados abertos, o de BS. Abaixo: a mesma coisa para a temperatura. Em ambos os casos, os pontos que representam o método de MC foram deslocados para a direita a fim de facilitar a identificação das barras de erro.
3.3 Mapeamento por eclipse
Variáveis Cataclísmicas são sistemas binários onde uma estrela de tipo espectral tardio ( a secundária ) preenche o seu Lobo de Roche e transfere massa para uma companheira anã branca ( a primária ) via um disco de acréscimo. Em sistemas onde a taxa de acréscimo de matéria é alta, o disco domina a emissão no ótico e no ultra-violeta [11]. Se o observador estiver num ângulo 8
Figura 4. Espectro dos dados sintéticos ( pontos ) gerado com o método de MC para uma razão S/R=200 a partir do espectro de atmosfera estelar observado. A linha cheia mostra o espectro de melhor ajuste.
adequado, ou seja, o mais próximo possível do plano da órbita, então haverá eclipses à medida que uma estrela orbitar a outra. Usando técnicas de mapeamento por eclipse [1,6] é possível reconstruir a distribuição bidimensional de intensidade do disco de acréscimo baseando-se na informação unidimensional da curva de luz do eclipse. Nesta seção utilizamos curvas de luz sintéticas para testar a propagação de erros em experimentos de mapeamento por eclipses usando os métodos de MC e BS. Para as simulações foi alterado um programa já existente, que cria curvas de luz sintéticas a partir da geometria do sistema binário e de uma distribuição de brilho selecionada. Adotamos uma distribuição de brilho para um modelo de disco estacionário e oticamente espesso, no qual a temperatura efetiva varia com a distância R ao centro do disco seguindo a lei T ∝ R−3/4 . Um dos ˙ parâmetros de entrada neste modelo é a taxa de acréscimo de massa (M). ˙ O programa alterado sorteia um novo valor de M para cada fase ao longo ˙ e amplitude da curva de luz usando uma distribuição gaussiana de média M σ ˙ e produz uma nova distribuição de brilho para cada fase da curva de luz. M Curvas de luz sintéticas são então geradas simulando o eclipse da distribuição ˙ variável. Desta forma é inserido ruído artificial nas curvas de brilho com M de luz, controlado pelo parâmetro σ ˙ . Este procedimento foi utilizado para M gerar curvas de luz para o mesmo conjunto de S/R das seções anteriores. 9
Figura 5. Acima: a razão entre o valor médio do erro propagado no fator de escala e o valor médio do fator de escala. Abaixo: a mesma coisa para a temperatura. A notação é a mesma da Fig. 3.
O conjunto de curvas gerado foi utilizado para testar a habilidade dos dois ˙ Para isso ajustou-se métodos em recuperar corretamente a incerteza em M. cada uma das curvas de luz usando o programa prida [2] para produzir o mapa da distribuição de brilho do disco de acréscimo que melhor descreve a curva de luz de entrada. Em seguida, foi obtido o perfil radial de temperatura ˙ pela comparação com de brilho. A partir dele é possível calcular o valor de M modelos de distribuição radial de temperatura de discos estacionários ( e.g. [7] ). Em analogia ao descrito nas seções 3.1 e 3.2, o procedimento é repetido várias 10
vezes de maneira que, para cada S/R são obtidos conjuntos de mapas sintéticos de distribuição de brilho. Os mapas são combinados de modo a produzir, para cada pixel, uma intensidade média e um desvio padrão com relação à ˙ são obtidos das distribuições radiais de temperatura média. Valores de M ˙ e correspondentes de brilho, resultando igualmente em valores médios de M desvios padrão para cada S/R testado. ˙ podemos testar a Neste caso, como nós controlamos o erro inserido em M, capacidade dos métodos utilizados para estimar a barra de erro em recuperar ˙ Se os métodos de propagação de erro o valor correto do erro inserido em M. funcionarem corretamente devemos ter uma reta y = x nos nossos gráficos de erro propagado contra o erro inserido. É importante ressaltar que o teste que estamos aplicando agora é mais completo do que os das seções anteriores (3.2 e 3.1), pois naqueles casos nós não podemos saber qual é o erro propagado que deveria ser retornado, o que não acontece agora. Assim, ao contrário de nos limitarmos à simples comparação entre os dois métodos, agora somos capazes de afirmar se os métodos funcionam ou não ( pelo menos para esse caso ). ˙ em função do erro inserido. Nota-se A Fig. 6 mostra o erro propagado em M que, assim como nos dois casos anteriores, há uma região de intersecção ( ao nível de 1σ ) dos resultados dos dois métodos. No entanto, é possível afirmar que neste problema em particular o MC é estatisticamente mais confiável do que o BS, uma vez que o primeiro fica bem mais perto da relação esperada y = x. Ao contrário dos dois experimentos anteriores, não é possível afirmar que o BS apresenta barras de erro maiores que o MC. Numa análise descuidada, seria possível, inclusive, afirmar o contrário. No entanto, como os dois métodos apresentam valores médios diferentes para o erro, é preciso fazer uma análise considerando-se a razão barra de erro por valor médio. Neste caso, percebe-se que os dois métodos são consistentes entre si, apresentando uma estabilidade comparável. Cabe também notar que, em contraste com o arguido por Dhillon & Watson [3], o método de MC resulta em erros propagados ligeiramente menores que o BS no limite das altas razões S/R. Efeito similar pode ser visto no caso do ajuste espectral (Fig. 5). Em altas razões S/R, as pequenas variações nos dados produzidos pelo MC tornam-se menos importantes do que o aumento nas barras de erro efetuado pelo BS. Este aumento equivale a diminuir o número de dados disponíveis para o ajuste, tornando a determinação da grandeza desejada mais incerta. Por outro lado, no limite das baixas razões S/R ( . 20 ), ˙ enquanto o MC ainda o BS subestima perceptivelmente a incerteza em M, fornece resultados corretos. Estas simulações mostram que o MC é mais confiável que o BS para a propagação de erros em experimentos de mapeamento por eclipse. 11
Figura 6. Gráfico do erro propagado na taxa de acréscimo de massa (Macr ) versus o erro original do modelo. A notaçao é a mesma das Figs. 3 e 5. A linha cheia mostra a relação y = x esperada, as linhas tracejadas e traço-ponto mostram o ajuste da melhor reta ao nível de 1σ para o MC e o BS, respectivamente.
4
Conclusão
A aplicação dos dois métodos aos problemas descritos na seção 3 deixou claro, que diferentemente do que Dhillon & Watson [3] afirmaram, não há evidência de superioridade do método de BS em relação ao MC, pelo menos para os exemplos testados. Nos dois primeiros casos ( seções 3.1 e 3.2 ) mostrou-se que se o método de MC superestima o erro propagado, então o BS o faz da mesma forma. No entanto, o BS apresenta uma variação dos resultados maior que o MC, o que se reflete no tamanho das barras de erros das figuras 3 e 5. Por fim, na seção 3.3 usou-se uma abordagem um pouco diferente que permitiu mostrar que nenhum dos métodos superestima o erro propagado. Em particular, neste caso é possível afirmar que o método de MC segue mais de perto a relação esperada. Nosso trabalho também mostrou que os dois métodos consistentemente produzem erros propagados mais incertos e menos confiáveis quando as incertezas nos dados de entrada são maiores. 12
Os diferentes resultados obtidos para os casos da propagação de erro em experimentos diretos ( ajuste espectral ) e indireto ( imageamento indireto, com vários passos intermediários entre os dados e a grandeza a ser determinada ) indicam que a performance comparativa dos dois métodos pode variar caso a caso, e que simulações específicas dos dois métodos são necessárias para avaliar o desempenho de cada método em cada problema específico. Este trabalho foi parcialmente apoiado pelo CNPq através do processo 300354/96-7 e pelo programa PIBIC/CNPq.
Referências [1] Baptista, R., 2001. Astro tomography, Indirect Imaging Methods in Observational Astronomy, Springer-Verlag, Berlin, p. 307-331 [2] Baptista, R., Steiner, J.E., 1993. A&A, 277, 331 [3] Dhillon, V.S., Watson, C.A., 2001. Astro tomography: Indirect Imaging Methods in Observational Astronomy, Springer-Verlag, Berlin, p. 94 [4] Efron, 1982. The jack knife, the Bootstrap and other Resampling Plans, SIAM, Philadelphia [5] Gunn, J.E., Stryker, L.L., 1983, ApJS, 52, 121 [6] Horne, K., 1985. MNRAS, 213, 129 [7] Horne, K., Cook, M. C., 1985. MNRAS, 214, 307 [8] Kurucz, R.L, 1979. ApJS, 40, 1 [9] Lupton, R., 1993. Statistics in Theory and Practice, Princeton University Press, Princeton [10] Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P., 1992. Numerical Recipes in C, 2nd edition, Cambridge University Press, Cambridge [11] Warner B., 1995. Cataclysmic variable stars. Cambridge University Press, Cambridge
13
