Calculate Akaike and Bayesian information criteria (AIC and BIC) for …

…the single and double Sersic fits.

docs/notebooks/doublesersic.ipynb

@@ -200,6 +200,8 @@
     "print('sersic_ellip =', morph.sersic_ellip)\n",
     "print('sersic_theta =', morph.sersic_theta)\n",
     "print('sersic_chi2_dof =', morph.sersic_chi2_dof)\n",
+    "print('sersic_aic =', morph.sersic_aic)\n",
+    "print('sersic_bic =', morph.sersic_bic)\n",
     "print()\n",
     "print('DOUBLE SERSIC')\n",
     "print('doublesersic_xc =', morph.doublesersic_xc)\n",
@@ -215,6 +217,8 @@
     "print('doublesersic_ellip2 =', morph.doublesersic_ellip2)\n",
     "print('doublesersic_theta2 =', morph.doublesersic_theta2)\n",
     "print('doublesersic_chi2_dof =', morph.doublesersic_chi2_dof)\n",
+    "print('doublesersic_aic =', morph.doublesersic_aic)\n",
+    "print('doublesersic_bic =', morph.doublesersic_bic)\n",
     "print()\n",
     "print('OTHER')\n",
     "print('sky_mean =', morph.sky_mean)\n",
@@ -231,6 +235,8 @@
    "source": [
     "Clearly, the fitted double Sersic model is consistent with the \"true\" light distribution that we originally created (n1 = 5, n2 = 1, etc.) and the reduced chi-squared statistic (doublesersic_chi2_dof) is close to 1, indicating a good fit without overfitting. On the other hand, the *single* Sersic fit has a reduced chi-squared statistic much larger than 1, indicating a poor fit (as expected).\n",
     "\n",
+    "We also calculate the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) for the two models, which again favor the *double* Sersic model as the statistically preferred one, since it returns much lower AIC and BIC values.\n",
+    "\n",
     "Also note that statmorph now returns an additional quality flag:\n",
     "\n",
     "- ``flag_doublesersic`` : indicates the quality of the double Sersic fit. Like ``flag`` and ``flag_sersic``, it can take the following values: 0 (good), 1 (suspect), 2 (bad), and 4 (catastrophic)."

statmorph/statmorph.py

@@ -15,10 +15,11 @@
 import skimage.feature
 import skimage.segmentation
 from astropy.utils import lazyproperty
-from astropy.stats import sigma_clipped_stats
+from astropy.stats import (sigma_clipped_stats, akaike_info_criterion_lsq,
+                           bayesian_info_criterion_lsq)
 from astropy.modeling import models, fitting, Fittable2DModel, Parameter
-from astropy.utils.exceptions import (
-    AstropyUserWarning, AstropyDeprecationWarning)
+from astropy.utils.exceptions import (AstropyUserWarning,
+                                      AstropyDeprecationWarning)
 from astropy.convolution import convolve
 import photutils
 
@@ -78,6 +79,8 @@
     'sersic_ellip',
     'sersic_theta',
     'sersic_chi2_dof',
+    'sersic_aic',
+    'sersic_bic',
     'sky_mean',
     'sky_median',
     'sky_sigma',
@@ -104,6 +107,8 @@
     'doublesersic_ellip2',
     'doublesersic_theta2',
     'doublesersic_chi2_dof',
+    'doublesersic_aic',
+    'doublesersic_bic',
 ]
 
 
@@ -2614,17 +2619,18 @@ def _sersic_model(self):
         # The sky background noise is already included in the weightmap:
         fit_weights[locs] = 1.0 / weightmap[locs]
         # Number of "valid" pixels
-        num_validpixels = np.sum(locs)
+        self._sersic_num_validpixels = np.sum(locs)
 
         # Calculate number of free parameters and degrees of freedom
-        num_freeparam = sersic_init.parameters.size  # 7 for Sersic2D
+        self._sersic_num_freeparam = sersic_init.parameters.size  # 7
         for kwarg in ['fixed', 'tied']:
             if kwarg in self._sersic_model_args:
                 for param, value in self._sersic_model_args[kwarg].items():
                     if value:
-                        num_freeparam -= 1
-        assert num_freeparam >= 0
-        self._sersic_num_dof = num_validpixels - num_freeparam
+                        self._sersic_num_freeparam -= 1
+        assert self._sersic_num_freeparam >= 0
+        self._sersic_num_dof = (self._sersic_num_validpixels
+                                - self._sersic_num_freeparam)
         if self._sersic_num_dof <= 0:
             warnings.warn('[sersic] Not enough data for fit.',
                           AstropyUserWarning)
