Как я могу улучшить интеграцию и параметризацию свернутых распределений?

Я пытаюсь решить для параметров гамма-распределения, которое свернуто как с нормальным, так и с логнормальным распределением. Я могу экспериментально получить параметры как для нормальной, так и для логнормальной составляющих, поэтому я просто хочу найти параметры гаммы.

Я попытался 3 подхода к этой проблеме:

1) создание свернутых случайных наборов данных (т. е. rnorm()+rlnorm()+rgamma()) и использование регрессии наименьших квадратов на линейных или логарифмических гистограммах данных (не показано, но было очень смещено ГСЧ и вообще не оптимизировалось).

2) численное интегрирование сверточных функций «грубой силой» (пример кода № 1)

3) подходы к численному интегрированию с пакетом distr. (пример кода #2)

У меня был ограниченный успех со всеми тремя подходами. Важно отметить, что эти подходы, по-видимому, хорошо работают для «номинальных» значений гамма-параметров, но все они начинают давать сбои, когда k (форма) низка, а тета (масштаб) высока — именно здесь находятся мои экспериментальные данные. пожалуйста, найдите примеры ниже.

Прямая числовая интеграция

# make the functions
f.N <- function(n) dnorm(n, N[1], N[2])
f.L <- function(l) dlnorm(l, L[1], L[2])
f.G <- function(g) dgamma(g, G[1], scale=G[2])

# make convolved functions
f.Z <- function(z) integrate(function(x,z) f.L(z-x)*f.N(x), -Inf, Inf, z)$value # L+N
f.Z <- Vectorize(f.Z)
f.Z1 <- function(z) integrate(function(x,z) f.G(z-x)*f.Z(x), -Inf, Inf, z)$value # G+(L+N)
f.Z1 <- Vectorize(f.Z1)

# params of Norm, Lnorm, and Gamma
N <- c(0,5)
L <- c(2.5,.5)
G <- c(2,7) # this distribution is the one we ultimately want to solve for.
# G <- c(.5,10) # 0<k<1
# G <- c(.25,5e4) # ballpark params of experimental data 

# generate some data
set.seed(1)
rN <- rnorm(1e4, N[1], N[2])
rL <- rlnorm(1e4, L[1], L[2])
rG <- rgamma(1e4, G[1], scale=G[2])
Z <- rN + rL
Z1 <- rN + rL + rG

# check the fit
hist(Z,freq=F,breaks=100, xlim=c(-10,50), col=rgb(0,0,1,.25))
hist(Z1,freq=F,breaks=100, xlim=c(-10,50), col=rgb(1,0,0,.25), add=T)
z <- seq(-10,50,1)
lines(z,f.Z(z),lty=2,col="blue", lwd=2) # looks great... convolution performs as expected.
lines(z,f.Z1(z),lty=2,col="red", lwd=2) # this works perfectly so long as k(shape)>=1

# I'm guessing the failure to compute when shape 0 < k < 1 is due to 
# numerical integration problems, but I don't know how to fix it.
integrate(dgamma, -Inf, Inf, shape=1, scale=1) # ==1
integrate(dgamma, 0, Inf, shape=1, scale=1) # ==1
integrate(dgamma, -Inf, Inf, shape=.5, scale=1) # !=1
integrate(dgamma, 0, Inf, shape=.5, scale=1) # != 1

# Let's try to estimate gamma anyway, supposing k>=1
optimFUN <- function(par, N, L) {
  print(par)
  -sum(log(f.Z1(Z1[1:4e2])))
}
f.G <- function(g) dgamma(g, par[1], scale=par[2])
fitresult <- optim(c(1.6,5), optimFUN, N=N, L=L)
par <- fitresult$par
lines(z,f.Z1(z),lty=2,col="green3", lwd=2) # not so great... likely better w/ more data, 
# but it is SUPER slow and I observe large step sizes.

Попытка свертки через дистрибутив

# params of Norm, Lnorm, and Gamma
N <- c(0,5)
L <- c(2.5,.5)
G <- c(2,7) # this distribution is the one we ultimately want to solve for.
# G <- c(.5,10) # 0<k<1
# G <- c(.25,5e4) # ballpark params of experimental data 

# make the distributions and "convolvings'
dN <- Norm(N[1], N[2])
dL <- Lnorm(L[1], L[2])
dG <- Gammad(G[1], G[2])
d.NL <- d(convpow(dN+dL,1))
d.NLG <- d(convpow(dN+dL+dG,1)) # for large values of theta, no matter how I change 
# getdistrOption("DefaultNrFFTGridPointsExponent"), grid size is always wrong.

# Generate some data
set.seed(1)
rN <- r(dN)(1e4)
rL <- r(dL)(1e4)
rG <- r(dG)(1e4)
r.NL <- rN + rL
r.NLG <- rN + rL + rG

# check the fit
hist(r.NL, freq=F, breaks=100, xlim=c(-10,50), col=rgb(0,0,1,.25))
hist(r.NLG, freq=F, breaks=100, xlim=c(-10,50), col=rgb(1,0,0,.25), add=T)
z <- seq(-10,50,1)
lines(z,d.NL(z), lty=2, col="blue", lwd=2) # looks great... convolution performs as expected.
lines(z,d.NLG(z), lty=2, col="red", lwd=2) # this appears to work perfectly 
# for most values of K and low values of theta

# this is looking a lot more promising... how about estimating gamma params?
optimFUN <- function(par, dN, dL) {
  tG <- Gammad(par[1],par[2])
  d.NLG <- d(convpow(dN+dL+tG,1))
  p <- d.NLG(r.NLG)
  p[p==0] <- 1e-15 # because sometimes very low probabilities evaluate to 0... 
# ...and logs don't like that.
  -sum(log(p))
}
fitresult <- optim(c(1,1e4), optimFUN, dN=dN, dL=dL)
fdG <- Gammad(fitresult$par[1], fitresult$par[2])
fd.NLG <- d(convpow(dN+dL+fdG,1))
lines(z,fd.NLG(z), lty=2, col="green3", lwd=2) ## this works perfectly when ~k>1 & ~theta<100... but throws
## "Error in validityMethod(object) : shape has to be positive" when k decreases and/or theta increases 
## (boundary subject to RNG).

Можно ли ускорить интеграцию в примере 1? можно ли увеличить размер сетки в примере 2 (пакет дистрибутива)? как я могу решить проблему k‹1? могу ли я изменить масштаб данных таким образом, чтобы облегчить оценку при высоких значениях тета? Есть ли лучший способ все вместе? Помощь!