Harrison 3063/#36 - PMCL - Maps spx: Use reduced chi-squared

Divide by degrees of freedom (which is actually #points - #parameters) to get the reduced chi-squared

Author: pleasant-menlo

imap_l3_processing/cdf/config/imap_hi_l3_variable_attrs.yaml

@@ -431,7 +431,7 @@ ena_spectral_scalar_stat_uncert:
 ena_spectral_index_chisq:
   NAME: ena_spectral_index_chisq
   DATA_TYPE: CDF_REAL4
-  CATDESC: ENA spectral index chi-squared
+  CATDESC: ENA spectral index reduced chi-squared
   DEPEND_0: epoch
   DEPEND_1: energy
   DEPEND_2: longitude

imap_l3_processing/cdf/config/imap_lo_l3_variable_attrs.yaml

@@ -694,7 +694,7 @@ ena_spectral_scalar_stat_uncert:
 ena_spectral_index_chisq:
   NAME: ena_spectral_index_chisq
   DATA_TYPE: CDF_REAL4
-  CATDESC: ENA spectral index chi-squared
+  CATDESC: ENA spectral index reduced chi-squared
   DEPEND_0: epoch
   DEPEND_1: energy
   DEPEND_2: longitude

imap_l3_processing/cdf/config/imap_ultra_l3_variable_attrs.yaml

@@ -431,7 +431,7 @@ ena_spectral_scalar_stat_uncert:
 ena_spectral_index_chisq:
   NAME: ena_spectral_index_chisq
   DATA_TYPE: CDF_REAL4
-  CATDESC: ENA spectral index chi-squared
+  CATDESC: ENA spectral index reduced chi-squared
   DEPEND_0: epoch
   DEPEND_1: energy
   DEPEND_2: longitude

imap_l3_processing/maps/spectral_fit.py

@@ -167,7 +167,7 @@ def fit_arrays_to_power_law(fluxes: np.ndarray, uncertainties: np.ndarray, energ
                     scalar_coefficients[i] = a
                     gamma_errors[i] = gamma_error
                     scalar_coefficient_errors[i] = a_error
-                    chisqs[i] = fit.fnorm
+                    chisqs[i] = fit.fnorm / fit.dof
         output_gammas[epoch, 0] = gammas.reshape(fluxes.shape[2:])
         output_gamma_errors[epoch, 0] = gamma_errors.reshape(fluxes.shape[2:])
         output_scalar_coefficients[epoch, 0] = scalar_coefficients.reshape(fluxes.shape[2:])

tests/maps/test_spectral_fit.py

@@ -55,9 +55,9 @@ def test_finds_best_fit(self):
                 np.testing.assert_array_almost_equal( result_chi_square ,np.zeros((1,1, *spacial_dimension_shape)))
 
     def test_finds_chi_squared_for_fit(self):
-        energies = np.array([1, 10, 100])
-        flux_data = np.array([99, 10, 1])
-        errors = np.array([1,1, 1e-3])
+        energies = np.array([1, 10, 30, 100])
+        flux_data = np.array([99, 10, 3, 1])
+        errors = np.array([1,1, 0.1, 1e-3])
 
         cases = [
             ("rectangular", (1, 1)),
@@ -73,7 +73,8 @@ def test_finds_chi_squared_for_fit(self):
                 fitted_flux = result_A * np.power(energies.reshape((-1,*spacial_dimension_shape)), -result_gamma)
                 residual = fitted_flux - flux
                 chisquared = np.sum(np.square(residual/uncertainty))
-                np.testing.assert_array_almost_equal(result_chi_square, chisquared)
+                reduced_chisquared = chisquared / (len(energies) - 2)
+                np.testing.assert_array_almost_equal(result_chi_square, reduced_chisquared)
 
     def test_spectral_fit_map(self):
         epoch = np.array([datetime.now()])
@@ -332,8 +333,8 @@ def test_finds_best_fit_with_ibex_data(self):
                                                                               [52.494673, 44.64706]]]]))
         np.testing.assert_array_almost_equal(scalar_coefficient_errors, np.array([[[[32.245348, 26.075043],
          [19.584518, 20.080801]]]]))
-        np.testing.assert_array_almost_equal(chisq, np.array([[[[1.33863 , 1.174914],
-         [1.034745, 1.426575]]]]))
+        np.testing.assert_array_almost_equal(chisq, np.array([[[[0.44621 , 0.391638],
+         [0.344915, 0.475525]]]]))
 
     def test_finds_best_fit_with_zeros_in_flux_and_not_uncertainty(self):
         energies = np.geomspace(1, 1e10, 23)
@@ -354,7 +355,7 @@ def test_finds_best_fit_with_zeros_in_flux_and_not_uncertainty(self):
                                                      np.array([1.472697]).reshape(1, 1, *spacial_dimension_shape))
                 np.testing.assert_array_almost_equal(scalar_coefficients,
                                                      np.array([1.251819]).reshape(1, 1, *spacial_dimension_shape))
-                np.testing.assert_array_almost_equal(chisq, np.array([46.587911]).reshape(1, 1, *spacial_dimension_shape))
+                np.testing.assert_array_almost_equal(chisq, np.array([2.218472]).reshape(1, 1, *spacial_dimension_shape))
 
     def test_returns_nan_when_only_one_point_is_valid(self):
         energies = np.geomspace(1, 1e10, 5)