Finding a right way of sampling 1/X knowing X follows the Moschopoulos distribution (sum of Gamma distribution with different (shape/rate parameters)

Question

I can generate, with COGA R library (with rcoga function), a sample from a random variable following the Moschopoulos distribution (https://arxiv.org/pdf/1806.04059.pdf).

For reminder, Moschopoulos distribution is followed by a sum of Gamma random variables with different (shape/rate) parameters.

We are interested in computing the variance of an observable $$ O=\frac{\sum_{\ell=1}^{N} \sum_{m=-\ell}^{\ell} a_{\ell m}^{2}}{\sum_{\ell=1}^{N} \sum_{m=-\ell}^{\ell}\left(a_{\ell m}^{\prime}\right)^{2}} $$

where $\left(a_{\ell m}, \ell \in\{1, \cdots, N\},|m| \leq \ell\right)$ and $\left(a_{\ell m}^{\prime}, \ell \in\{1, \cdots, N\},|m| \leq \ell\right)$ are independent random variables, with $a_{\ell m} \sim \mathcal{N}\left(0, C_{\ell}\right)$ for each $|m| \leq \ell$ and $a_{\ell m}^{\prime} \sim \mathcal{N}\left(0, C_{\ell}^{\prime}\right)$ for each $|m| \leq \ell .$

Let us consider firstly the numerator of observable $O$: supposing we have independent values $a_{\ell, m} \sim \mathrm{N}\left(0, C_{\ell}\right)$, then the quantity of interest to you has the following distribution: $$ \begin{aligned} Q_{N} & \equiv \sum_{\ell=1}^{N} \sum_{m=-\ell}^{\ell} a_{\ell, m}^{2} \\ &=\sum_{\ell=1}^{N} \sum_{m=-\ell}^{\ell} C_{\ell} \cdot\left(\frac{a_{\ell, m}}{\sqrt{C_{\ell}}}\right)^{2} \\ & \sim \sum_{\ell=1}^{N} \sum_{m=-\ell}^{\ell} C_{\ell} \cdot \mathrm{ChiSq}(1) \\ &=\sum_{\ell=1}^{N} C_{\ell} \sum_{m=-\ell}^{\ell} \mathrm{ChiSq}(1) \\ &=\sum_{\ell=1}^{N} C_{\ell} \cdot \mathrm{ChiSq}(2\ell+1) \end{aligned} $$

Then we can write for the numerator with Gamma distribution and (shape/rate) convention :

$$\sum_{\ell=1}^N C_\ell \cdot \text{ChiSq}(2 \ell + 1)=\sum_{\ell=1}^N \, \text{Gamma}\Big((2 \ell + 1)/2,\dfrac{1}{2\,C_\ell}\Big)$$

This is done easily in R language by the following scritpt (I generate a sammple for each of 5 bins from data) :

library(coga)

my_data <- read.delim("Array_total_WITH_Shot_Noise.txt", header = FALSE, sep = " ")

array_2D <- array(my_data)

nRed <- 5
nRow <- NROW(my_data)

z_ph <- c(0.9595, 1.087, 1.2395, 1.45, 1.688)
b_sp <- c(1.42904922, 1.52601862, 1.63866958, 1.78259615, 1.91956918)
b_ph <- c(sqrt(1+z_ph))
ratio_squared <- (b_sp/b_ph)^2

y1 <- array(0, dim=c(nRow))
y2 <- array(0, dim=c(nRed,nRow))

y2_sp <- array(0, dim=c(nRed,nRow))
y2_ph <- array(0, dim=c(nRed,nRow))

nSample = 100000
y3 <- array(0, dim=c(nSample,nRed))
y3_1 <- array(0, dim=c(nSample,nRed))
y3_2 <- array(0, dim=c(nSample,nRed))

# Shape parameter for Gamma distribution 
y1 <- (2*(array_2D[,1])+1)/2

for (i in 2:6) {
# (shape/rate) convention : 
y2_sp[i-1,] <- 0.5/(array_2D[, i] * (b_sp[i-1]/b_ph[i-1])^2)
y2_ph[i-1,] <- 0.5/(array_2D[, i])
}

# Sampling for numerator
# (shape/rate) convention : 
for (i in 1:nRed) {
  y3_1[,i] <- rcoga(nSample, y1, y2_sp[i,])
  # variance of numerator
  print(paste0('var_numerator = ', var(y3[,i]), ', sigma = ', sd(y3[,i])))
}

The result of standard deviation is about 1e-5, which seems to be acceptable.

Now, I am faced to the main issue of this post, that is to say, to sample the denominator ( $\sum_{\ell=1}^{N} \sum_{m=-\ell}^{\ell}\left(a_{\ell m}^{\prime}\right)^{2}$) :

if I sample a random variable following a distribution and after that, I compute the variance of the inverse of sample (1/X) : will I have the same result if I sample directly the random (1/X) and compute after the variance of the sample ? I say that since I don't know a priori the distribution of 1/X even if I know the distribution of X.

This is not an evident answer. Indeed, by transfer theorem, if we have $X$ follows $f(x)$ and $Y=1/X$ , what is the distribution expression of $Y$ :

We must have : $$f(x) \text{d}x = g(y) \text{d}dy$$

$$\Rightarrow g(y) = f(x) \dfrac{1}{y^2} = f(\dfrac{1}{y}) \dfrac{1}{y^2}$$

How to implement this new PDF $g(y)$ (Probability Density Function) with rcoga (Moschopoulos distribution) ?

Especially, how to express $\dfrac{1}{y}$ into rcoga distribution and how to handle the factorization by $\dfrac{1}{y^2}$ ?

Any help is welcome to do the things correctly and use in a right way the function rcoga by including the $\dfrac{1}{y}$ randow variable within the PDF and add factor $\dfrac{1}{y^2}$ to respect the right distribution expression of $Y=1/X$ random variable.

Up to now, to compute the variance of the observable $O$, I simply did one sampling for numerator and one sampling for denominator and performed the division between both :

# Monte-Carlo for X and Y
# (shape/rate) convention : 
for (i in 1:nRed) {
  y3_1[,i] <- rcoga(nSample, y1, y2_sp[i,])
  y3_2[,i] <- rcoga(nSample, y1, y2_ph[i,])
  y3[,i] <- y3_1[,i] / y3_2[,i] 
}

# print mean, variance and standard deviation
for (i in 1:nRed) {
  print(paste0('mean_fid = ', ratio_squared[i]))
  print(paste0('mean_exp = ', mean(y3[,i])))
  print(paste0('var = ', var(y3[,i]), ', sigma = ', sd(y3[,i])))
}

With this method, the numerator and denominator have a standard deviation about ~ 1e-5 but their ratio has a larger standard deviation (~ 1e-2).

That's why I would like to try to make decrease the standard deviation of observable $O$ in case if I made an error with the approach adopted that I explain from the beginning of this post, that is to say, sampling numerator and denominator without taking into account the distribution of denominator that should be followed by a $1/Y$ random variable.

Don't hesitate to ask me for further precisions if I have been unclear.

PS: the input file Array_total_WITH_Shot_Noise.txt is available on input file

Answer 1

That the distribution of $Y=1/X$ is the same as the distribution of the transform of $X$ by the function\begin{align*} f~: ~&\mathbb R \longmapsto \mathbb R\\&\ x \longmapsto 1/x\end{align*} is a truism as this is how the new random variable $Y$ is defined.

In measure theoretic terms, $X$ is itself a function from a Borelian space to $\mathbb R$:\begin{align*} X~: ~&\Omega \longmapsto \mathbb R\\&\,\omega \longmapsto X(\omega)\tag{1}\end{align*} and the function $Y$ is hence the composition of the functions $X$ and $f$:$$ Y = f\circ X$$

This identity therefore applies to simulation in that simulating $X_1,\ldots,X_n$ and taking $Y_1=1/X_1,\ldots,Y_n=1/X_n$ produces a genuine sample of $Y$'s from the correct distribution. One can for instance think of $\omega$ in (1) as a realisation of a Uniform $(0,1)$ variate used to simulate $X=X(\omega)$ and thus $Y=1/X(\omega)$.

As an example, take a Pareto distribution $$X\sim \frac{3}{x^4}\mathbb I_{(1,\infty)}(x)$$ which has a finite variance. Then $Y=1/X$ enjoys a Beta distribution $$Y\sim 3y^2\mathbb I_{(0,1)}(x)\tag{2}$$ endowed with a finite variance as well. If generating $Y$ as the inverse of a Pareto,

y=1/(1-runif(nSample)^{-1/3})

one produces a sample distributed according to (2).

That the variance of the ratio is larger than the variances of both numerator and denominator is a feature of that distribution, not of a mistake in the implementation. It does happen that a random variable $X$ has a finite variance, while the inverse random variable $1/X$ has an infinite or even no variance. Consider the example of a standard Normal variate, for instance.

In the Pareto example,

> var(y)
[1] 0.03734464

corresponds to the Beta $\mathcal B(3,1)$ variance:

> 3*1/(3+1)^2/(3+1+1)
[1] 0.0375