calc_sbp_excess.py: add function "estimate_excess_error()" using Monte Carlo

1 files changed, 123 insertions, 33 deletions

diff --git a/calc_sbp_excess.py b/calc_sbp_excess.py
index 2b472ca..6bf91b2 100755
--- a/calc_sbp_excess.py
+++ b/calc_sbp_excess.py
@@ -3,11 +3,14 @@
 #
 # Aaron LI
 # Created: 2016-04-26
-# Updated: 2016-04-26
+# Updated: 2016-05-06
+#
+# Changelog:
+# 2016-05-06:
+#   * Add function `estimate_excess_error()` to estimate uncertainty
+#   * update according to `sbp_fit` renamed to `fit_sbp`
+#   * PEP8 fixes
 #
-# TODO:
-#   * to estimate errors for the excess value and ratio (Monte Carlo:
-#     repeatedly fit the randomized SBP data and calculate excess)
 #
 
 """
@@ -19,9 +22,6 @@ NOTE:
 * excess ratio: excess_value / brightness_observed
 """
 
-__version__ = "0.1.0"
-__date__    = "2016-04-26"
-
 
 import sys
 import argparse
@@ -31,51 +31,125 @@ from collections import OrderedDict
 import numpy as np
 from configobj import ConfigObj
 
-from sbp_fit import make_model
+from fit_sbp import make_model, make_sbpfit
 
+__version__ = "0.2.0"
+__date__ = "2016-05-06"
 
-def calc_excess(data, model):
+
+def calc_excess(data, fitted_model):
     """
     Calculate the central brightness excess value and ratio with respect
     to the fitted SBP model (single-beta).
 
+    TODO/XXX:
+      * whether to interpolate the SBP?
+
     Arguments:
       * data: raw 4-column SBP data
-      * model: fitted SBP model
+      * fitted_model: fitted SBP model
     """
-    radii            = data[:, 0]
-    radii_err        = data[:, 1]
-    brightness       = data[:, 2]
-    brightness_model = model.f(radii)
-    rin   = radii - radii_err
-    rout  = radii + radii_err
+    radii = data[:, 0]
+    radii_err = data[:, 1]
+    brightness = data[:, 2]
+    brightness_model = fitted_model.f(radii)
+    rin = radii - radii_err
+    rout = radii + radii_err
     areas = np.pi * (rout**2 - rin**2)
-    bsum_obs   = np.sum(brightness * areas)
+    bsum_obs = np.sum(brightness * areas)
     bsum_model = np.sum(brightness_model * areas)
     excess_value = bsum_obs - bsum_model
     excess_ratio = excess_value / bsum_obs
     excess = {
-            "brightness_obs":   bsum_obs,
-            "brightness_model": bsum_model,
-            "excess_value":     excess_value,
-            "excess_ratio":     excess_ratio,
+        "brightness_obs":   bsum_obs,
+        "brightness_model": bsum_model,
+        "excess_value":     excess_value,
+        "excess_ratio":     excess_ratio,
     }
     return excess
 
 
+def estimate_excess_error(data, sbpfit, mctimes, verbose=False):
+    """
+    Estimate the uncertainty of central excess value by Monte Carlo method.
+
+    XXX/TODO:
+      * whether also consider the uncertainty of R500?
+
+    Arguments:
+      * data: 4-column SBP data (radius, r_err, brightness, brightness_err)
+      * sbpfit: `SbpFit` object used to perform SBP fitting
+      * mctimes: number of Monte Carlo iterations
+    """
+    brightness = data[:, 2]
+    brightness_err = data[:, 3]
+    params = sbpfit.dump_params()
+    ev_results = []
+    er_results = []
+    if verbose:
+        print("Estimating excess uncertainty by Monte Carlo " +
+              "(%d times) ..." % mctimes, end="", flush=True)
+    for i in range(mctimes):
+        if verbose and i % 100 == 0:
+            print("%d..." % i, end="", flush=True)
+        # randomize SBP data
+        brightness_rand = np.random.normal(brightness, scale=brightness_err)
+        sbpdata_rand = data.copy()
+        sbpdata_rand[:, 2] = brightness_rand
+        # load randomized data and perform SBP fit
+        sbpfit.reset(keep_source=True)
+        sbpfit.load_data(sbpdata_rand, keep_mask=True)
+        sbpfit.load_params(params)
+        sbpfit.fit()
+        model_rand = sbpfit.get_model()
+        excess = calc_excess(data=sbpdata_rand, fitted_model=model_rand)
+        ev_results.append(excess["excess_value"])
+        er_results.append(excess["excess_ratio"])
+    if verbose:
+        print("DONE!", flush=True)
+    ev_mean = np.mean(ev_results)
+    ev_std = np.std(ev_results)
+    ev_q25, ev_median, ev_q75 = np.percentile(ev_results, q=(25, 50, 75))
+    er_mean = np.mean(er_results)
+    er_std = np.std(er_results)
+    er_q25, er_median, er_q75 = np.percentile(er_results, q=(25, 50, 75))
+    results = {
+        "excess_value_mean":   ev_mean,
+        "excess_value_median": ev_median,
+        "excess_value_q25":    ev_q25,
+        "excess_value_q75":    ev_q75,
+        "excess_value_std":    ev_std,
+        "excess_ratio_mean":   er_mean,
+        "excess_ratio_median": er_median,
+        "excess_ratio_q25":    er_q25,
+        "excess_ratio_q75":    er_q75,
+        "excess_ratio_std":    er_std,
+    }
+    return results
+
+
 def main():
     # default arguments
     default_outfile = "excess.json"
