[R-br] Integração numérica de valores tabulados

[Alexandre Castagna]

Apenas uma correção... a velocidade da função de integração na minha
máquina é 0.02 segundos e combinado com a de otimização é 0.2 segundos.


Em 13 de agosto de 2014 22:27, Alexandre Castagna <
alexandre.castagna em gmail.com> escreveu:

> Prezados,
>
> Preciso aplicar repetidas vezes uma integral em um distribuição de
> probabilidade e por isso preciso de um método muito rápido (são milhares de
> repetições). Porém não existe uma função matemática que descreva a
> distribuição, os valores estão tabulados.
>
> A função de integração que construí é muito lenta (0.2 segundos na minha
> máquina); queria reduzir essa velocidade em 10 ou 100 vezes. Apresento a
> função com os valores tabulados:
>
>     angles <- log(c(0.001, 0.01, 0.1, 0.12589, 0.15849, 0.19953, 0.25119,
> 0.31623, 0.39811,
>     0.50119, 0.63096, 0.79433, 1, 1.2589, 1.5849, 1.9953, 2.5119, 3.1623,
>     3.9811, 5.0119, 6.3096, 7.9433, 10, 15, 20, 25, 30, 35, 40, 45, 50,
> 55,
>     60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130,
> 135,
>     140, 145, 150, 155, 160, 165, 170, 175, 180)*pi/180)
>
>     mc.petold <- log(c(869105, 39180, 1.76661e3, 1.29564e3, 9.50172e2,
> 6.99092e2, 5.13687e2,
>     3.76373e2, 2.76318e2, 2.18839e2, 1.44369e2, 1.02241e2, 7.16082e1,
>     4.95803e1, 3.39511e1, 2.28129e1, 1.51622e1, 1.00154e1, 6.57957,
> 4.29530,
>     2.80690, 1.81927, 1.15257, 4.89344e-1, 2.44424e-1, 1.47151e-1,
> 8.60848e-2,
>     5.93075e-2, 4.20985e-2, 3.06722e-2, 2.27533e-2, 1.69904e-2,
> 1.31254e-2,
>     1.04625e-2, 8.48826e-3, 6.97601e-3, 5.84232e-3, 4.95306e-3,
> 4.29232e-3,
>     3.78161e-3, 3.40405e-3, 3.11591e-3, 2.91222e-3, 2.79696e-3,
> 2.68568e-3,
>     2.57142e-3, 2.47603e-3, 2.37667e-3, 2.32898e-3, 2.31308e-3,
> 2.36475e-3,
>     2.50584e-3, 2.66183e-3, 2.83472e-3, 3.03046e-3, 3.09206e-3,
> 3.15366e-3))
>
> cum.beta.ph.p.petzold <- function(psi)
>    {
>     Integral <- integrate(function(x) { exp(approx(angles,
>                  mc.petold, xout = log(x))$y)*sin(x) },
>                  lower = 1.745329e-05, upper = psi)$value*2*pi
>     return(Integral)
>    }
>
> cum.beta.ph.p.petzold(psi = pi)
>
> Onde psi deve ser entrado em radianos. A integral de "0" a pi é apenas
> aproximadamente 1, por conta de uma simplificação dos valores e pelo limite
> inferior de integração em 0.001º (1.745329e-05 rad), já que o valor é
> infinito em 0. A normalização para 1 é fácil de ser realizada.
>
> A função é aplicada milhares de vezes em um programa de simulação com
> números aleatórios, onde para tirar uma amostra aleatória dessa
> distribuição, eu uso:
>
> randon.n <- runif(1)
> optimize( function(x) abs(cum.beta.ph.p.petzold(x)-randon.n), interval =
> c(0, pi), tol = 1e-3)$minimum
>
> Agradeço se alguém puder oferecer uma sugestão.
>
> Obrigado,
> --
> Alexandre Castagna Mourão e Lima
> Laboratório de Radioecologia e Mudanças Globais - UERJ
> Rua São Francisco Xavier, 524, Maracanã
> Pav. Haroldo Lisboa da Cunha, Sub-solo
>



-- 
Alexandre Castagna Mourão e Lima
Laboratório de Radioecologia e Mudanças Globais - UERJ
Rua São Francisco Xavier, 524, Maracanã
Pav. Haroldo Lisboa da Cunha, Sub-solo
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20140813/b1c04d83/attachment.html>