Statistical error propagation

[ ]:
from radiocalibrationtoolkit import *
from IPython.display import Markdown as md
[ ]:
# some global plot settings
plt.rcParams["axes.labelweight"] = "bold"
plt.rcParams["font.weight"] = "bold"
plt.rcParams["font.size"] = 16
plt.rcParams["legend.fontsize"] = 14

plt.rcParams["xtick.major.width"] = 2
plt.rcParams["ytick.major.width"] = 2

plt.rcParams["xtick.major.size"] = 5
plt.rcParams["ytick.major.size"] = 5

plt.rcParams["xtick.labelsize"] = 14
plt.rcParams["ytick.labelsize"] = 14
[ ]:
# create class instances
galactic_map_inst = GlobalSkyModel2016(freq_unit="MHz")
antenna_inst = AntennaPattern("./antenna_setup_files/SALLA_EW.xml")

# some constants
update_antenna_conventions = {
    "shift_phi": -90,
    "flip_theta": True,
    "flip_phi": False,
    "in_degrees": True,
    "add_invisible_sky": True,
}
LATITUDE = -35.206667
NSIDE = 64

[ ]:
# first part of this examples will use this frequency and LST
frequency_MHz = 45
lst = 22
[ ]:
# create maps
rotation_parameters = create_rotation_parameters(lst, LATITUDE)
rotator = Rotator(coord=["G", "C"], rot=rotation_parameters)

antenna_map_inst = antenna_inst.convert2hp(
    frequency=frequency_MHz, quantity="absolute", **update_antenna_conventions
)
antenna_map = antenna_map_inst.get_map(rotator=rotator)

galactic_map = hp.ma(
    hp.pixelfunc.ud_grade(galactic_map_inst.generate(frequency_MHz), NSIDE)
).copy()
galactic_map.mask = create_local_mask(NSIDE, rotation_parameters)

Propagation to effective temperature and voltage squared spectral denisty

[ ]:
# set up the relative standard deviations on galactic emission and the antenna gain
rstd_galactic_map = 0.5
rstd_antenna_map = 0.

Error propagation from sky temperature and antenna gain to effective temperature

[ ]:
# analytical calculation
rstd_teff = np.sqrt(rstd_galactic_map**2 + 4*rstd_antenna_map**2)

\begin{equation} \rho_{T\mathrm{eff}}(t,f,\alpha,\delta) = \frac{\sigma_{T\mathrm{eff}}(t,f,\alpha,\delta)}{T_{\mathrm{eff}}(t,f,\alpha,\delta)} = \sqrt{ \rho_{T}^2(t,\alpha,\delta) + 4 \rho_H^2 (t,f,\alpha,\delta) } \end{equation}

[ ]:
# experimental test by distorting the galactic and antenna maps
r = (
    distort_array(galactic_map, rstd_galactic_map)
    * distort_array(antenna_map, rstd_antenna_map) ** 2
) / (galactic_map * antenna_map**2)
[ ]:
# compare
print('Experimentally propagated: {:.2f}'.format(np.std(r)))
print('Analytical calculation: ', rstd_teff)
[ ]:
# just a plot comparison of how the distorted effective temperature by 50% looks like
teff = galactic_map * antenna_map**2
teff_distored = distort_array(teff, rstd=rstd_teff)

fontsize = {
    "xlabel": 22,
    "ylabel": 22,
    "title": 14,
    "xtick_label": 22,
    "ytick_label": 22,
    "cbar_label": 22,
    "cbar_tick_label": 22,
}

for m in [teff, teff_distored]:
    projview(
        m,
        cmap="jet",
        min=3500,
        max=100000,
        norm="log",
        graticule=True,
        graticule_labels=True,
        unit="[$m^2$K]",
        xlabel="RA",
        ylabel="DEC",
        cb_orientation="vertical",
        latitude_grid_spacing=30,
        xtick_label_color="white",
        title='',
        override_plot_properties={'cbar_shrink':1},
        fontsize=fontsize
    )

Integrated effective temperature and the voltage square spectral density

  • this is done in loop 100x

