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

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 = " ")

# 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 :

 "mean_fid = 1.04219529123889"
 "mean_exp = 1.04221682330884"
 "numerator : var = 0.0121217744959571, sigma = 0.110098930494157"
 "denominator : var = 0.0111600755430205, sigma = 0.10564125871562"
 "var = 4.00480985775818e-05, sigma = 0.00632835670435713"
 "mean_fid = 1.11582790061653"
 "mean_exp = 1.1158459102749"
 "numerator : var = 0.0100610443302181, sigma = 0.100304757266134"
 "denominator : var = 0.00808064436739286, sigma = 0.0898924043921001"
 "var = 4.58718033043056e-05, sigma = 0.0067728726035786"
 "mean_fid = 1.19903460255297"
 "mean_exp = 1.19905585384901"
 "numerator : var = 0.0080798197066825, sigma = 0.0898878173429664"
 "denominator : var = 0.00561998564225798, sigma = 0.0749665634950541"
 "var = 5.30041218254675e-05, sigma = 0.00728039297191213"
 "mean_fid = 1.29699960571217"
 "mean_exp = 1.29702702713326"
 "numerator : var = 0.00569332745066, sigma = 0.0754541413751426"
 "denominator : var = 0.00338435913663655, sigma = 0.0581752450500774"
 "var = 6.20670352876997e-05, sigma = 0.00787826346904568"
 "mean_fid = 1.37081318333552"
 "mean_exp = 1.37083587027734"
 "numerator : var = 0.00476304996098467, sigma = 0.069014853191068"
 "denominator : var = 0.00253473732103918, sigma = 0.050346174840192"
 "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}$$ ?

1. 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

1. 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.

• Hi: It seems that you have a few questions. I can't answer the first because I don't see where the 2l + 1 is coming from. For the second, you can't write that if Y = X_1 + X_2 then pdf of Y = pdf of X_1 + pdf of X_2 because it's not necessarily true. Take the uniform(0,1) case. The sum of 2 uniforms is not uniform. So, equality doesn't hold when one is talking about pdf's. The idea is that sometimes the sum can follow another distribution. Suppose you have n rvs and each of them is N(0,1). Then, if you square each of them and sum them, that sum will have a chi squared n distribution. Sep 3, 2021 at 15:14
• @mlofton . Thanks for yout quick answer. The $(2\ell+1)$ comes from the summing $\sum_{m=-\ell}^{\ell}$. I think I can't include the factor $(2\ell+1)$ in $\chi^2$ since there is a dependence of $\ell$ with $C_{\ell}$ but I post here to check this reasoning. Best regards
– user226073
Sep 3, 2021 at 15:44
• I'm still not totally clear on the first part of your question but, in general, don't use equalities with distributions. Think of it as "distributed as" rather than an equality. I don't know what your prob-stat background is but, if you take a math stat sequence ( 2 semesters ) that's calculus based, you'll learn a lot about these distributional concepts and the notion of sums of rv's. For example, the sum of 2 uniforms has a distribution that is a triangle. As you sum more and more of them, it becomes more and more hump shaped like the normal distribution. I hope that helps a little. Sep 4, 2021 at 2:15
• Just one more thing. If you're one of those people who can do self-study ( and don't need the positives of a classroom setting and homework etc ) and have a decent calculus background, then the text by casella and berger ( I forget the name. it's on amazon ) is probably a good recommendation for the kinds of concepts being discussed here. Good luck. Sep 4, 2021 at 2:19
• @mlofton. Thans for your advices. As you can see, I used "$\stackrel{d}{=}$" symbol to mean that it follows a given distribution. Do you agree with this notation ? Best regards
– user226073
Sep 4, 2021 at 2:22

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

• Thanks for your detailled answer. Why couldn't we write after your last equation : $= \sum_{\ell=1}^N C_\ell \cdot \text{ChiSq}(2 \ell + 1)=\sum_{\ell=1}^N C_\ell \cdot \text{Gamma}((2 \ell + 1)/2,1/2)$ with shape/rate convention ?
– user226073
Sep 5, 2021 at 9:52
• You can certainly write it that way if you want; it doesn't simplify things at all, but that is a valid form too.
– Ben
Sep 5, 2021 at 9:56
• If I can write $\sum_{\ell=1}^N C_\ell \cdot \text{ChiSq}(2 \ell + 1)=\sum_{\ell=1}^N C_\ell \cdot \text{Gamma}((2 \ell + 1)/2,1/2)$, I can include the scale factor $C_\ell$ into Gamma. Then, I should be able to write : $\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)$ with (shape/rate) notation for Gamma distrbution, shouldn't I ?
– user226073
Sep 5, 2021 at 11:38
• Yes, you can do that, but you can't bring the summation of scales into the chi-squared.
– Ben
Sep 5, 2021 at 20:40
• Moschopoulos (1985) studied the general form for the distribution of a sum of independent gamma random variables, allowing different shape parameters. So naturally, that corresponds to this case. I've not heard anyone call this the Moschopoulos distribution, but sure, those are the same case.
– Ben
Sep 13, 2021 at 22:32