Merge pull request #57 from pcubillos/release

arxiv release

MCcubed/VERSION.py

@@ -4,4 +4,4 @@
 # MC3 Version:
 MC3_VER =  2  # Major version
 MC3_MIN =  2  # Minor version
-MC3_REV = 14  # Revision
+MC3_REV = 15  # Revision

MCcubed/plots/mcplots.py

@@ -15,7 +15,7 @@
 sys.path.append(os.path.dirname(os.path.realpath(__file__)) + '/../lib')
 import binarray as ba
 
-def trace(posterior, Zchain=None, title=None, parname=None, thinning=1,
+def trace(posterior, Zchain=None, parname=None, thinning=1,
           burnin=0, fignum=-10, savefile=None, fmt="."):
   """
   Plot parameter trace MCMC sampling
@@ -26,8 +26,6 @@ def trace(posterior, Zchain=None, title=None, parname=None, thinning=1,
      An MCMC posterior sampling with dimension: [nsamples, npars].
   Zchain: 1D integer ndarray
      the chain index for each posterior sample.
-  title: String
-     Plot title.
   parname: Iterable (strings)
      List of label names for parameters.  If None use ['P0', 'P1', ...].
   thinning: Integer
@@ -78,8 +76,6 @@ def trace(posterior, Zchain=None, title=None, parname=None, thinning=1,
   # Make the trace plot:
   plt.figure(fignum, figsize=(8,8))
   plt.clf()
-  if title is not None:
-    plt.suptitle(title, size=16)
 
   plt.subplots_adjust(left=0.15, right=0.95, bottom=0.10, top=0.90,
                       hspace=0.15)
@@ -104,7 +100,7 @@ def trace(posterior, Zchain=None, title=None, parname=None, thinning=1,
     plt.savefig(savefile)
 
 
-def pairwise(posterior, title=None, parname=None, thinning=1,
+def pairwise(posterior, parname=None, thinning=1,
              fignum=-11, savefile=None, nbins=35, nlevels=20,
              absolute_dens=False):
   """
@@ -114,8 +110,6 @@ def pairwise(posterior, title=None, parname=None, thinning=1,
   ----------
   posterior: 2D ndarray
      An MCMC posterior sampling with dimension: [nsamples, nparameters].
-  title: String
-     Plot title.
   parname: Iterable (strings)
      List of label names for parameters.  If None use ['P0', 'P1', ...].
   thinning: Integer
@@ -172,8 +166,6 @@ def pairwise(posterior, title=None, parname=None, thinning=1,
 
   fig = plt.figure(fignum, figsize=(8,8))
   plt.clf()
-  if title is not None:
-    plt.suptitle(title, size=16)
 
   # Plot:
   h = 1 # Subplot index
@@ -222,7 +214,7 @@ def pairwise(posterior, title=None, parname=None, thinning=1,
     plt.savefig(savefile)
 
 
-def histogram(posterior, title=None, parname=None, thinning=1, fignum=-12,
+def histogram(posterior, parname=None, thinning=1, fignum=-12,
                savefile=None, percentile=None, pdf=None, xpdf=None):
   """
   Plot parameter marginal posterior distributions
@@ -232,8 +224,6 @@ def histogram(posterior, title=None, parname=None, thinning=1, fignum=-12,
   posterior: 1D or 2D float ndarray
      An MCMC posterior sampling with dimension [nsamples] or
      [nsamples, nparameters].
-  title: String
-     Plot title.
   parname: Iterable (strings)
      List of label names for parameters.  If None use ['P0', 'P1', ...].
   thinning: Integer
@@ -294,25 +284,17 @@ def histogram(posterior, title=None, parname=None, thinning=1, fignum=-12,
     ncolumns = (npars+2)/3 + (npars+2)%3  # (Trust me!)
 
   histheight = np.amin((2 + 2*(nrows), 8))
-  if nrows == 1:
-    bottom = 0.25
-  else:
-    bottom = 0.15
-
   plt.figure(fignum, figsize=(8, histheight))
   plt.clf()
-  plt.subplots_adjust(left=0.1, right=0.95, bottom=bottom, top=0.9,
-                      hspace=0.4, wspace=0.1)
-
-  if title is not None:
-    a = plt.suptitle(title, size=16)
+  plt.subplots_adjust(left=0.1, right=0.95, bottom=0.18, top=0.95,
+                      hspace=0.55, wspace=0.1)
 
   maxylim = 0  # Max Y limit
   for i in np.arange(npars):
     ax = plt.subplot(nrows, ncolumns, i+1)
-    a  = plt.xticks(size=fs, rotation=90)
+    a  = plt.xticks(size=fs-1.5, rotation=90)
     if i%ncolumns == 0:
-      a = plt.yticks(size=fs)
+      a = plt.yticks(size=fs-1.5)
     else:
       a = plt.yticks(visible=False)
     plt.xlabel(parname[i], size=fs)
@@ -436,7 +418,7 @@ def RMS(binsz, rms, stderr, rmslo, rmshi, cadence=None, binstep=1,
     plt.savefig(savefile)
 
 
-def modelfit(data, uncert, indparams, model, nbins=75, title=None,
+def modelfit(data, uncert, indparams, model, nbins=75,
              fignum=-22, savefile=None, fmt="."):
   """
   Plot the binned dataset with given uncertainties and model curves
@@ -455,8 +437,6 @@ def modelfit(data, uncert, indparams, model, nbins=75, title=None,
     Model of data.
   nbins:  Integer
     Number of bins in the output plot.
-  title:  String
-    If not None, plot title.
   fignum:  Integer
     The figure number.
   savefile:  Boolean
@@ -485,8 +465,6 @@ def modelfit(data, uncert, indparams, model, nbins=75, title=None,
 
   # Data and Model:
   a = plt.axes([0.15, 0.35, 0.8, 0.55])
-  if title is not None:
-    p = plt.title(title, size=fs)
   p = plt.errorbar(binindp, bindata, binuncert, fmt='ko', ms=4,
                    label='Binned Data')
   p = plt.plot(indparams, model, "b", lw=2, label='Best Fit')

MCcubed/rednoise/prayer.py

@@ -14,7 +14,7 @@
 from .. import fit   as mf
 
 
-def prayer(configfile, nprays=0, savefile=None):
+def prayer(configfile=None, nprays=0, savefile=None):
   """
   Implement a prayer-bead method to estimate parameter uncertainties.
 
@@ -34,6 +34,11 @@ def prayer(configfile, nprays=0, savefile=None):
   for god's sake!
   """
 
+  print("Believing in a prayer bead is a mere act of faith, "
+        "please don't use it\nfor published articles (Cubillos et al. 2016).")
+  return None
+
+  # Here's the code.
   config = ConfigParser.SafeConfigParser()
   config.read([configfile])
   cfgsec = "MCMC"

docs/getstarted.rst

@@ -66,8 +66,8 @@ binaries, run:
      import mccubed as mc3
      help(mc3.mcmc)
 
-Example 1 (Interactive)
------------------------
+Example 1 Interactive
+---------------------
 
 The following example (`demo01 <https://github.com/pcubillos/MCcubed/blob/master/examples/demo01/demo01.py>`_) shows a basic MCMC run with ``MC3`` from
 the Python interpreter.
@@ -190,20 +190,17 @@ The plots sub-package provides the plotting functions:
 .. code-block:: python
 
    # Plot best-fitting model and binned data:
-   mc3.plots.modelfit(data, uncert, x, y, title="Best-fitting Model",
-                      savefile="quad_bestfit.png")
+   mc3.plots.modelfit(data, uncert, x, y, savefile="quad_bestfit.png")
    # Plot trace plot:
    parname = ["constant", "linear", "quadratic"]
-   mc3.plots.trace(posterior, Zchain, title="Fitting-parameter Trace Plots",
-                   parname=parname, savefile="quad_trace.png")
+   mc3.plots.trace(posterior, Zchain, parname=parname, savefile="quad_trace.png")
 
    # Plot pairwise posteriors:
-   mc3.plots.pairwise(posterior, title="Pairwise posteriors", parname=parname,
-                      savefile="quad_pairwise.png")
+   mc3.plots.pairwise(posterior, parname=parname, savefile="quad_pairwise.png")
 
-   # Plot marginal posterior histograms:
-   mc3.plots.histogram(posterior, title="Marginal posterior histograms",
-                 parname=parname, savefile="quad_hist.png", percentile=0.683)
+   # Plot marginal posterior histograms (with 68% highest-posterior-density credible regions):
+   mc3.plots.histogram(posterior, parname=parname, savefile="quad_hist.png",
+                       percentile=0.683)
 
 .. image:: ./quad_bestfit.png
    :width: 50%
@@ -222,8 +219,8 @@ The plots sub-package provides the plotting functions:
           MCMC run (see `File Outputs
           <http://pcubillos.github.io/MCcubed/tutorial.html#file-outputs>`_).
 
-Example 2 (Shell Run)
----------------------
+Example 2: Shell Run
+--------------------
 
 The following example
 (`demo02 <https://github.com/pcubillos/MCcubed/blob/master/examples/demo02/>`_)

docs/index.rst

@@ -56,6 +56,7 @@ development of this package.
 - Nate Lust (UCF)
 - `AJ Foster <http://aj-foster.com>`_ (UCF)
 - Madison Stemm (UCF)
+- Tom Loredo (Cornell)
 - Kevin Stevenson (UCF)
 - Chris Campo (UCF)
 - Matt Hardin (UCF)
@@ -70,6 +71,7 @@ Documentation
    getstarted
    mctutorial
    fittutorial
+   timeaveraging
    contributing
    license
 
@@ -86,8 +88,8 @@ Be Kind
 
 
 Please cite this paper if you found ``MC3`` useful for your research:
-  `Cubillos et al. 2016: On the Correlated Noise Analyses Applied to
-  Exoplanet Light Curves`_, submitted.
+  `Cubillos et al. 2016: On the Correlated-noise Analyses Applied to
+  Exoplanet Light Curves <https://arxiv.org/abs/1610.01336>`_, ApJ.
 
 We welcome your feedback, but do not necessarily guarantee support.
 Please send feedback or inquiries to:

docs/mctutorial.rst

@@ -52,16 +52,16 @@ The following code block shows an example for an MC3 configuration file:
   # The data and uncertainties:
   data      = data.npz
 
-MCMC-run Configuration
-----------------------
+MCMC Run
+--------
 
 This example describes the basic MCMC argument configuration.
 The following sub-sections make up a script meant to be run from the Python
 interpreter.  The complete example script is located at `tutorial01 <https://github.com/pcubillos/MCcubed/blob/master/examples/tutorial01/tutorial01.py>`_.
 
 
-Data and Data Uncertainties
-^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Data
+^^^^
 
 The ``data`` argument (required) defines the dataset to be fitted.
 This argument can be either a 1D float ndarray or the filename (a string)
@@ -307,8 +307,8 @@ argument:
 
 .. _mcchains:
 
-MCMC Chains Configuration
-^^^^^^^^^^^^^^^^^^^^^^^^^
+MCMC Config
+^^^^^^^^^^^
 
 The following arguments set the MCMC chains configuration:
 
@@ -384,8 +384,8 @@ uncertainties by a common scale factor.
    chisqscale = False  # Scale the data uncertainties such that red. chisq = 1
 
 
-Gelman-Rubin Convergence Test
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Convergence Test
+^^^^^^^^^^^^^^^^
 
 The ``grtest`` argument (optional, boolean, default=False) is a flag that
 indicates MC3 to run the Gelman-Rubin convergence test for the MCMC sample of
@@ -513,11 +513,13 @@ When run from a pyhton interactive session, ``MC3`` will return six arrays:
   include the values for all model parameters, including fixed and
   shared parameters, whereas ``posterior`` includes only
   the free parameters.  Be careful with the dimesions.
+..
+
+   Resume a previous MC3 Run
+   ^^^^^^^^^^^^^^^^^^^^^^^^^
 
-Resume a previous MC3 Run
-^^^^^^^^^^^^^^^^^^^^^^^^^
+   TBD
 
-TBD
 
 Inputs from Files
 -----------------

docs/timeaveraging.rst

@@ -0,0 +1,70 @@
+.. _timeaveraging:
+
+Time Averaging
+==============
+
+The ``MCcubed.rednoise.binrms`` routine computes the binned RMS array
+(as function of bin size) used in the the time-averaging procedure.
+Given a (model-data) residuals array.  The routine returns the RMS of
+the binend data (:math:`{\rm rms}_N`), the lower and upper RMS
+uncertainties, the extrapolated RMS for Gaussian (white) noise
+(:math:`\sigma_N`), and the bin-size array (:math:`N`).
+
+This function uses an asymptotic approximation to compute the RMS
+uncertainties (:math:`\sigma_{\rm rms} = \sqrt{{\rm rms}_N / 2M}`) for
+number of bins :math:`M> 35`.  For smaller values of :math:`M`
+(equivalently, large bin size) this routine computes the errors from
+the posterior PDF of the RMS (an inverse-gamma distribution).  See
+[Cubillos2016]_.
+
+
+
+Example
+^^^^^^^
+
+.. code-block:: python
+
+  import numpy as np
+  import matplotlib.pyplot as plt
+  import MCcubed as mc3   # Add path to mc3 if necessary
+  plt.ion()
+
+  # Generate residuals signal:
+  N = 1000
+  # White-noise signal:
+  white = np.random.normal(0, 5, N)
+  # (Sinusoidal) time-correlated signal:
+  red = np.sin(np.arange(N)/(0.1*N))*np.random.normal(1.0, 1.0, N)
+
+  # Plot the time-correlated residuals signal:
+  plt.figure(0)
+  plt.clf()
+  plt.plot(white+red, ".k")
+  plt.ylabel("Residuals", fontsize=14)
+
+  # Compute the residuals rms-vs-binsize:
+  maxbins = N/5
+  rms, rmslo, rmshi, stderr, binsz = mc3.rednoise.binrms(white+red, maxbins)
+
+  # Plot the rms with error bars along with the Gaussian standard deviation curve:
+  plt.figure(-6)
+  plt.clf()
+  plt.errorbar(binsz, rms, yerr=[rmslo, rmshi], fmt="k-", ecolor='0.5', capsize=0, label="Data RMS")
+  plt.loglog(binsz, stderr, color='red', ls='-', lw=2, label="Gaussian std.")
+  plt.xlim(1,200)
+  plt.legend(loc="upper right")
+  plt.xlabel("Bin size", fontsize=14)
+  plt.ylabel("RMS", fontsize=14)
+
+.. image:: ./time-correlated_signal.png
+   :width: 50%
+
+.. image:: ./rms-vs-binsize.png
+   :width: 50%
+
+References
+^^^^^^^^^^
+
+.. [Cubillos2016] `Cubillos et al. (2016): On the Correlated-noise
+  Analyses Applied to Exoplanet Light Curves
+  <https://arxiv.org/abs/1610.01336>`_