[ ]:
def propagate_errors(rstd_teff, rstd_galactic_map, rstd_antenna_map):
    # calculate the effective temperature
    teff = galactic_map * antenna_map**2

    ratios_teff = np.array([])
    for i in range(100):
        # teff_distored = distort_array(teff, rstd=rstd_teff)
        teff_distored = (
            distort_array(galactic_map, rstd_galactic_map)
            * distort_array(antenna_map, rstd_antenna_map) ** 2
        )
        ratios_teff = np.append(
            ratios_teff, integrate_hpmap(teff) / integrate_hpmap(teff_distored)
        )

    # do the same with voltage square density, this should give the same results as before (more-less)
    # this is to check that there the extra constants have no effect
    ratios_v2sp = np.array([])
    v2sp = voltage_squared_spectral_density(antenna_map, galactic_map, frequency_MHz)

    for i in range(100):
        galactic_map_distorted = distort_array(galactic_map, rstd_galactic_map)
        antenna_map_distorted = distort_array(antenna_map, rstd_antenna_map)

        v2sp_distorted = voltage_squared_spectral_density(
            antenna_map_distorted, galactic_map_distorted, frequency_MHz
        )

        ratios_v2sp = np.append(ratios_v2sp, v2sp / v2sp_distorted)

    # figure
    fig, ax = plt.subplots()
    bins = np.histogram_bin_edges(np.concatenate([ratios_teff, ratios_v2sp]), bins=20)
    ax.set_title(
        "Error on sky temperature: {:.2f}, error on antenna gain: {:.2f}".format(
            rstd_galactic_map, rstd_antenna_map
        )
    )

    ax.hist(
        ratios_teff,
        bins=bins,
        alpha=0.5,
        label=r"T$_{{\mathrm{{eff}}}}$: Mean={:.2f}, Std={:.1e}".format(
            np.mean(ratios_teff), np.std(ratios_teff)
        ),
    )
    ax.hist(
        ratios_v2sp,
        bins=bins,
        alpha=0.5,
        label=r"V$^2_f$: Mean={:.2f}, Std={:.1e}".format(
            np.mean(ratios_v2sp), np.std(ratios_v2sp)
        ),
    )
    ax.legend()

From the histograms bellow we conclude two things:

  • there is virtually no effect of the constants on the relative error propagation (as expected)

  • the error on the antenna gains is causing a bias, i.e. the whole distributions is shifted

[ ]:
propagate_errors(rstd_teff=0.5, rstd_galactic_map=0.3, rstd_antenna_map=0.2)
[ ]:
propagate_errors(rstd_teff=0.5, rstd_galactic_map=0.5, rstd_antenna_map=0.)
[ ]:
propagate_errors(rstd_teff=0.5, rstd_galactic_map=0., rstd_antenna_map=0.25)

Analytical calculations of the error propagation to the integrated sky temperature

We attempt to calculate the following equation other assumption of isotropic effective temperature distorted by 50%

\begin{equation} \rho_{G\mathrm{sky}} (t,f) = \sqrt{ \frac{\sum_{\delta} \sin^2{\delta} }{N_{\alpha} \left( \sum_{\delta} \sin{\delta} \right)^2} } \rho_{T\mathrm{eff}}(t,f,\alpha,\delta) \end{equation}

[ ]:
# check what is the size of the theta and phi when integrating (should be 1000x500)
# so, for half sphere it is 1000x250
# check the value of the propagated error at 250 then
PHI, THETA, grid_map = hpmap2grid(galactic_map)
phi_delta = abs(np.diff(PHI[0, :])[0])
phi_theta = abs(np.diff(THETA[:, 0])[0])
print(PHI[0, :].size, THETA[:, 0].size)
[ ]:
propagated_err_analytical = np.array([])

Ns = [10, 100, 250, 500, 1000, 10000, 50000, 100000]
for N in Ns:
    sum_1 = 0
    sum_2 = 0
    for x in np.linspace(0, np.pi/2, N):
        sum_1 += np.sin(x)**2
        sum_2 += np.sin(x)

    propagated_err_analytical = np.append(propagated_err_analytical, np.sqrt(sum_1/sum_2**2 /(N*4))*rstd_teff)

print(propagated_err_analytical)
[ ]:
fig, ax = plt.subplots()
ax.plot(Ns, np.log10(propagated_err_analytical))
ax.scatter(Ns[2], np.log10(propagated_err_analytical)[2], c='r', label='our used resolution')
ax.set_xlabel('N$_{\\theta}$')
ax.set_ylabel('propagated log$_{10}$($\\rho$)')
ax.legend()

Propagation to the power

[ ]:
# check the propagation to the finally calculated power at a range of LST and frequency values
lst_range = np.asarray(list(range(24))) + 0.5
freq_Mhz_range = range(30, 81, 1)

