Температурный спектр мощности задней формы и здоровье - PullRequest
Новая
       68

Температурный спектр мощности задней формы и здоровье

0 голосов
/ 27 марта 2020

Я пытаюсь повторить эксперимент, приведенный в разделе 4.1 в следующей статье: https://iopscience.iop.org/article/10.1088/0004-637X/697/1/258/meta

Однако, используя те же настройки, я получаю разные задние формы с незначительной частью масса усекается до 0.

Вот мой код. Правильно ли я использую все здоровые части?

В следующем примере я взял NSIDE = 512, L_MAX = 1000 и среднеквадратическое значение шума sqrt (40) мкК. опечатка, потому что я не могу иметь SNR 1 при l = 550 для любого из лучей WMAP - и использовать размер луча Q-диапазона с FWHM 0,51 градуса. Это дает мне SNR 1 вокруг l = 522.

Формула для апостериорного задана уравнением (22) вышеупомянутой статьи. Тем не менее, я предполагаю, что есть несколько опечаток: C_ {l} должен быть заменен на C_ {b}, и перед каждым C_ {b} должен быть коэффициент 2 * pi / (l * (l + 1)) в случае мы биннинг, как дано уравнение (21). Я делаю исправления в своем коде. 

import numpy as np
import healpy as hp
from scipy import stats

L_MAX_SCALARS = 1000
NSIDE = 512
Npix = 12*NSIDE**2
w = 4*np.pi/Npix
noise_var = 40
fwhm_deg = 0.51

#Loading the Cls generated using Class with cosmo params:
#NAMES = ["n_s", "omega_b", "omega_cdm", "100*theta_s", "ln10^{10}A_s", "tau_reio"]
#Values= np.array([0.9665, 0.02242, 0.11933, 1.04101, 3.047, 0.0561])
data = np.load(path_data, allow_pickle=True)
data = data.item()
Cls = data["cls_"]

#Generating beam and map
bl_gauss = hp.gauss_beam(fwhm=(np.pi/180)*fwhm_deg, lmax=L_MAX_SCALARS)
d = hp.synfast(Cls, lmax=L_MAX_SCALARS, nside=NSIDE, fwhm=(np.pi/180)*fwhm_deg, new = True) + np.sqrt(noise_var)*np.random.normal(size=Npix)

# Compute the likelihood of each (possibly) binned power spectrum component: (product of) inverse gamma.
def compute_likelihood(x, alphas, betas, locs, bl_gauss, l_start, l_end):
    return np.prod(stats.invgamma.pdf(x, a=alphas[l_start:l_end], scale=betas[l_start:l_end], loc=locs[l_start:l_end])
                  /bl_gauss[l_start:l_end]**2)


#Trace the posterior of each (possibly) binned power spectrum component
def trace_posterior(d, bins, max_range):
    l_start, l_end = bins
    scale_dl = np.array([l*(l+1)/(2*np.pi) for l in range(L_MAX_SCALARS+1)])
    power_spec = hp.anafast(d, lmax=L_MAX_SCALARS, alm=False)
    power_spec *= scale_dl/bl_gauss**2
    #Compute the parameters of each of the inverse gamma:
    betas = np.array([(2*l+1)*pow_sp/2 for l, pow_sp in enumerate(power_spec)])
    alphas = np.array([(2*l-1)/2 for l in range(L_MAX_SCALARS+1)])
    locs = - w*noise_var*scale_dl/(bl_gauss**2)

    #Getting the normalising constant
    norm, err = scipy.integrate.quad(compute_likelihood, a=0, b=np.inf, args=(alphas, betas ,locs, bl_gauss, l_start, l_end))

    steps = max_range/10000
    xs = np.arange(0, max_range ,steps)
    y = []
    for x in xs:
        result = compute_likelihood(x, alphas, betas, locs, bl_gauss, l_start, l_end)
        y.append(result)

    return np.array(y), xs, norm


bins = [855, 1001]
y, x, norm = trace_posterior(d, bins, 20000)


plt.plot(x, y/norm)
plt.title("Binned posterior ell range " + str(bins[0]) +"-"+str(bins[1]))
plt.show()

Этот код дает мне следующий апостериор для бин-спектра мощности [855-1000]: enter image description here

что противоречит рисунку 4 статьи ... Я понятия не имею, почему я наблюдаю такую ​​разницу.

Правильно ли я использую все здоровые части? Как я создаю луч или генерирую данные? Или вычислить наблюдаемый спектр мощности?

Спасибо,

Габриэль.

...