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
Best Answer
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,
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,
corresponds to the Beta $\mathcal B(3,1)$ variance: