- 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}$ ?
- 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
- 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
andy3_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.