@@ -2746,6 +2752,32 @@ def sersic_chi2_dof(self):
 
         return self._sersic_chi2 / self._sersic_num_dof
 
+    @lazyproperty
+    def sersic_aic(self):
+        """
+        Akaike Information Criterion (AIC) of the fitted model.
+        """
+        _ = self._sersic_model
+        if self._sersic_chi2 == -99.0:
+            return -99.0
+
+        return akaike_info_criterion_lsq(
+            self._sersic_chi2, self._sersic_num_freeparam,
+            self._sersic_num_validpixels)
+
+    @lazyproperty
+    def sersic_bic(self):
+        """
+        Bayesian Information Criterion (BIC) of the fitted model.
+        """
+        _ = self._sersic_model
+        if self._sersic_chi2 == -99.0:
+            return -99.0
+
+        return bayesian_info_criterion_lsq(
+            self._sersic_chi2, self._sersic_num_freeparam,
+            self._sersic_num_validpixels)
+
     ###########################
     # DOUBLE SERSIC MODEL FIT #
     ###########################
@@ -2919,17 +2951,18 @@ def _doublesersic_model(self):
         # The sky background noise is already included in the weightmap:
         fit_weights[locs] = 1.0 / weightmap[locs]
         # Number of "valid" pixels
-        num_validpixels = np.sum(locs)
+        self._doublesersic_num_validpixels = np.sum(locs)
 
         # Calculate number of free parameters and degrees of freedom
-        num_freeparam = doublesersic_init.parameters.size  # 12 for DoubleSersic2D
+        self._doublesersic_num_freeparam = doublesersic_init.parameters.size  # 12
         for kwarg in ['fixed', 'tied']:
             if kwarg in self._doublesersic_model_args:
                 for param, value in self._doublesersic_model_args[kwarg].items():
                     if value:
-                        num_freeparam -= 1
-        assert num_freeparam >= 0
-        self._doublesersic_num_dof = num_validpixels - num_freeparam
+                        self._doublesersic_num_freeparam -= 1
+        assert self._doublesersic_num_freeparam >= 0
+        self._doublesersic_num_dof = (self._doublesersic_num_validpixels
+                                      - self._doublesersic_num_freeparam)
         if self._doublesersic_num_dof <= 0:
             warnings.warn('[doublesersic] Not enough data for fit.',
                           AstropyUserWarning)
@@ -3096,6 +3129,32 @@ def doublesersic_chi2_dof(self):
 
         return self._doublesersic_chi2 / self._doublesersic_num_dof
 
+    @lazyproperty
+    def doublesersic_aic(self):
+        """
+        Akaike Information Criterion (AIC) of the fitted model.
+        """
+        _ = self._doublesersic_model
+        if self._doublesersic_chi2 == -99.0:
+            return -99.0
+
+        return akaike_info_criterion_lsq(
+            self._doublesersic_chi2, self._doublesersic_num_freeparam,
+            self._doublesersic_num_validpixels)
+
+    @lazyproperty
+    def doublesersic_bic(self):
+        """
+        Bayesian Information Criterion (BIC) of the fitted model.
+        """
+        _ = self._doublesersic_model
+        if self._doublesersic_chi2 == -99.0:
+            return -99.0
+
+        return bayesian_info_criterion_lsq(
+            self._doublesersic_chi2, self._doublesersic_num_freeparam,
+            self._doublesersic_num_validpixels)
+
 
 def source_morphology(image, segmap, **kwargs):
     """

statmorph/tests/test_statmorph.py

@@ -328,6 +328,8 @@ def setup_class(self):
             'sersic_ellip': 0.05083866217150,
             'sersic_theta': 2.47831542907976,
             'sersic_chi2_dof': 1.3376238276749,
+            'sersic_aic': 7640.9506975345,
+            'sersic_bic': 7698.176779495,
             'doublesersic_xc': 81.833668324998,
             'doublesersic_yc': 80.569118306697,
             'doublesersic_amplitude1': 538.85377533706,
@@ -341,6 +343,8 @@ def setup_class(self):
             'doublesersic_ellip2': 0.06669592350554,
             'doublesersic_theta2': 2.4997429173919,
             'doublesersic_chi2_dof': 1.1574114573895,
+            'doublesersic_aic': 3848.3566991179,
+            'doublesersic_bic': 3946.4585539074,
             'sky_mean': 3.48760604858398,
             'sky_median': -2.68543863296509,
             'sky_sigma': 150.91754150390625,