2 approaches for Monte-Carlo : weighted sum of $\chi^2$ distribution and Moschopoulos distribution with Gamma distribution

Question

1. If I take as definition of $a_{lm}$ following a normal distribution with mean equal to zero and $C_\ell=\langle a_{lm}^2 \rangle=\text{Var}(a_{lm})$, and if I have a sum of $\chi^2$, can I write the 2 lines below (We use $\stackrel{d}{=}$ to denote equality in distribution).

Important remark : $C_{\ell}$ depends on $\ell$ : $C_{\ell} = C_{\ell}(\ell)$

\begin{aligned} \sum_{\ell=1}^{N} \sum_{m=-\ell}^{\ell} a_{\ell m}^{2} & \stackrel{d}{=} \sum_{\ell=1}^{N}\,C_\ell \chi^{2}\left(2\ell+1)\right) \end{aligned}

For the moment, on a coding technical point of view, I did :

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

# Number of redshift
nRed <- 5
# Number of multipoles for each of nRed redshift
nRow <- NROW((my_data))

nSample_var <- 100*nRow
nSample_mc <- 1000

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

array_2D <- array(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

# Multipoles
y1 <- array(0, dim=c(nRow))
y1 <- (2*(array_2D[,1])+1)

# Compute degree of freedom
for (i in 2:6) {
y2_sp[,i-1] <- array_2D[, i] * ratio_squared[i-1]
y2_ph[,i-1] <- array_2D[, i]
}

color_curve <- c("red", "green", "grey", "black")

# Declare arrays
y3_1<-array(0, dim=c(nSample_var*nSample_mc,nRed));y3_2<-array(0, dim=c(nSample_var*nSample_mc,nRed));
y3 <- array(0, dim=c(nSample_var*nSample_mc,nRed))

# For random
set.seed(123)
for (i in 1:nRed) {
    y3_1[,i] <- y2_sp[,i] * replicate(nSample_mc, rchisq(nSample_var,df=sum(y1)))
    y3_2[,i] <- y2_ph[,i] * replicate(nSample_mc, rchisq(nSample_var,df=sum(y1)))
    y3[,i] <- y3_1[,i]/y3_2[,i]
  print(paste0('mean_fid = ', ratio_squared[i]))
  print(paste0('mean_exp = ', mean(y3[,i])))
  print(paste0('numerator : var = ', var(y3_1[,i]), ', sigma = ', sd(y3_1[,i])))
  print(paste0('denominator : var = ', var(y3_2[,i]), ', sigma = ', sd(y3_2[,i])))
  print(paste0('var = ', var(y3[,i]), ', sigma = ', sd(y3[,i])))
  }

and I get mean computed in good agreement with theoritical values :

[1] "mean_fid = 1.04219529123889"
[1] "mean_exp = 1.04221682330884"
[1] "numerator : var = 0.0121217744959571, sigma = 0.110098930494157"
[1] "denominator : var = 0.0111600755430205, sigma = 0.10564125871562"
[1] "var = 4.00480985775818e-05, sigma = 0.00632835670435713"
[1] "mean_fid = 1.11582790061653"
[1] "mean_exp = 1.1158459102749"
[1] "numerator : var = 0.0100610443302181, sigma = 0.100304757266134"
[1] "denominator : var = 0.00808064436739286, sigma = 0.0898924043921001"
[1] "var = 4.58718033043056e-05, sigma = 0.0067728726035786"
[1] "mean_fid = 1.19903460255297"
[1] "mean_exp = 1.19905585384901"
[1] "numerator : var = 0.0080798197066825, sigma = 0.0898878173429664"
[1] "denominator : var = 0.00561998564225798, sigma = 0.0749665634950541"
[1] "var = 5.30041218254675e-05, sigma = 0.00728039297191213"
[1] "mean_fid = 1.29699960571217"
[1] "mean_exp = 1.29702702713326"
[1] "numerator : var = 0.00569332745066, sigma = 0.0754541413751426"
[1] "denominator : var = 0.00338435913663655, sigma = 0.0581752450500774"
[1] "var = 6.20670352876997e-05, sigma = 0.00787826346904568"
[1] "mean_fid = 1.37081318333552"
[1] "mean_exp = 1.37083587027734"
[1] "numerator : var = 0.00476304996098467, sigma = 0.069014853191068"
[1] "denominator : var = 0.00253473732103918, sigma = 0.050346174840192"
[1] "var = 6.91620589647813e-05, sigma = 0.00831637294526774"

Do you think this approach is correct compared to the distribution following by random variable $Z=\sum_{\ell=1}^{N} \sum_{m=-\ell}^{\ell} a_{\ell m}^{2}$ ?

2. I don't understand why I can't reproduce the same results by reasoning on the Moschopoulos distribution which describes the law followed by a sum of Gamma distribution.

Indeed, The random variable $Y_{l}=\sum_{m=-\ell}^{m=\ell}\,a_{lm}^2$ follows a Gamma distribution $\text{Gamma}((2\ell+1)/2,2C_{\ell})$.

AND

The sum of $N$ random variables $Y_{l}$ following each one a $\text{Gamma}((2\ell+1)/2,2C_{\ell})$ distribution follows a Moschopoulos distribution ((from paper The computation of Moschopoulos on https://www.ism.ac.jp/editsec/aism/pdf/037_3_0541.pdf)

?

Do you agree with this formulation ? if not, there are subtilities that I have not yet grasped.

Technically, If I take for Gamma distribution the notation (shape/rate), is it enough to inverse the scale parameter to get an equivalent expression with the other notation (shape/scale) ? so, I could write in this case :

The sum of $N$ random variables $Y_{l}$ following each one a $\text{Gamma}((2\ell+1)/2,1/(2C_{\ell}))$ distribution follows a Moschopoulos distribution

4. The goal of all this computation is to estimate numerically the variance of the ratio $X/Y$ with $X$ and $Y$ following a Moschopoulos distribution (respectively represented by y3_1 and y3_2 vectors)

I tried to perfom this with coga library in R language :

# For random
set.seed(123)
for (i in 1:5) {
  y3_1 <- c(rcoga(nSample, y1, y2_sp[i,]))
  y3_2 <- c(rcoga(nSample, y1, y2_ph[i,]))
  # Ratio of 2 samples
  y3 = y3_1/y3_2
  y4[,i] <- y3
  print(paste0('var1 = ', var(y3_1), ', sigma = ', sd(y3_1)))
  print(paste0('var2 = ', var(y3_2), ', sigma = ', sd(y3_2)))
  print(paste0('var = ', var(y4[,i]), ', sigma = ', sd(y4[,i])))
} 

I don't understand why I have a such difference with the case 1) considered above (weighted sum of $\chi^2$ distribution) :

Here what I get with sum of Gamma distribution (Moschopoulos distribution) :

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

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

nSample_var <- 1000000
nSample_mc <- 1000000

y1 <- (2*(array_2D[,1])+1)/2

for (i in 2:6) {
y2_sp[,i-1] <- 1/(2 * array_2D[, i] * ratio_squared[i-1])
y2_ph[,i-1] <- 1/(2 * array_2D[, i])
}

# Declare arrays
y3_1<-array(0, dim=c(nSample_var,nRed));y3_2<-array(0, dim=c(nSample_var,nRed));
y3<-array(0, dim=c(nSample_var,nRed));

for (i in 1:nRed) {
  y3_1 <- c(rcoga(nSample_var, y1, y2_sp[,i]))
  y3_2 <- c(rcoga(nSample_var, y1, y2_ph[,i]))
  # Ratio of samples
  y3[,i] <- y3_1/y3_2
  # Variance of ratio
  print(paste0('mean_fid = ', ratio_squared[i]))
  print(paste0('mean_exp = ', mean(y3[,i])))
  print(paste0('numerator : var = ', var(y3_1), ', sigma = ', sd(y3_1)))
  print(paste0('denominator : var = ', var(y3_2), ', sigma = ', sd(y3_2)))
  print(paste0('var = ', var(y3[,i]), ', sigma = ', sd(y3[,i])))
  print(paste0('mean_y3 = ', mean(y3[,i])))
}

Any help to try to reproduce the reesults of case 1) with Moschopoulos distribution is welcome.

Answer 1

The main problem here is that your notation is a mess, and this is getting in the way of determining the result. In particular, when writing the chi-squared distribution, you need to separate the degrees-of-freedom parameter from any scaling parameter rather than mixing these together. Here is a better way to set things out. Supposing we have independent values $a_{\ell,m} \sim \text{N}(0, C_\ell)$, then the quantity of interest to you has the following distribution:

$$\begin{align} Q_N &\equiv \sum_{\ell=1}^N \sum_{m=-\ell}^\ell a_{\ell,m}^2 \\[6pt] &= \sum_{\ell=1}^N \sum_{m=-\ell}^\ell C_\ell \cdot \bigg( \frac{a_{\ell,m}}{\sqrt{C_\ell}} \bigg)^2 \\[6pt] &\sim \sum_{\ell=1}^N \sum_{m=-\ell}^\ell C_\ell \cdot \text{ChiSq}(1) \\[6pt] &= \sum_{\ell=1}^N C_\ell \sum_{m=-\ell}^\ell \text{ChiSq}(1) \\[6pt] &= \sum_{\ell=1}^N C_\ell \cdot \text{ChiSq}(2 \ell + 1). \\[6pt] \end{align}$$

The result is a weighted sum of chi-squared random variables with respective degrees-of-freedom value $3, 5, 7, ..., 2N+1$. You cannot simplify this further without specifying a form for the scaling values $C_1,...,C_N$; a weighted sum of chi-squared random variables with different degrees-of-freedom does not reduce to a gamma random variable. So no, you cannot bring the weighted sum inside the degrees-of-freedom parameter for the chi-squared distribution. The general distribution for a weighted sum of chi-squared random variables is complicated, and its density function does not have a closed form representation. Bodenham and Adams (2015) examine approximations to this distribution.

(There is another way to tell that your asserted distribution is wrong. Observe that the quantity of interest involves a summation over the index $\ell$ --- because this is an index variable, the quantity of interest is not a function of $\ell$. Since this value does not enter into the function anywhere, your asserted distribution, which is a function of $\ell$, must be wrong. It is good practice to learn to recognise when a variable in an equation is merely an index or integrand value that gets summed or integrated out of the function.)