ZH: (WIP) fitting random phasor sum to antenna phases

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,
                     )
 

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)}"