hw_file_path = "./antenna_setup_files/HardwareProfileList_realistic.xml"
hw_dict = read_hw_file(hw_file_path)
[ ]:
def propagated_errors_2_power(rstd_antenna_map=0.2, rstd_galactic_map=0.3):

    lst_range = np.asarray(list(range(24))) + 0.5
    freq_Mhz_range = range(30, 81, 1)
    impedance_func = hw_dict["IImpedance"][
        "antenna_EW"
    ]

    power_density_DF = calculate_power_spectral_density(
        antenna_inst=antenna_inst,
        galactic_map_inst=galactic_map_inst,
        lst_range=lst_range,
        freq_Mhz_range=freq_Mhz_range,
        latitude=LATITUDE,
        update_antenna_conventions=update_antenna_conventions,
        impedance_func=impedance_func,
        distort_antenna_map=rstd_antenna_map,
        distort_galactic_map=rstd_galactic_map
    )

    power_DF = integrate_spectral_density(power_density_DF, integrated_MHz_bands=power_density_DF.columns.values)

    hw_reponse_1 = dB2PowerAmp(
        hw_dict["RResponse"]["LNA"](power_DF.columns)
    )
    hw_reponse_2 = dB2PowerAmp(
        hw_dict["RResponse"]["digitizer"](power_DF.columns)
    )
    hw_reponse_3 = dB2PowerAmp(
        hw_dict["RResponse"]["cable_fromLNA2digitizer"](power_DF.columns)
    )
    hw_reponse_4 = dB2PowerAmp(
        hw_dict["RResponse"][
            "impedance_matching_EW"
        ](power_DF.columns)
    )


    power_in_HW_DF = power_DF.multiply(
        hw_reponse_1 * hw_reponse_2 * hw_reponse_3 * hw_reponse_4
    )

    # to piko and round
    power_in_HW_DF_distorted = (power_in_HW_DF*1e+12).round(3)


    power_DF_normal = pd.read_csv('./simulated_power_datasets/Salla_EW_GSM16.csv', index_col=0)
    power_DF_normal.columns = power_DF_normal.columns.astype(float)

    ratio_DF = power_DF_normal/power_in_HW_DF_distorted.values

    # figures

    fig = px.imshow(
        ratio_DF.T - 1, width=600, aspect="cube", color_continuous_scale="jet"
    )
    fig.update_layout(
    #    title="<b></b>",
        xaxis=dict(title="<b>LST</b>", tickprefix="<b>", ticksuffix="</b>", dtick=2),
        yaxis=dict(
            title="<b>frequency [MHz]</b>",
            tickprefix="<b>",
            ticksuffix="</b>",
            range=(30, 80),
            tick0=0,
            dtick=10,
            autorange=False,
        ),
        coloraxis=dict(
            colorbar=dict(
                title=dict(
                    text="<b>power [pW]</b>",
                    side="right",
                ),
                tickprefix="<b>",
                ticksuffix="</b>",
            ),
        ),
        font=dict(
            # family=font,
            size=20,
            color="black",
        ),
    )
    fig.update_layout(
        coloraxis=dict(
            colorbar=dict(title=dict(text="<b>zundistorted/distorted power - 1</b>", side="right")),
            # cmin=0,
            # cmax=40,
        )
    )
    fig.show()

    # fig, ax = plt.subplots()
    # ax.plot(ratio_DF.mean(axis=0).values)
    # ax.set_xlabel('frequency [MHz]')

    # fig, ax = plt.subplots()
    # ax.plot(ratio_DF.mean(axis=1).values)
    # ax.set_xlabel('LST')

    # ratios_power = ratio_DF.values.flatten()

    # fig, ax = plt.subplots()

    # ax.hist(
    #     ratios_power,
    #     bins=20,
    #     alpha=0.5,
    #     label=r"Mean={:.2f}, Std={:.1e}".format(
    #         np.mean(ratios_power), np.std(ratios_power)
    #     ),
    # )

    # ax.legend()
    # ax.set_xlabel(r'$\frac{\mathrm{undistorted}}{\mathrm{distorted}}$ power')
    # ax.set_ylabel('entries')
    return ratio_DF
[ ]:
ratio_DF_1 = propagated_errors_2_power(rstd_antenna_map=0.2, rstd_galactic_map=0.3)
[ ]:
ratio_DF_2 = propagated_errors_2_power(rstd_antenna_map=0.25, rstd_galactic_map=0.)
[ ]:
ratio_DF_3 = propagated_errors_2_power(rstd_antenna_map=0., rstd_galactic_map=0.5)
[ ]:
ratio_DF_2aa = propagated_errors_2_power(rstd_antenna_map=0.20, rstd_galactic_map=0.)
[ ]:
ratio_DF_2a = propagated_errors_2_power(rstd_antenna_map=0.15, rstd_galactic_map=0.)
[ ]:
ratio_DF_2b = propagated_errors_2_power(rstd_antenna_map=0.10, rstd_galactic_map=0.)
[ ]:
ratio_DF_2c = propagated_errors_2_power(rstd_antenna_map=0.05, rstd_galactic_map=0.)
[ ]:
fig, ax = plt.subplots()

