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

[Alexandre Castagna]

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
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20140813/55ac1417/attachment.html>