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)
[ ]:
# 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%
[ ]:
# 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}
| 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
And from the color maps we learned that there is no LST or frequency correlation after propagating the error to the power.