Hi,

I’m trying to fit data that represents the difference between latent random variable that are gamma distributed (with the same rate). I’m obviously not succeeding.

I generate my data like this in r. Basically I want to fit `Y`

, assuming I know `alpha2`

, and `beta`

.

```
rgammadiff <- function(N, alpha1=1, alpha2=1, beta=1){
rgamma(N, alpha1, beta= beta) - rgamma(N, alpha2, beta = beta)
}
N <- 10
alpha1 <- 500
alpha2 <- 10
beta <- 1
Y<- rgammadifft(N = N, alpha1, alpha2, beta)
```

I’m following equation (4) from Klar (2014), where the density at z is:

f(z) = c \cdot e^{\beta z} \int_z^\infty x^{\alpha_1-1} (x-z)^{\alpha_2-1} e^{-2\beta x} dx if z >0

where c = \beta^{\alpha_1+\alpha_2}/ (\Gamma(\alpha_1) \Gamma(\alpha2))

I’m ignoring z<0 because in my real data, z is always positive (I guess I need to fit a truncated difference of gammas eventually, but I’m postponing that problem).

And my code in stan is as follows:

```
functions {
real pos(real x, // Function argument
real xc, // Complement of function argument
real[] theta, // parameters
real[] x_r, // data (real)
int[] x_i) { // data (integer)
real alpha1 = theta[1];
real alpha2 = theta[2];
real beta = theta[3];
real z = x_r[1];
return x^(alpha1 - 1) * (x - z)^(alpha2 - 1)*exp(-2 * beta * x);
}
real gammadifft_lpdf(data real y, real alpha1, real alpha2, real beta) {
int x_i[0];
real logc = log(beta)*(alpha1 + alpha2) - lgamma(alpha1) - lgamma(alpha2);
real logint = log(integrate_1d(pos,
y,
positive_infinity(),
{alpha1, alpha2, beta},
{y}, x_i));
return logc + beta * y + logint ;
}
}
data {
int N;
vector[N] Y;
real alpha2;
real beta;
}
parameters {
real<lower=0> alpha1;
}
model {
target += lognormal_lpdf(alpha1 | log(500),1);
for(i in 1:N)
target += gammadifft_lpdf(Y[i] | alpha1, alpha2, beta);
}
```

I’m just trying to make it work, so as I said I assume I know everything except `alpha1`

:

```
lsdata <- list(Y=Y, alpha2 = alpha2, beta = beta, N = length(Y) )
library(cmdstanr)
cmdstan_version()
cmd <- cmdstan_model("gammadifft_lpdf.stan")
fit_cmd <- cmd$sample(
data = lsdata,
seed = 123,
chains = 1,
cores = 4)
```

And I get the following error ` Exception: Error in function boost::math::quadrature::exp_sinh<double>::integrate: The exp_sinh quadrature evaluated your function at a singular point and returned -nan. Please ensure your function evaluates to a finite number over its entire domain. (in '/tmp/RtmpzVP2ah/model-4a7f71da3259.stan', line 13, column 4 to column 57)`

.

(Line 13 is the return of function `pos`

.)

Is it even possible to do this in Stan? The paper of Klar (2014) also mentions an alternative formula by means of Whittaker’s W function, a confluent hypergeometric function (Olver et al., 2010, 13.14.3). Would it make more sense to use this function (somehow from some library from C++)?

(By the way, the R function `dcoga2dim`

(package `coga`

) uses the C++ function `gsl_sf_hyperg_1F1`

to calculate the convolution of two gammas, but I couldn’t find enough documentation to understand it.)

Thanks!!

Bruno