1.) I am interested in computing the variance of this observable $O$ involving the coefficients of spherical harmonics $a_{\ell m}$ and the $C_{\ell}$ which is the variance of an $a_{\ell m}$ :
$$O=\sum_{\ell=\ell_{\min }}^{\ell_{\max }} \sum_{m=-\ell}^{\ell} a_{\ell m}^{2}$$
where $\left(a_{\ell m}\right)$, and $\left(\ell \in\{1, \cdots, N\},|m| \leq \ell\right)$, with $a_{\ell m} \sim \mathcal{N}\left(0, C_{\ell}\right)$ for each $|m| \leq \ell$ .
I recall the properties of a few basic distributions. I use $\sim$ to denote equality in distribution. We have :
- $\mathcal{N}(0, C)^{2} \sim C \chi^{2}(1)=\Gamma\left(\frac{1}{2}, 2 C\right)$,
- $\langle\Gamma(k, \theta)\rangle=k \theta$ and $\operatorname{Var}(\Gamma(k, \theta))=k \theta^{2}$, and
- $\sum_{i=1}^{N} \Gamma\left(k_{i}, \theta\right)=\Gamma\left(\sum_{i=1}^{N} k_{i}, \theta\right)$ for independent summands.
2.) Distribution followed by tUsing points 1 and 3 again, we obtain :
$$
\sum_{m=-\ell}^{\ell}\left(a_{\ell m}\right)^{2}=\sum_{m=-\ell}^{\ell} C_{\ell} \cdot\left(\frac{a_{\ell, m}}{\sqrt{C_{\ell}}}\right)^{2}
$$
$$
\begin{aligned}
&\sim \sum_{m=-\ell}^{\ell} C_{\ell} \operatorname{ChiSq}(1) \\
&\sim C_{\ell} \sum_{m=-\ell}^{\ell} \operatorname{ChiSq}(1) \\
&\sim C_{\ell} \operatorname{ChiSq}(2 \ell+1)
\end{aligned}
$$
$\sim C_{\ell} \operatorname{Gamma}((2 \ell+1) / 2,2)$
$\sim \operatorname{Gamma}\left((2 \ell+1) / 2,2 C_{\ell}\right)$,
We have taken the convention (shape,scale) parameters for Gamma distribution. Given the fact that we consider the random variable :
$$
\sum_{\ell=\ell_{\min }}^{\ell_{\max }} \sum_{m=-\ell}^{\ell} a_{\ell m}^{2}
$$
This sum of random variables $\sum_{\ell=\ell_{\min }}^{\ell_{\max }} \sum_{m=-\ell}^{\ell} a_{\ell m}^{2}$ follows a Moschopoulos distribution : it represents the sum of random variables each one following a Gamma distribution with different shape and scale parameters.
Question 1) Is this reasoning correct ? I mean, I don't know if I can write this step :
$$
\begin{aligned}
&\sim \sum_{m=-\ell}^{\ell} C_{\ell} \operatorname{ChiSq}(1) \\
&\sim C_{\ell} \sum_{m=-\ell}^{\ell} \operatorname{ChiSq}(1) \\
&\sim C_{\ell} \operatorname{ChiSq}(2 \ell+1)
\end{aligned}
$$
Question 2) :
In relationship of Question 1) above, to compute the variance of $O$, I have implemented 2 ways (or called 2 methods) (in R language) with nRed is a bin redshift (like an index of sample) :
# Weighted sum of Chi squared distribution
y3<-array(0,dim=c(nSample_var,nRed));
for (i in 1:nRed) {
for (j in 1:nRow) {
# Try to summing all the random variable
y3[,i] <- y3[,i] + Cl[j,i] * rchisq(nSample_var,df=L[j])
}
}
And another method using coga library (convolution of Gamma distribution, see https://github.com/ChaoranHu/coga ) :
y1 <- array(0, dim=c(nRow))
# y2 is the Cl assimilated to scale parameter
# (shape/rate) convention :
for (i in 1:nRed) {
y3[,i] <- rcoga(nSample, y1, y2)
}
Both methods gives the quasi-same variance, mean and skewness (which is close to 0, so quasi-gaussianity).
What it disturbs me is the first method with weighted squared : how can I justify that sum of Chisquared distribution has a PDF equal to the sum of each Chisquared(1) distribution, weighted by $C_\ell$ ?
Clasically, when we add random variables, the PDF is the convolution of all these random variables, isn't it ?
If you look at the first code snippet, I simply add on $\ell$ the chisquared distribution.
So it seems weird that I get the same results than with coga function (convolution of Gamma distribution).
Where is the error for the second code snippet ? (if there is an error).
Any track/suggestion is welcome.
UPDATE 1:
I think this derivation below is wrong :
$$
\begin{aligned}
\sum_{m=-\ell}^{\ell}\left(a_{\ell m}\right)^{2}&\sim \sum_{m=-\ell}^{\ell} C_{\ell} \operatorname{ChiSq}(1) \\
&\sim C_{\ell} \sum_{m=-\ell}^{\ell} \operatorname{ChiSq}(1) \\
\end{aligned}
$$
I should write instead :
$$
\begin{aligned}
\sum_{m=-\ell}^{\ell}\left(a_{\ell m}\right)^{2}&\sim C_{\ell}\sum_{m=-\ell}^{\ell} \operatorname{ChiSq}(1) \\
&\sim C_{\ell} \operatorname{ChiSq}((2\ell+1) \\
\end{aligned}
$$
Shouldn't I ?
But I have still doubts if I can write with pdf the probability density function of random variable $X$ (here $X=\sum_{m=-\ell}^{\ell}\left(a_{\ell m}\right)^{2}$) :
$$C_\ell X \sim C_\ell\,\text{Distribution}(X)\quad(1)$$
with $\text{Distribution}(X)=\operatorname{ChiSq}((2\ell+1)$.
As you can see, it is not still very clear.
Any clarification is welcome.
UPDATE 2:
I think that @J.Delaney is perfectly right, $C_\ell X \sim C_\ell\,\text{Distribution}(X)$ doesn't make sense.
But once we are in the code and we want to generate sample, can we do like I did with R language :
L <- 2*(array_l)+1
nSample_var <- 100000
y3<-array(0,dim=c(nSample_var))
for (j in 1:nRow) {
# Try to summing all the random variable
y3 <- y3 + Cl[j] * rchisq(nSample_var,df=L[j])
}
?
with Cl[j] $\equiv C_{\ell,j}$
As you can see, I perform an increment sum with $C_\ell$ in front of the function rchisq which is the $\chi^2$ function of distribution : is it a right way to compute samples for random variable $O$ ?
If someone could translate this coding part and a mathematical expression of how are done things from a distribution point of view, this would be fine.
Best Answer
Regarding your question 1) : the result is correct, although your notation is a bit sloppy. first you have:
$$a_{\ell m}^2 /C_{\ell} \sim \chi^2_1$$
By definition, a sum of $k$ independent $\chi^2_1$ random variables has a $\chi^2_k$ distribution (notice the assumption of independence: you have to assume that the $a_{\ell m}$'s are independent) so indeed $$ \sum_{m=-\ell}^\ell a_{\ell m}^2 /C_{\ell} \sim \chi^2_{2\ell +1}$$
Now if $X \sim \chi^2_{2\ell +1}=\Gamma(\ell+\frac{1}{2},2)$ then $C_{\ell}X \sim \Gamma(\ell+\frac{1}{2},2C_{\ell})$ by the scaling property of the Gamma distribution, therefore
$$ \sum_{m=-\ell}^\ell a_{\ell m}^2 \sim \Gamma(\ell+\frac{1}{2},2C_{\ell})$$
as you derived. as for question 2) : you are interested in the variable
$$O = \sum_{\ell=\ell_{min}}^{\ell_{max}} Y_{\ell} $$
where $Y_{\ell} \sim \Gamma(\ell+\frac{1}{2},2C_{\ell})$. Since $O$ is a sum of independent Gamma random variables, it's distribution is indeed a convolution of Gamma distribution. However you can sample that distribution by generating samples from the distributions of the $Y_{\ell}$'s and summing them. This is what you do in your first method
You are generating the $Y_{\ell}$'s using the $\chi^2_{2\ell +1}$ distribution since as we saw before scaling a $\chi^2_{2\ell +1}$ random variable by a constant $C_\ell$ gives you the distribution of $Y_{\ell}$.
In the second method you are using the
cogalibrary that can calculate the convolution of Gamma distributions, however again you are using it to sample the distribution and not actually calculate the pdf (which you can do using thepcogafunction). The simplest method to sample this distribution is probably just by generating the individual Gamma random variables and summing them, so there's a good chance that under the hood thercogafunction is doing exactly what you do in your first method.