2016-03-24 4 views
-1
fun1 = function(y,mu=mu0,lsig=lsig0) { 
    res = 1/(exp(-y)+1)^2 * 1/sqrt(2*pi)/exp(lsig) * exp(-(y-mu)^2/2/exp(lsig)^2) 
    return(res) 
} 

fun4 = function(para=c(mu1,lsig1)) { 
    mu1 = para[1] 
    lsig1 = para[2] 
    res = n1 * log(noze(integrate(fun1,-Inf,Inf,mu=mu1,lsig=lsig1)$value)) + 
    n3 * log(noze(integrate(fun2,-Inf,Inf,mu=mu1,lsig=lsig1)$value)) + 
    n2 * log(noze(integrate(fun3,-Inf,Inf,mu=mu1,lsig=lsig1)$value)) 
    return(-res) 
} 

noze = function(x) { 
    if (x < 1e-100) { x = 1e-100 } 
    return(x) 
} 

optim(c(0.5,2),fun4,method="L-BFGS-B",lower=c(-5,-3),upper=c(3.5,3.5))$par 

Мне нужно найти два параметра функции funx, в котором используется интеграл 'fun1.' («fun2» и «fun3» немного отличаются от «fun1»)Интеграл, вероятно, расходящийся

Я столкнулся с ошибкой «Ошибка в интеграции (fun1, -Inf, Inf, -3.9538, -3): Интеграл, вероятно,

Используя диаграмму рассеяния, я обнаружил, что fun1 близок к нулю почти везде, кроме (-4.2, -3.7). Таким образом, интегрирование для этого интервала дает (приблизительно) правильный интеграл.

> integrate(fun1,-4.2,-3.6,-3.9538,-3) 
0.0003558953 with absolute error < 3e-11 

Это может быть подтверждено с помощью близлежащих значений параметров

> integrate(fun1,-Inf,Inf,-3.9538,-3.1) 
0.0003555906 with absolute error < 2.6e-05 
> integrate(fun1,-Inf,Inf,-3.9538,-2.85) 
0.0003564842 with absolute error < 3.7e-06 

Если интервал слишком велик, это дает неправильным интеграл.

> integrate(fun1,-5,5,-3.9538,-3) 
0.0003558953 with absolute error < 2.3e-08 
> integrate(fun1,-15,15,-3.9538,-3) 
3.492547e-11 with absolute error < 6.5e-11 
> integrate(fun1,-30,30,-3.9538,-3) 
1.980146e-41 with absolute error < 3.4e-41 
> integrate(fun1,-50,50,-3.9538,-3) 
0 with absolute error < 0 
> integrate(fun1,-Inf,Inf,-3.9538,-3) 
Error in integrate(fun1, -Inf, Inf, -3.9538, -3) : 
    the integral is probably divergent 

Если я должен интегрировать только один раз, я могу найти интервал, где «fun1» является достаточно большим и интегрировать только для этого интервала. Но проблема в том, что я использую функцию оптимизации, которая пытается найти различные параметры для поиска минимизатора «fun4».

Использование (-Inf, Inf) дает ошибку и достаточно широкий интервал дает неверные интегралы.

Есть ли хороший способ решить эту проблему?

+0

Вы не должны объединить подразделения. Уравнение не однозначно. Например, '1/2/3' может быть либо' (1/2)/3', либо '1/(2/3)', с разными результатами. 'R' анализирует код слева направо, поэтому вы получаете' (1/2)/3', но неясно, действительно ли это то, что вы намереваетесь. Вместо этого используйте скобки или конвертируйте в умножения (например, '1/(2/3) = (1/2) * 3'). – RHertel

+0

Ваша функция f4 относится к n1, n2 и n3, но они не определены в вашем коде. – G5W

ответ

-1

сверток с гауссовым ядром может быть решен с помощью Gauss-Hermite integration, и есть R пакет для этого: https://cran.r-project.org/web/packages/gaussquad/gaussquad.pdf

Некоторого тест-код:

library(gaussquad) 

n.quad <- 128 # integration order 

# get the particular weights/abscissas as data frame with 2 observables and n.quad observations 
rule <- ghermite.h.quadrature.rules(n.quad, mu = 0.0)[[n.quad]] 

# test function - integrate 1 over exp(-x^2) from -Inf to Inf 
# should get sqrt(pi) as an answer 
f <- function(x) { 
    1.0 
} 

q <- ghermite.h.quadrature(f, rule) 
print(q - sqrt(pi))