add (back) student t likelihood (#61); new example to demonstrate it

examples/51Peg/kima_setup.cpp

@@ -8,6 +8,7 @@ const bool trend = false;
 const int degree = 0;
 const bool multi_instrument = false;
 const bool known_object = false;
+const bool studentt = false;
 
 RVmodel::RVmodel():fix(false),npmax(1)
 {

examples/API/kima_setup.cpp

@@ -31,6 +31,9 @@ const bool multi_instrument = false;
 /// include a (better) known extra Keplerian curve? (KO mode!)
 const bool known_object = false;
 
+/// use a Student-t distribution for the likelihood (instead of Gaussian)
+const bool studentt = false;
+
 /// whether the model includes hyper-priors for the orbital period and semi-amplitude
 const bool hyperpriors = false;
 

examples/BL2009/kima_setup.cpp

@@ -8,6 +8,7 @@ const bool trend = false;
 const int degree = 0;
 const bool multi_instrument = false;
 const bool known_object = false;
+const bool studentt = false;
 
 RVmodel::RVmodel():fix(false),npmax(1)
 {

examples/CoRoT7/kima_setup.cpp

@@ -8,6 +8,7 @@ const bool trend = false;
 const int degree = 0;
 const bool multi_instrument = false;
 const bool known_object = false;
+const bool studentt = false;
 
 RVmodel::RVmodel():fix(false),npmax(5)
 {

examples/activity_correlations/kima_setup.cpp

@@ -8,6 +8,7 @@ const bool trend = false;
 const int degree = 0;
 const bool multi_instrument = false;
 const bool known_object = false;
+const bool studentt = false;
 
 RVmodel::RVmodel():fix(false),npmax(2)
 {

examples/default_priors/kima_setup.cpp

@@ -8,6 +8,7 @@ const bool trend = false;
 const int degree = 0;
 const bool multi_instrument = false;
 const bool known_object = false;
+const bool studentt = false;
 
 RVmodel::RVmodel():fix(true),npmax(1)
 {}

examples/many_planets/kima_setup.cpp

@@ -8,6 +8,7 @@ const bool trend = false;
 const int degree = 0;
 const bool multi_instrument = false;
 const bool known_object = false;
+const bool studentt = false;
 
 RVmodel::RVmodel():fix(false),npmax(10)
 {

examples/multi_instrument/kima_setup.cpp

@@ -8,6 +8,7 @@ const bool trend = false;
 const int degree = 0;
 const bool multi_instrument = true;
 const bool known_object = false;
+const bool studentt = false;
 
 RVmodel::RVmodel():fix(false),npmax(1)
 {

examples/studentT/Makefile

@@ -0,0 +1,2 @@
+KIMA_DIR = ../..
+include $(KIMA_DIR)/examples.mk

examples/studentT/OPTIONS

@@ -0,0 +1,13 @@
+# File containing parameters for DNest4
+# Put comments at the top, or at the end of the line.
+2	# Number of particles
+1000	# new level interval
+500	# save interval
+10	# threadSteps: number of steps each thread does independently before communication
+0	# maximum number of levels
+10	# Backtracking scale length (lambda)
+100	# Strength of effect to force histogram to equal push (beta)
+5000	# Maximum number of saves (0 = infinite)
+    # (optional) samples file
+    # (optional) sample_info file
+    # (optional) levels file

examples/studentT/README.md

@@ -0,0 +1,29 @@
+This example showcases the use of a Student t likelihood for more robustness
+regarding outlier RV measurements. The folder contains a simulated dataset with
+known outliers and a planetary signal. 
+
+The `kima_setup.cpp` file sets the main options for the model, including
+
+```c++
+const bool studentt = true;
+```
+
+This introduces an additional parameter in the model, the degrees of freedom of
+the Student t likelihood (called `nu` in the code).
+
+The number of planets in the model is fixed to 1. This example uses default
+priors for all parameters. The default prior for `nu` is a log-uniform
+distribution between 2 and 1000. 
+Note that the variance of the t-distribution is infinite for 1 < `nu` <= 2.
+Larger values for `nu` make the t-distribution approach a Gaussian distribution.
+
+To compile and run, type
+
+```
+kima-run
+```
+
+With the default options, this example should take about a minute to finish.
+
+It's interesting to set `studentt = false` in order to check the nasty effect of
+the outliers on the Gaussian distribution.

examples/studentT/create_data.py

@@ -0,0 +1,59 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.stats import t as Tstudent, norm
+
+from pykima.keplerian import keplerian
+
+# reproducible
+np.random.seed(43)
+
+vsys = 11.1  # (m/s)
+P = 32.1  # (days)
+K = 4.6  # (m/s)
+ecc = 0.1
+w = 0.3
+
+
+def create_d1():
+    N = 68
+    # from 01/03/2010 to 24/01/2019, because why not
+    t = np.sort(np.random.uniform(55256, 58507, N))
+
+    rv = np.full_like(t, vsys)
+
+    Tp = t[0] + 10
+    planet = keplerian(t, P, K, ecc, w, Tp, 0.0)
+    rv += planet
+
+    srv = np.random.uniform(0.6, 0.9, t.size)
+    jit = 0.3
+
+    err = norm(0, np.hypot(np.median(srv), jit)).rvs(t.size)
+    # err = Tstudent(df=2.1, loc=0, scale=np.hypot(np.median(srv), jit)).rvs(t.size)
+    rv += err
+
+    # outliers!!!
+    rv[10] += 16.1
+    rv[15] += 20
+    rv[56] -= 13
+
+    header = 'bjd\tvrad\tsvrad\n---\t----\t-----'
+    np.savetxt('d1.txt', np.c_[t, rv, srv],
+               header=header, comments='', fmt='%6.3f')
+
+    return t, rv, srv, (P, K, ecc, w, Tp, vsys, jit)
+
+
+def plot_dataset():
+    t, rv, srv, pars = create_d1()
+
+    tt = np.linspace(t[0], t[-1], 5000)
+    pl = keplerian(tt, *pars[:6])
+
+    _, ax = plt.subplots(1, 1)
+    ax.set(xlabel='Time [days]', ylabel='RV [m/s]')
+
+    ax.errorbar(t, rv, srv, fmt='o')
+    ax.plot(tt, pl, 'k', zorder=-1)
+
+    plt.show()

examples/studentT/kima_setup.cpp

@@ -0,0 +1,31 @@
+#include "kima.h"
+
+const bool obs_after_HARPS_fibers = false;
+const bool GP = false;
+const bool MA = false;
+const bool hyperpriors = false;
+const bool trend = false;
+const int degree = 0;
+const bool multi_instrument = false;
+const bool known_object = false;
+const bool studentt = true;  // use a Student's t-distribution for the likelihood
+
+RVmodel::RVmodel():fix(true),npmax(1)
+{
+    // all default priors
+    // the default prior for the degrees of freedom (nu) of the t-distribution:
+    // nu_prior = make_prior<LogUniform>(2, 1000);
+}
+
+
+int main(int argc, char** argv)
+{
+    datafile = "d1.txt";
+
+    load(datafile, "ms", 2);
+
+    Sampler<RVmodel> sampler = setup<RVmodel>(argc, argv);
+    sampler.run(50);
+
+    return 0;
+}

src/RVmodel.cpp

@@ -8,6 +8,7 @@
 #include <limits>
 #include <fstream>
 #include <chrono>
+#include <cmath>
 #include <time.h>
 
 using namespace std;
@@ -76,6 +77,9 @@ void RVmodel::setPriors()  // BUG: should be done by only one thread!
             throw std::logic_error("When known_object=true, please set all priors: KO_Pprior, KO_Kprior, KO_eprior, KO_phiprior, KO_wprior");
     }
 
+    if (studentt)
+        nu_prior = make_prior<LogUniform>(2, 1000);
+
 }
 
 
@@ -146,6 +150,10 @@ void RVmodel::from_prior(RNG& rng)
         KO_w = KO_wprior->generate(rng);
     }
 
+    if (studentt)
+        nu = nu_prior->generate(rng);
+
+
     calculate_mu();
 
     if(GP) calculate_C();
@@ -707,6 +715,10 @@ double RVmodel::perturb(RNG& rng)
                 Jprior->perturb(extra_sigma, rng);
             }
 
+            if (studentt)
+                nu_prior->perturb(nu, rng);
+
+
             if (known_object)
             {
                 remove_known_object();
@@ -876,31 +888,44 @@ double RVmodel::log_likelihood() const
     }
     else
     {
-        // The following code calculates the log likelihood
-        // in the case of a t-Student model
-        //  for(size_t i=0; i<y.size(); i++)
-        //  {
-        //      var = sig[i]*sig[i] + extra_sigma*extra_sigma;
-        //      logL += gsl_sf_lngamma(0.5*(nu + 1.)) - gsl_sf_lngamma(0.5*nu)
-        //          - 0.5*log(M_PI*nu) - 0.5*log(var)
-        //          - 0.5*(nu + 1.)*log(1. + pow(y[i] - mu[i], 2)/var/nu);
-        //  }
-
-        // The following code calculates the log likelihood
-        // in the case of a Gaussian likelihood
-        double var, jit;
-        for(size_t i=0; i<N; i++)
-        {
-            if(multi_instrument)
+        if (studentt){
+            // The following code calculates the log likelihood 
+            // in the case of a t-Student model
+            double var, jit;
+            for(size_t i=0; i<N; i++)
             {
-                jit = jitters[obsi[i]-1];
-                var = sig[i]*sig[i] + jit*jit;
+                if(multi_instrument)
+                {
+                    jit = jitters[obsi[i]-1];
+                    var = sig[i]*sig[i] + jit*jit;
+                }
+                else
+                    var = sig[i]*sig[i] + extra_sigma*extra_sigma;
+
+                logL += std::lgamma(0.5*(nu + 1.)) - std::lgamma(0.5*nu)
+                        - 0.5*log(M_PI*nu) - 0.5*log(var)
+                        - 0.5*(nu + 1.)*log(1. + pow(y[i] - mu[i], 2)/var/nu);
             }
-            else
-                var = sig[i]*sig[i] + extra_sigma*extra_sigma;
 
-            logL += - halflog2pi - 0.5*log(var)
-                    - 0.5*(pow(y[i] - mu[i], 2)/var);
+        }
+
+        else{
+            // The following code calculates the log likelihood
+            // in the case of a Gaussian likelihood
+            double var, jit;
+            for(size_t i=0; i<N; i++)
+            {
+                if(multi_instrument)
+                {
+                    jit = jitters[obsi[i]-1];
+                    var = sig[i]*sig[i] + jit*jit;
+                }
+                else
+                    var = sig[i]*sig[i] + extra_sigma*extra_sigma;
+
+                logL += - halflog2pi - 0.5*log(var)
+                        - 0.5*(pow(y[i] - mu[i], 2)/var);
+            }
         }
 
     }
@@ -959,22 +984,25 @@ void RVmodel::print(std::ostream& out) const
     if(GP)
     {
         if (kernel == standard)
-            out<<eta1<<'\t'<<eta2<<'\t'<<eta3<<'\t'<<eta4<<'\t';
+            out << eta1 << '\t' << eta2 << '\t' << eta3 << '\t' << eta4 << '\t';
         else
-            out<<eta1<<'\t'<<eta2<<'\t'<<eta3<<'\t';
+            out << eta1 << '\t' << eta2 << '\t' << eta3 << '\t';
     }
 
     if(MA)
-        out<<sigmaMA<<'\t'<<tauMA<<'\t';
+        out << sigmaMA << '\t' << tauMA << '\t';
 
     if(known_object) // KO mode!
         out << KO_P << "\t" << KO_K << "\t" << KO_phi << "\t" << KO_e << "\t" << KO_w << "\t";
 
-
     planets.print(out);
 
-    out<<' '<<staleness<<' ';
-    out<<background;
+    out << ' ' << staleness << ' ';
+
+    if (studentt)
+        out << '\t' << nu << '\t';
+
+    out << background;
 }
 
 string RVmodel::description() const
@@ -1036,7 +1064,11 @@ string RVmodel::description() const
     if (planets.get_max_num_components()>0)
         desc += "P   K   phi   ecc   w   ";
 
-    desc += "staleness   vsys";
+    desc += "staleness   ";
+    if (studentt)
+        desc += "nu   ";
+
+    desc += "vsys";
 
     return desc;
 }
@@ -1065,6 +1097,7 @@ void RVmodel::save_setup() {
     fout << "degree: " << degree << endl;
     fout << "multi_instrument: " << multi_instrument << endl;
     fout << "known_object: " << known_object << endl;
+    fout << "studentt: " << studentt << endl;
     fout << "indicator_correlations: " << data.indicator_correlations << endl;
     fout << "indicators: ";
     for (auto f: data.indicator_names){
@@ -1098,6 +1131,8 @@ void RVmodel::save_setup() {
         fout << "fiber_offset_prior: " << *fiber_offset_prior << endl;
     if (multi_instrument)
         fout << "offsets_prior: " << *offsets_prior << endl;
+    if (studentt)
+        fout << "nu_prior: " << *nu_prior << endl;
 
     if (GP){
         fout << endl << "[priors.GP]" << endl;