labels = [
    r"$\rho_{ant}$=0.20, $\rho_{sky}$=0.30",
    r"$\rho_{ant}$=0.00, $\rho_{sky}$=0.50",
    r"$\rho_{ant}$=0.25, $\rho_{sky}$=0.00",
    r"$\rho_{ant}$=0.20, $\rho_{sky}$=0.00",
    r"$\rho_{ant}$=0.15, $\rho_{sky}$=0.00",
    r"$\rho_{ant}$=0.10, $\rho_{sky}$=0.00",
    r"$\rho_{ant}$=0.05, $\rho_{sky}$=0.00",
]

bins = np.histogram_bin_edges(
    np.concatenate(
        [
            ratio_DF_1.values.flatten(),
            ratio_DF_3.values.flatten(),
            ratio_DF_2.values.flatten(),
            ratio_DF_2aa.values.flatten(),
            ratio_DF_2a.values.flatten(),
            ratio_DF_2b.values.flatten(),
            ratio_DF_2c.values.flatten(),
        ]
    ),
    bins=100,
)
bins -= 1
results = []

for i, df in enumerate(
    [
        ratio_DF_1,
        ratio_DF_3,
        ratio_DF_2,
        ratio_DF_2aa,
        ratio_DF_2a,
        ratio_DF_2b,
        ratio_DF_2c,
    ]
):
    ratios_power = df.values.flatten() - 1
    ax.hist(ratios_power, bins=bins, alpha=0.5, label=labels[i])
    stats = r"$\sigma$={:.3f}, $\mu$={:.3f}".format(
        np.abs(np.std(ratios_power)), round(np.mean(ratios_power), 3)+0
    )
    # print(labels[i] +',  '+ stats)
    temp = labels[i] + ",  " + stats
    numbers = re.findall(r"[-+]?\d*\.?\d+(?:[eE][-+]?\d+)?", temp)
    # print(numbers)
    results.append([float(i) for i in numbers])

ax.legend()
ax.set_xlabel(r"$\frac{\mathrm{(undistorted \: power)}}{\mathrm{(distorted \: power)}}$ - 1")
ax.set_ylabel("entries")

s = (
    r"""| $\rho_{{ant}}$ | $\rho_{{sky}}$ | $\rho$  | $b$     |
|----------------|----------------|---------|---------|
| {:.2f}         | {:.2f}         | {:.1e}  | {:.2f}  |
| {:.2f}         | {:.2f}         | {:.1e}  | {:.2f}  |
| {:.2f}         | {:.2f}         | {:.1e}  | {:.2f}  |
| {:.2f}         | {:.2f}         | {:.1e}  | {:.2f}  |
| {:.2f}         | {:.2f}         | {:.1e}  | {:.2f}  |
| {:.2f}         | {:.2f}         | {:.1e}  | {:.2f}  |
| {:.2f}         | {:.2f}         | {:.1e}  | {:.2f}  |"""
).format(*np.asarray(results).flatten())

md(print("RESULTS:"))
md(s)
[106]:
results_df = pd.DataFrame(results, columns=['$\\rho_{ant}$', '$\\rho_{sky}$', '$\\rho$', '$b$'])*100
s = results_df.style.format('{:.1f}').hide(axis='index')
latex_table = s.to_latex()
print(latex_table)
display(s)
\begin{tabular}{rrrr}
$\rho_{ant}$ & $\rho_{sky}$ & $\rho$ & $b$ \\
20.0 & 30.0 & 0.3 & -3.8 \\
0.0 & 50.0 & 0.3 & 0.0 \\
25.0 & 0.0 & 0.3 & -5.9 \\
20.0 & 0.0 & 0.2 & -3.9 \\
15.0 & 0.0 & 0.2 & -2.2 \\
10.0 & 0.0 & 0.1 & -1.0 \\
5.0 & 0.0 & 0.1 & -0.2 \\
\end{tabular}

$\rho_{ant}$ $\rho_{sky}$ $\rho$ $b$
20.0 30.0 0.3 -3.8
0.0 50.0 0.3 0.0
25.0 0.0 0.3 -5.9
20.0 0.0 0.2 -3.9
15.0 0.0 0.2 -2.2
10.0 0.0 0.1 -1.0
5.0 0.0 0.1 -0.2

From the distributions above we see that:

  • the error on the antenna gain is indeed causing a bias

  • the propagated error is of order of \(10^{-3}\) just like was analytically calculated

And from the color maps we learned that there is no LST or frequency correlation after propagating the error to the power.