+    default_mctimes = 1000
 
     parser = argparse.ArgumentParser(
             description="Calculate the central brightness excess value",
             epilog="Version: %s (%s)" % (__version__, __date__))
     parser.add_argument("-V", "--version", action="version",
-            version="%(prog)s " + "%s (%s)" % (__version__, __date__))
+                        version="%(prog)s " +
+                                "%s (%s)" % (__version__, __date__))
+    parser.add_argument("-v", "--verbose", dest="verbose",
+                        required=False, action="store_true",
+                        help="show verbose information")
+    parser.add_argument("-m", "--mctimes", dest="mctimes", required=False,
+                        type=int, default=default_mctimes,
+                        help="number of MC times to estimate excess error " +
+                             "(default: %d)" % default_mctimes)
     parser.add_argument("config", help="Config file for SBP fitting")
     parser.add_argument("outfile", nargs="?", default=default_outfile,
-            help="Output json file to save the results " + \
-                 "(default: %s)" % default_outfile)
+                        help="Output json file to save the results " +
+                             "(default: %s)" % default_outfile)
     args = parser.parse_args()
 
     config = ConfigObj(args.config)
@@ -83,7 +157,7 @@ def main():
     config_model = config[modelname]
 
     model = make_model(config, modelname=modelname)
-    print("SBP model: %s" % model.name, file=sys.stderr)
+    print("SBP model: %s" % model.long_name, file=sys.stderr)
 
     sbpfit_outfile = config.get("outfile")
     sbpfit_outfile = config_model.get("outfile", sbpfit_outfile)
@@ -94,17 +168,33 @@ def main():
     for pname, pvalue in sbpfit_results["params"].items():
         model.set_param(name=pname, value=pvalue[0])
 
+    sbpfit = make_sbpfit(config, model=model)
+
     sbpdata = np.loadtxt(config["sbpfile"])
-    excess = calc_excess(data=sbpdata, model=model)
+    excess = calc_excess(data=sbpdata, fitted_model=model)
+    excess_err = estimate_excess_error(data=sbpdata, sbpfit=sbpfit,
+                                       mctimes=args.mctimes,
+                                       verbose=args.verbose)
 
     excess_data = OrderedDict([
-            ("name",             config["name"]),
-            ("obsid",            int(config["obsid"])),
-            ("model",            modelname),
-            ("brightness_obs",   excess["brightness_obs"]),
-            ("brightness_model", excess["brightness_model"]),
-            ("excess_value",     excess["excess_value"]),
-            ("excess_ratio",     excess["excess_ratio"]),
+        ("name",                config["name"]),
+        ("obsid",               int(config["obsid"])),
+        ("model",               modelname),
+        ("brightness_obs",      excess["brightness_obs"]),
+        ("brightness_model",    excess["brightness_model"]),
+        ("excess_value",        excess["excess_value"]),
+        ("excess_value_mean",   excess_err["excess_value_mean"]),
+        ("excess_value_median", excess_err["excess_value_median"]),
+        ("excess_value_q25",    excess_err["excess_value_q25"]),
+        ("excess_value_q75",    excess_err["excess_value_q75"]),
+        ("excess_value_std",    excess_err["excess_value_std"]),
+        ("excess_ratio",        excess["excess_ratio"]),
+        ("excess_ratio_mean",   excess_err["excess_ratio_mean"]),
+        ("excess_ratio_median", excess_err["excess_ratio_median"]),
+        ("excess_ratio_q25",    excess_err["excess_ratio_q25"]),
+        ("excess_ratio_q75",    excess_err["excess_ratio_q75"]),
+        ("excess_ratio_std",    excess_err["excess_ratio_std"]),
+        ("mc_times",            args.mctimes),
     ])
     excess_data_json = json.dumps(excess_data, indent=2)
     print(excess_data_json)