From 103bde61f821cff05280cc4ff07980458194f786 Mon Sep 17 00:00:00 2001 From: Eric Teunis de Boone Date: Thu, 13 Apr 2023 12:34:14 +0200 Subject: [PATCH] ZH: (WIP) fitting random phasor sum to antenna phases --- .../bd_antenna_phase_deltas.py | 2 +- airshower_beacon_simulation/lib/figlib.py | 88 ++++++++++++++++--- 2 files changed, 78 insertions(+), 12 deletions(-) diff --git a/airshower_beacon_simulation/bd_antenna_phase_deltas.py b/airshower_beacon_simulation/bd_antenna_phase_deltas.py index 62d4fb6..4b54b85 100755 --- a/airshower_beacon_simulation/bd_antenna_phase_deltas.py +++ b/airshower_beacon_simulation/bd_antenna_phase_deltas.py @@ -206,7 +206,7 @@ if __name__ == "__main__": figsize=figsize, hist_kwargs=hist_kwargs, fit_gaussian=plot_residuals, - fit_ricianphase=plot_residuals, + fit_randomphasesum=plot_residuals, return_fit_info = True, ) diff --git a/airshower_beacon_simulation/lib/figlib.py b/airshower_beacon_simulation/lib/figlib.py index 36f1c61..81a943a 100644 --- a/airshower_beacon_simulation/lib/figlib.py +++ b/airshower_beacon_simulation/lib/figlib.py @@ -1,9 +1,35 @@ import matplotlib.pyplot as plt import numpy as np -import scipy.stats as stats +from scipy import stats +from scipy import special +from scipy import optimize from itertools import zip_longest +def expectation(x,pdfx): + dx = x[1]-x[0] + return np.sum(x*pdfx*dx) + +def variance(x,pdfx): + mu = expectation(x,pdfx) + dx = x[1]-x[0] + return np.sum((x**2*pdfx*dx))-mu**2 + +def random_phase_sum_distribution(theta, sigma, s=1): + theta = np.asarray(theta) + ct = np.cos(theta) + st = np.sin(theta) + k = s/sigma + pipi = 2*np.pi + return (np.exp(-k**2/2)/pipi) + ( + (pipi**-0.5)*k*np.exp(-(k*st)**2/2)) * ( + (1.+special.erf(k*ct*2**-0.5))*ct/2) + +def gaussian_phase_distribution(theta, sigma, s=1): + theta = np.asarray(theta) + k=s/sigma + return (2*np.pi)**-0.5*k*np.exp(-(k*theta)**2/2) + def phase_comparison_figure( measured_phases, true_phases, @@ -16,6 +42,7 @@ def phase_comparison_figure( legend_on_scatter=True, secondary_axis='time', fit_gaussian=False, + fit_randomphasesum=False, mean_snr=None, return_fit_info=False, **fig_kwargs @@ -29,14 +56,19 @@ def phase_comparison_figure( default_text_kwargs = dict(fontsize=14, verticalalignment='top') default_sc_kwargs = dict(alpha=0.6, ls='none') + do_hist_plot = hist_kwargs is not False + if hist_kwargs is False: + hist_kwargs = {} + + do_scatter_plot = sc_kwargs is not False + if sc_kwargs is False: + sc_kwargs = {} + fig_kwargs = {**default_fig_kwargs, **fig_kwargs} hist_kwargs = {**default_hist_kwargs, **hist_kwargs} text_kwargs = {**default_text_kwargs, **text_kwargs} sc_kwargs = {**default_sc_kwargs, **sc_kwargs} - do_hist_plot = hist_kwargs is not False - do_scatter_plot = sc_kwargs is not False - fig, axs = plt.subplots(0+do_hist_plot+do_scatter_plot, 1, **fig_kwargs) if not hasattr(axs, '__len__'): @@ -67,10 +99,14 @@ def phase_comparison_figure( text_kwargs=text_kwargs, hist_kwargs={**hist_kwargs, **dict(label='Measured', color=colors[0], ls='solid')}, mean_snr=mean_snr, + fit_distr=[], ) if fit_gaussian: - this_kwargs['fit_distr'] = 'gaussian' + this_kwargs['fit_distr'].append('gaussian') + + if fit_randomphasesum: + this_kwargs['fit_distr'].append('randomphasesum') _, fit_info = fitted_histogram_figure( measured_phases, @@ -126,7 +162,7 @@ def fitted_histogram_figure( text_kwargs = {**default_text_kwargs, **text_kwargs} if ax is None: - fig, ax = plt.subplots(1,1, **fig_kwargs) + fig, ax = plt.subplots(1, 1, **fig_kwargs) else: fig = ax.get_figure() @@ -139,6 +175,8 @@ def fitted_histogram_figure( min_x = min(amplitudes) max_x = max(amplitudes) + bin_centers = bins[:-1] + np.diff(bins) / 2 + dx = bins[1] - bins[0] scale = len(amplitudes) * dx @@ -146,6 +184,9 @@ def fitted_histogram_figure( for distr in fit_distr: fit_params2text_params = lambda x: x + fit_ys = None + fit_params = None + cdf = None if 'rice' == distr: name = "Rice" @@ -166,19 +207,44 @@ def fitted_histogram_figure( fit_params2text_params = lambda x: (x[0]+x[1]/2,) + elif 'randomphasesum' == distr: + name = "RandPhaseS" + param_names = [ "$\\sigma$", 's'] + pdf = random_phase_sum_distribution + + bounds = ((0,0.9999), (np.inf,1)) + fit_params, pcov = optimize.curve_fit(pdf, bin_centers, counts, bounds=bounds) + fit_ys = pdf( xs, *fit_params) + + fit_params2text_params = lambda x: (x[1], x[0]) + + elif 'gaussphase' == distr: + name = 'GaussPhase' + param_names = [ "$\\sigma$", 's'] + pdf = gaussian_phase_distribution + + + bounds = ((0,0.9999), (np.inf,1)) + fit_params, pcov = optimize.curve_fit(pdf, bin_centers, counts, bounds=bounds) + fit_ys = pdf( xs, *fit_params) + + fit_params2text_params = lambda x: (x[1], x[0]) + else: raise ValueError('Unknown distribution function '+distr) label = name +"(" + ','.join(param_names) + ')' - fit_params = distr_func.fit(amplitudes) - fit_ys = distr_func.pdf(xs, *fit_params) + if fit_ys is None: + fit_params = distr_func.fit(amplitudes) + fit_ys = scale * distr_func.pdf(xs, *fit_params) + cdf = distr_func.cdf - ax.plot(xs, fit_ys*scale, label=label) + ax.plot(xs, fit_ys, label=label) chisq_strs = [] - if calc_chisq: - ct = np.diff(distr_func.cdf(bins, *fit_params))*np.sum(counts) + if calc_chisq and cdf: + ct = np.diff(cdf(bins, *fit_params))*np.sum(counts) c2t = stats.chisquare(counts, ct, ddof=len(fit_params)) chisq_strs = [ f"$\\chi^2$/dof = {c2t[0]: .2g}/{len(fit_params)}"