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