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

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)

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

• I do not understand: if you can simulate $X$ exactly, then you can simulate $1/X$ exactly. The variance you obtain is thus exactly the variance of the rv you do simulate. Sep 17, 2021 at 19:55
• Why "translate" it at all? Generate $X$, then invert it to get $1/X$. Sep 17, 2021 at 20:09
• Your question is difficult to get through. It is unclear what you are trying to do and what the problem is. I suggest you narrow it down and make it smaller. What I currently do not get is why you need the pdf of Y=1/X. You are sampling so you do not need an explicit PDF. You can generate a sample of X's and then compute a sample of Y's from it. Sep 17, 2021 at 20:14
• The answer to your core question: "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 ? " is "yes", and that fact is the basis of the various responses in comments. Sep 17, 2021 at 22:37
• invrcoga=function(n,...)1/rcoga(n,...) Sep 18, 2021 at 9:21

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)
 0.03734464


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

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