In [1]:

    

from __future__ import division

import pandas as pd
from exosyspop.populations import BinaryPopulation
from exosyspop.populations import TRILEGAL_BGBinaryPopulation
from exosyspop.populations import KeplerBinaryPopulation, PoissonPlanetPopulation
from exosyspop.populations import KeplerPowerLawBinaryPopulation
from exosyspop.populations import PopulationMixture

targets = pd.read_hdf('targets.h5')
bgstars = pd.read_hdf('bgstars.h5')

# Sanitize dtypes of targets DataFrame
for c in targets.columns:
    if targets[c].dtype == object:
        targets.loc[:,c] = targets.loc[:,c].astype(str)


import logging
rootLogger = logging.getLogger()


    

    




/u/tdm/anaconda/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')
WARNING:root:progressbar not imported

In [1]:

    

pop = KeplerPowerLawBinaryPopulation.load('plaw_pop')
pop.set_params(period_min=20, period_max=1200, beta=-0.95, fB=0.14)
catalog = pop.observe(new=True, regr_trap=True)


    

    




---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-1-465a672ac200> in <module>()
----> 1 pop = KeplerPowerLawBinaryPopulation.load('plaw_pop')
      2 pop.set_params(period_min=20, period_max=1200, beta=-0.95, fB=0.14)
      3 catalog = pop.observe(new=True, regr_trap=True)

NameError: name 'KeplerPowerLawBinaryPopulation' is not defined

In [3]:

    

from exosyspop.survey import DetectionThreshold, DetectionRamp

eff = DetectionRamp(6,16)


    

In [24]:

    

import sys
sys.path.append('..')

from simpleabc.simple_abc import Model, basic_abc, pmc_abc
from scipy.stats import gaussian_kde, entropy, anderson_ksamp, uniform
import numpy as np

class PopulationModel(Model):
    """
    Test model for stellar binary population where parameters are fB, beta
    """
    def __init__(self, poplist, eff=None):
        self.poplist = poplist     
        self.eff = eff
        self.period_min = poplist.params['period_min']
        self.period_max = poplist.params['period_max']
        
        self.distance_norms = None
        
    bounds = [(0,1), (-1.5,0), (0,5),(0,5)]
    prior = [uniform(0,1), uniform(-1.5, 1.5), uniform(0,5), uniform(0,5)]
        
    def draw_theta(self):
        """ Draw parameters from prior
        """
        return [p.rvs() for p in self.prior]
        
    def generate_data(self, theta):
        """Generates synthetic catalog
        """
        fB, beta, beta_a, beta_b = theta
        self.poplist.set_params(fB=fB, beta=beta, beta_a=beta_a, beta_b=beta_b)
        try:
            return self.poplist.observe(new=True, 
                                        regr_trap=True).observe(self.eff)
        except:
            print('Error!  theta={}'.format(theta))
            
    def summary_stats(self, data):
        """Computes summary statistics from data
        """
        if data is None:
            return [np.nan]*3
        
        N = len(data)
        
        try:
            Pmin, Pmax = np.log(self.period_min), np.log(self.period_max)
            Pgrid = np.linspace(Pmin, Pmax, 1000)
            if N > 1:
                k = gaussian_kde(np.log(data.period.values))
                p = k(Pgrid)
            else:
                p = np.ones(len(Pgrid))*1./(Pmax - Pmin)
        except ValueError:
            print(data.period.values)
            raise
        
        phase_sec = data.phase_sec.dropna().values
        
        return p, N, phase_sec
        
    def d_period(self, summary_stats, summary_stats_synth):
        p1, _, _ = summary_stats
        p2, _, _ = summary_stats_synth
        try:
            len(p1)
            len(p2)
        except:
            return np.inf
        kl_period = entropy(p1, p2)
        return kl_period
    
    def Ndist(self, N1, N2):
        if N1==0. or N2==0.:
            dist = 1
        else:
            dist = max(1 - 1.*N1/N2, 1-1*N2/N1)
        return dist
            
    def d_N(self, summary_stats, summary_stats_synth):
        _, N1, _ = summary_stats
        _, N2, _ = summary_stats_synth
        return self.Ndist(N1, N2) 
        
    def d_fsec(self, summary_stats, summary_stats_synth):
        _, N1, phase_sec1 = summary_stats
        _, N2, phase_sec2 = summary_stats_synth
        
        f_sec1 = len(phase_sec1)/float(N1)
        f_sec2 = len(phase_sec2)/float(N2)
        
        return np.absolute(f_sec1 - f_sec2)

    def d_phase(self, summary_stats, summary_stats_synth, nbins=11):
        _, _, phase_sec1 = summary_stats
        _, _, phase_sec2 = summary_stats_synth

        try:
            len(phase_sec2)
        except:
            return np.inf
        
        k1 = gaussian_kde(phase_sec1)
        k2 = gaussian_kde(phase_sec2)

        phs = np.linspace(0,1,100)
        pdf1 = k1(phs)
        pdf2 = k2(phs)
        
        return entropy(pdf1, pdf2)
    
    distance_functions = ('d_period', 'd_N', 'd_fsec', 'd_phase')

    def null_distance_test(self, theta_0, N=100):
        data1 = [self.generate_data(theta_0) for i in range(N)]
        data2 = [self.generate_data(theta_0) for i in range(N)]
        
        ds = []
        for dfn in self.distance_functions:
            fn = getattr(self, dfn)
            ds.append([fn(self.summary_stats(data1[i]),
                         self.summary_stats(data2[i])) for i in range(N)])
        
        self._null_stds = np.array([np.std(d) for d in ds])
        self.distance_norms = self._null_stds / self._null_stds[0]
        
    def distance_function(self, summary_stats, summary_stats_synth):
        """Computes distance
        """
        ds = []
        for dfn in self.distance_functions:
            fn = getattr(self, dfn)
            ds.append(fn(summary_stats, summary_stats_synth))
        
        return np.sum([d / self.distance_norms[i] for i,d in enumerate(ds)])
            
        #return d1 + d2 * (0.015/0.072)  #renormalized based on null test
        #return (d1 + d2 + d3 + d4)/4.


    

In [25]:

    

model = PopulationModel(PopulationMixture([pop]))
theta_0 = 0.14, -0.95, 0.8, 2.0
data = model.generate_data(theta_0)
model.set_data(data)


    

In [26]:

    

model.null_distance_test(theta_0, N=100)


    

In [27]:

    

model.distance_norms


    

    
Out[27]:





array([ 1.        ,  3.14920833,  1.97371597,  1.99408711])

In [8]:

    

posterior = basic_abc(model, data, min_samples=100, epsilon=0.4, verbose=True)


    

    




10 samples accepted (10.0%).
24 samples accepted (12.0%).
35 samples accepted (11.7%).
41 samples accepted (10.2%).
53 samples accepted (10.6%).
67 samples accepted (11.2%).
Error!  fB=0.000101185769932, beta=-0.810974815254
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-4-aa75c292b951>", line 127, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-4-aa75c292b951>", line 66, in d_period
    kl_period = entropy(p1, p2)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/_distn_infrastructure.py", line 2498, in entropy
    if len(qk) != len(pk):
TypeError: len() of unsized object

81 samples accepted (11.6%).
92 samples accepted (11.5%).
Error!  fB=0.00359319682528, beta=-0.332810440468
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-4-aa75c292b951>", line 127, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-4-aa75c292b951>", line 66, in d_period
    kl_period = entropy(p1, p2)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/_distn_infrastructure.py", line 2498, in entropy
    if len(qk) != len(pk):
TypeError: len() of unsized object

In [28]:

    

pmc_posterior = pmc_abc(model, data, epsilon_0=0.3, min_samples=200, steps=20, verbose=True)


    

    




Starting step 0, epsilon=0.3
Error!  theta=[0.0027675345111886074, -0.04173890078481546, 0.08395916653485547, 4.353831790534039]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

10 samples accepted (10.0%).
15 samples accepted (7.5%).
22 samples accepted (7.3%).
31 samples accepted (7.8%).
35 samples accepted (7.0%).
44 samples accepted (7.3%).
Error!  theta=[0.005375463660796331, -0.5558756286084341, 2.1497139562112384, 3.8248304953532153]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

46 samples accepted (6.6%).
56 samples accepted (7.0%).
Error!  theta=[0.0022767346864469573, -0.6678270906681323, 3.277281100756543, 1.700603240582625]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

60 samples accepted (6.7%).
Error!  theta=[0.0015212031886308042, -1.4938851339439052, 2.118383016702849, 3.5906205226548273]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

65 samples accepted (6.5%).
72 samples accepted (6.5%).
82 samples accepted (6.8%).
Error!  theta=[0.000746935814917693, -0.39647852000587624, 0.7969539680702425, 1.3432517993791981]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

93 samples accepted (7.2%).
99 samples accepted (7.1%).
106 samples accepted (7.1%).
113 samples accepted (7.1%).
121 samples accepted (7.1%).
129 samples accepted (7.2%).
139 samples accepted (7.3%).
Error!  theta=[0.0016343664595199758, -1.4676699027887148, 1.094270225076, 1.4630669285886384]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

Error!  theta=[4.7489349074059106e-05, -0.6973740371735735, 4.930244133509744, 0.7572229060287711]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

146 samples accepted (7.3%).
151 samples accepted (7.2%).
160 samples accepted (7.3%).
173 samples accepted (7.5%).
178 samples accepted (7.4%).
Error!  theta=[0.00010512316139965616, -0.29189943278180897, 0.47236264345677126, 0.6523785136140364]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

191 samples accepted (7.6%).
199 samples accepted (7.7%).
Starting step 1, epsilon=0.277278499196
15 samples accepted (15.0%).
36 samples accepted (18.0%).
56 samples accepted (18.7%).
75 samples accepted (18.8%).
97 samples accepted (19.4%).
118 samples accepted (19.7%).
139 samples accepted (19.9%).
162 samples accepted (20.2%).
189 samples accepted (21.0%).
Starting step 2, epsilon=0.252217696824
25 samples accepted (25.0%).
54 samples accepted (27.0%).
85 samples accepted (28.3%).
102 samples accepted (25.5%).
Error!  theta=[ 0.00354989 -0.36527296  3.50009842  2.50991425]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

125 samples accepted (25.0%).
150 samples accepted (25.0%).
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 103, in d_phase
    k2 = gaussian_kde(phase_sec2)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 171, in __init__
    self.set_bandwidth(bw_method=bw_method)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 498, in set_bandwidth
    self._compute_covariance()
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 509, in _compute_covariance
    self._data_inv_cov = linalg.inv(self._data_covariance)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/linalg/basic.py", line 687, in inv
    raise LinAlgError("singular matrix")
LinAlgError: singular matrix

171 samples accepted (24.4%).
191 samples accepted (23.9%).
Starting step 3, epsilon=0.228895782674
24 samples accepted (24.0%).
48 samples accepted (24.0%).
70 samples accepted (23.3%).
Error!  theta=[  2.83210896e-03  -1.31448032e+00   1.00028494e+00   2.99016473e+00]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

87 samples accepted (21.8%).
108 samples accepted (21.6%).
130 samples accepted (21.7%).
156 samples accepted (22.3%).
Error!  theta=[  9.03585135e-05  -1.26792070e+00   8.42405242e-01   3.55590681e+00]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

178 samples accepted (22.2%).
Starting step 4, epsilon=0.204313660948
27 samples accepted (27.0%).
50 samples accepted (25.0%).
75 samples accepted (25.0%).
91 samples accepted (22.8%).
Error!  theta=[  4.52890052e-04  -1.33958680e+00   1.20729724e+00   9.30716535e-01]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

120 samples accepted (24.0%).
150 samples accepted (25.0%).
176 samples accepted (25.1%).
Starting step 5, epsilon=0.181857606856
Error!  theta=[ 0.00184027 -1.26971165  1.0193786   1.40290299]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

13 samples accepted (13.0%).
31 samples accepted (15.5%).
46 samples accepted (15.3%).
62 samples accepted (15.5%).
73 samples accepted (14.6%).
80 samples accepted (13.3%).
Error!  theta=[ 0.00436981 -0.81424867  1.67751932  2.45230749]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

92 samples accepted (13.1%).
104 samples accepted (13.0%).
120 samples accepted (13.3%).
Error!  theta=[ 0.00346171 -1.28696013  3.46136271  2.70422445]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

Error!  theta=[  6.86972947e-04  -1.43261058e+00   4.50221235e+00   4.97374802e+00]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

133 samples accepted (13.3%).
144 samples accepted (13.1%).
159 samples accepted (13.2%).
177 samples accepted (13.6%).
197 samples accepted (14.1%).
Starting step 6, epsilon=0.159525886702
12 samples accepted (12.0%).
29 samples accepted (14.5%).
47 samples accepted (15.7%).
Error!  theta=[  6.84766303e-04  -4.95944342e-01   3.67207372e+00   4.40259114e+00]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

64 samples accepted (16.0%).
76 samples accepted (15.2%).
94 samples accepted (15.7%).
101 samples accepted (14.4%).
113 samples accepted (14.1%).
130 samples accepted (14.4%).
142 samples accepted (14.2%).
157 samples accepted (14.3%).
170 samples accepted (14.2%).
187 samples accepted (14.4%).
Starting step 7, epsilon=0.143145748714
19 samples accepted (19.0%).
30 samples accepted (15.0%).
53 samples accepted (17.7%).
63 samples accepted (15.8%).
76 samples accepted (15.2%).
96 samples accepted (16.0%).
111 samples accepted (15.9%).
124 samples accepted (15.5%).
137 samples accepted (15.2%).
161 samples accepted (16.1%).
177 samples accepted (16.1%).
191 samples accepted (15.9%).
Starting step 8, epsilon=0.129159888604
10 samples accepted (10.0%).
27 samples accepted (13.5%).
44 samples accepted (14.7%).
56 samples accepted (14.0%).
67 samples accepted (13.4%).
76 samples accepted (12.7%).
89 samples accepted (12.7%).
106 samples accepted (13.2%).
122 samples accepted (13.6%).
136 samples accepted (13.6%).
150 samples accepted (13.6%).
163 samples accepted (13.6%).
179 samples accepted (13.8%).
194 samples accepted (13.9%).
Starting step 9, epsilon=0.118836005755
10 samples accepted (10.0%).
20 samples accepted (10.0%).
38 samples accepted (12.7%).
48 samples accepted (12.0%).
56 samples accepted (11.2%).
62 samples accepted (10.3%).
79 samples accepted (11.3%).
96 samples accepted (12.0%).
110 samples accepted (12.2%).
124 samples accepted (12.4%).
141 samples accepted (12.8%).
158 samples accepted (13.2%).
167 samples accepted (12.8%).
179 samples accepted (12.8%).
193 samples accepted (12.9%).
199 samples accepted (12.4%).
Starting step 10, epsilon=0.106017549056
9 samples accepted (9.0%).
19 samples accepted (9.5%).
33 samples accepted (11.0%).
41 samples accepted (10.2%).
54 samples accepted (10.8%).
66 samples accepted (11.0%).
79 samples accepted (11.3%).
97 samples accepted (12.1%).
106 samples accepted (11.8%).
122 samples accepted (12.2%).
132 samples accepted (12.0%).
139 samples accepted (11.6%).
148 samples accepted (11.4%).
159 samples accepted (11.4%).
168 samples accepted (11.2%).
174 samples accepted (10.9%).
186 samples accepted (10.9%).
Starting step 11, epsilon=0.0968442618082
7 samples accepted (7.0%).
14 samples accepted (7.0%).
31 samples accepted (10.3%).
41 samples accepted (10.2%).
47 samples accepted (9.4%).
54 samples accepted (9.0%).
60 samples accepted (8.6%).
72 samples accepted (9.0%).
83 samples accepted (9.2%).
87 samples accepted (8.7%).
99 samples accepted (9.0%).
109 samples accepted (9.1%).
121 samples accepted (9.3%).
128 samples accepted (9.1%).
Error!  theta=[  3.68095582e-04  -1.25009819e+00   1.08970961e+00   6.79293593e-01]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

136 samples accepted (9.1%).
143 samples accepted (8.9%).
149 samples accepted (8.8%).
158 samples accepted (8.8%).
165 samples accepted (8.7%).
Error!  theta=[  2.18758086e-03  -1.33313626e+00   3.33271768e+00   3.99739384e+00]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

170 samples accepted (8.5%).
185 samples accepted (8.8%).
194 samples accepted (8.8%).
Starting step 12, epsilon=0.0900013828249
5 samples accepted (5.0%).
16 samples accepted (8.0%).
26 samples accepted (8.7%).
33 samples accepted (8.2%).
37 samples accepted (7.4%).
43 samples accepted (7.2%).
47 samples accepted (6.7%).
57 samples accepted (7.1%).
64 samples accepted (7.1%).
73 samples accepted (7.3%).
77 samples accepted (7.0%).
80 samples accepted (6.7%).
84 samples accepted (6.5%).
95 samples accepted (6.8%).
Error!  theta=[  6.52666441e-05  -1.23190628e+00   1.79010103e+00   3.01113622e+00]
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 91, in d_fsec
    f_sec2 = len(phase_sec2)/float(N2)
TypeError: object of type 'float' has no len()

107 samples accepted (7.1%).
117 samples accepted (7.3%).
120 samples accepted (7.1%).
125 samples accepted (6.9%).
132 samples accepted (6.9%).
139 samples accepted (7.0%).
144 samples accepted (6.9%).
151 samples accepted (6.9%).
157 samples accepted (6.8%).
166 samples accepted (6.9%).
175 samples accepted (7.0%).
180 samples accepted (6.9%).
188 samples accepted (7.0%).
193 samples accepted (6.9%).
Starting step 13, epsilon=0.0845961674655
6 samples accepted (6.0%).
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 103, in d_phase
    k2 = gaussian_kde(phase_sec2)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 171, in __init__
    self.set_bandwidth(bw_method=bw_method)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 498, in set_bandwidth
    self._compute_covariance()
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 509, in _compute_covariance
    self._data_inv_cov = linalg.inv(self._data_covariance)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/linalg/basic.py", line 687, in inv
    raise LinAlgError("singular matrix")
LinAlgError: singular matrix

16 samples accepted (8.0%).
26 samples accepted (8.7%).
33 samples accepted (8.2%).
39 samples accepted (7.8%).
45 samples accepted (7.5%).
55 samples accepted (7.9%).
65 samples accepted (8.1%).
71 samples accepted (7.9%).
74 samples accepted (7.4%).
81 samples accepted (7.4%).
85 samples accepted (7.1%).
93 samples accepted (7.2%).
98 samples accepted (7.0%).
108 samples accepted (7.2%).
112 samples accepted (7.0%).
119 samples accepted (7.0%).
124 samples accepted (6.9%).
133 samples accepted (7.0%).
143 samples accepted (7.1%).
150 samples accepted (7.1%).
157 samples accepted (7.1%).
161 samples accepted (7.0%).
166 samples accepted (6.9%).
173 samples accepted (6.9%).
185 samples accepted (7.1%).
193 samples accepted (7.1%).
Starting step 14, epsilon=0.0787459007417
7 samples accepted (7.0%).
13 samples accepted (6.5%).
19 samples accepted (6.3%).
26 samples accepted (6.5%).
32 samples accepted (6.4%).
34 samples accepted (5.7%).
39 samples accepted (5.6%).
43 samples accepted (5.4%).
51 samples accepted (5.7%).
63 samples accepted (6.3%).
71 samples accepted (6.5%).
76 samples accepted (6.3%).
80 samples accepted (6.2%).
89 samples accepted (6.4%).
94 samples accepted (6.3%).
101 samples accepted (6.3%).
108 samples accepted (6.4%).
111 samples accepted (6.2%).
119 samples accepted (6.3%).
127 samples accepted (6.3%).
136 samples accepted (6.5%).
141 samples accepted (6.4%).
146 samples accepted (6.3%).
156 samples accepted (6.5%).
163 samples accepted (6.5%).
172 samples accepted (6.6%).
179 samples accepted (6.6%).
183 samples accepted (6.5%).
190 samples accepted (6.6%).
196 samples accepted (6.5%).
Starting step 15, epsilon=0.0727617534787
7 samples accepted (7.0%).
10 samples accepted (5.0%).
16 samples accepted (5.3%).
23 samples accepted (5.8%).
29 samples accepted (5.8%).
35 samples accepted (5.8%).
42 samples accepted (6.0%).
48 samples accepted (6.0%).
51 samples accepted (5.7%).
57 samples accepted (5.7%).
62 samples accepted (5.6%).
66 samples accepted (5.5%).
70 samples accepted (5.4%).
78 samples accepted (5.6%).
84 samples accepted (5.6%).
91 samples accepted (5.7%).
93 samples accepted (5.5%).
98 samples accepted (5.4%).
102 samples accepted (5.4%).
108 samples accepted (5.4%).
114 samples accepted (5.4%).
118 samples accepted (5.4%).
126 samples accepted (5.5%).
129 samples accepted (5.4%).
132 samples accepted (5.3%).
140 samples accepted (5.4%).
146 samples accepted (5.4%).
149 samples accepted (5.3%).
150 samples accepted (5.2%).
152 samples accepted (5.1%).
157 samples accepted (5.1%).
160 samples accepted (5.0%).
163 samples accepted (4.9%).
171 samples accepted (5.0%).
173 samples accepted (4.9%).
179 samples accepted (5.0%).
182 samples accepted (4.9%).
186 samples accepted (4.9%).
192 samples accepted (4.9%).
199 samples accepted (5.0%).
Starting step 16, epsilon=0.0685772072093
6 samples accepted (6.0%).
10 samples accepted (5.0%).
18 samples accepted (6.0%).
25 samples accepted (6.2%).
29 samples accepted (5.8%).
35 samples accepted (5.8%).
41 samples accepted (5.9%).
51 samples accepted (6.4%).
53 samples accepted (5.9%).
56 samples accepted (5.6%).
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 103, in d_phase
    k2 = gaussian_kde(phase_sec2)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 171, in __init__
    self.set_bandwidth(bw_method=bw_method)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 498, in set_bandwidth
    self._compute_covariance()
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 509, in _compute_covariance
    self._data_inv_cov = linalg.inv(self._data_covariance)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/linalg/basic.py", line 687, in inv
    raise LinAlgError("singular matrix")
LinAlgError: singular matrix

60 samples accepted (5.5%).
64 samples accepted (5.3%).
67 samples accepted (5.2%).
70 samples accepted (5.0%).
73 samples accepted (4.9%).
75 samples accepted (4.7%).
83 samples accepted (4.9%).
86 samples accepted (4.8%).
92 samples accepted (4.8%).
95 samples accepted (4.8%).
99 samples accepted (4.7%).
99 samples accepted (4.5%).
103 samples accepted (4.5%).
108 samples accepted (4.5%).
112 samples accepted (4.5%).
120 samples accepted (4.6%).
124 samples accepted (4.6%).
127 samples accepted (4.5%).
132 samples accepted (4.6%).
135 samples accepted (4.5%).
140 samples accepted (4.5%).
143 samples accepted (4.5%).
147 samples accepted (4.5%).
150 samples accepted (4.4%).
156 samples accepted (4.5%).
159 samples accepted (4.4%).
165 samples accepted (4.5%).
170 samples accepted (4.5%).
176 samples accepted (4.5%).
185 samples accepted (4.6%).
187 samples accepted (4.6%).
193 samples accepted (4.6%).
196 samples accepted (4.6%).
199 samples accepted (4.5%).
Starting step 17, epsilon=0.06362508346
4 samples accepted (4.0%).
5 samples accepted (2.5%).
8 samples accepted (2.7%).
11 samples accepted (2.8%).
13 samples accepted (2.6%).
20 samples accepted (3.3%).
24 samples accepted (3.4%).
26 samples accepted (3.2%).
30 samples accepted (3.3%).
31 samples accepted (3.1%).
34 samples accepted (3.1%).
37 samples accepted (3.1%).
40 samples accepted (3.1%).
42 samples accepted (3.0%).
44 samples accepted (2.9%).
46 samples accepted (2.9%).
49 samples accepted (2.9%).
55 samples accepted (3.1%).
58 samples accepted (3.1%).
61 samples accepted (3.0%).
63 samples accepted (3.0%).
66 samples accepted (3.0%).
71 samples accepted (3.1%).
75 samples accepted (3.1%).
78 samples accepted (3.1%).
78 samples accepted (3.0%).
79 samples accepted (2.9%).
84 samples accepted (3.0%).
87 samples accepted (3.0%).
88 samples accepted (2.9%).
89 samples accepted (2.9%).
94 samples accepted (2.9%).
96 samples accepted (2.9%).
98 samples accepted (2.9%).
100 samples accepted (2.9%).
102 samples accepted (2.8%).
103 samples accepted (2.8%).
107 samples accepted (2.8%).
111 samples accepted (2.8%).
111 samples accepted (2.8%).
115 samples accepted (2.8%).
117 samples accepted (2.8%).
119 samples accepted (2.8%).
121 samples accepted (2.8%).
124 samples accepted (2.8%).
128 samples accepted (2.8%).
131 samples accepted (2.8%).
134 samples accepted (2.8%).
137 samples accepted (2.8%).
142 samples accepted (2.8%).
145 samples accepted (2.8%).
145 samples accepted (2.8%).
149 samples accepted (2.8%).
153 samples accepted (2.8%).
155 samples accepted (2.8%).
160 samples accepted (2.9%).
166 samples accepted (2.9%).
170 samples accepted (2.9%).
176 samples accepted (3.0%).
180 samples accepted (3.0%).
180 samples accepted (3.0%).
184 samples accepted (3.0%).
185 samples accepted (2.9%).
190 samples accepted (3.0%).
192 samples accepted (3.0%).
194 samples accepted (2.9%).
197 samples accepted (2.9%).
Starting step 18, epsilon=0.0602894730466
2 samples accepted (2.0%).
4 samples accepted (2.0%).
8 samples accepted (2.7%).
10 samples accepted (2.5%).
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 103, in d_phase
    k2 = gaussian_kde(phase_sec2)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 171, in __init__
    self.set_bandwidth(bw_method=bw_method)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 498, in set_bandwidth
    self._compute_covariance()
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 509, in _compute_covariance
    self._data_inv_cov = linalg.inv(self._data_covariance)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/linalg/basic.py", line 687, in inv
    raise LinAlgError("singular matrix")
LinAlgError: singular matrix

12 samples accepted (2.4%).
14 samples accepted (2.3%).
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 103, in d_phase
    k2 = gaussian_kde(phase_sec2)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 171, in __init__
    self.set_bandwidth(bw_method=bw_method)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 498, in set_bandwidth
    self._compute_covariance()
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 509, in _compute_covariance
    self._data_inv_cov = linalg.inv(self._data_covariance)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/linalg/basic.py", line 687, in inv
    raise LinAlgError("singular matrix")
LinAlgError: singular matrix

15 samples accepted (2.1%).
18 samples accepted (2.2%).
20 samples accepted (2.2%).
22 samples accepted (2.2%).
24 samples accepted (2.2%).
24 samples accepted (2.0%).
26 samples accepted (2.0%).
27 samples accepted (1.9%).
28 samples accepted (1.9%).
32 samples accepted (2.0%).
37 samples accepted (2.2%).
39 samples accepted (2.2%).
40 samples accepted (2.1%).
46 samples accepted (2.3%).
49 samples accepted (2.3%).
53 samples accepted (2.4%).
54 samples accepted (2.3%).
57 samples accepted (2.4%).
59 samples accepted (2.4%).
60 samples accepted (2.3%).
62 samples accepted (2.3%).
64 samples accepted (2.3%).
69 samples accepted (2.4%).
69 samples accepted (2.3%).
70 samples accepted (2.3%).
73 samples accepted (2.3%).
78 samples accepted (2.4%).
79 samples accepted (2.3%).
81 samples accepted (2.3%).
83 samples accepted (2.3%).
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 103, in d_phase
    k2 = gaussian_kde(phase_sec2)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 171, in __init__
    self.set_bandwidth(bw_method=bw_method)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 498, in set_bandwidth
    self._compute_covariance()
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 509, in _compute_covariance
    self._data_inv_cov = linalg.inv(self._data_covariance)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/linalg/basic.py", line 687, in inv
    raise LinAlgError("singular matrix")
LinAlgError: singular matrix

86 samples accepted (2.3%).
89 samples accepted (2.3%).
89 samples accepted (2.3%).
90 samples accepted (2.2%).
92 samples accepted (2.2%).
96 samples accepted (2.3%).
97 samples accepted (2.3%).
97 samples accepted (2.2%).
100 samples accepted (2.2%).
102 samples accepted (2.2%).
103 samples accepted (2.2%).
104 samples accepted (2.2%).
106 samples accepted (2.2%).
107 samples accepted (2.1%).
109 samples accepted (2.1%).
Error generating data!
Traceback (most recent call last):
  File "../simpleabc/simple_abc.py", line 245, in basic_abc
    synthetic_summary_stats)
  File "<ipython-input-24-d7325be88d5b>", line 132, in distance_function
    ds.append(fn(summary_stats, summary_stats_synth))
  File "<ipython-input-24-d7325be88d5b>", line 103, in d_phase
    k2 = gaussian_kde(phase_sec2)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 171, in __init__
    self.set_bandwidth(bw_method=bw_method)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 498, in set_bandwidth
    self._compute_covariance()
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/stats/kde.py", line 509, in _compute_covariance
    self._data_inv_cov = linalg.inv(self._data_covariance)
  File "/u/tdm/anaconda/lib/python2.7/site-packages/scipy/linalg/basic.py", line 687, in inv
    raise LinAlgError("singular matrix")
LinAlgError: singular matrix

111 samples accepted (2.1%).
111 samples accepted (2.1%).
112 samples accepted (2.1%).
116 samples accepted (2.1%).
118 samples accepted (2.1%).
119 samples accepted (2.1%).
122 samples accepted (2.1%).
124 samples accepted (2.1%).
126 samples accepted (2.1%).
130 samples accepted (2.1%).
135 samples accepted (2.2%).
137 samples accepted (2.2%).
140 samples accepted (2.2%).
141 samples accepted (2.2%).
145 samples accepted (2.2%).
147 samples accepted (2.2%).
149 samples accepted (2.2%).
152 samples accepted (2.2%).
156 samples accepted (2.2%).
157 samples accepted (2.2%).
158 samples accepted (2.2%).
159 samples accepted (2.2%).
164 samples accepted (2.2%).
167 samples accepted (2.2%).
169 samples accepted (2.2%).
172 samples accepted (2.2%).
175 samples accepted (2.2%).
178 samples accepted (2.3%).
178 samples accepted (2.2%).
180 samples accepted (2.2%).
185 samples accepted (2.3%).
188 samples accepted (2.3%).
189 samples accepted (2.2%).
194 samples accepted (2.3%).
198 samples accepted (2.3%).
Starting step 19, epsilon=0.0551133442272
3 samples accepted (3.0%).
6 samples accepted (3.0%).
8 samples accepted (2.7%).
9 samples accepted (2.2%).
14 samples accepted (2.8%).
20 samples accepted (3.3%).
23 samples accepted (3.3%).
24 samples accepted (3.0%).
26 samples accepted (2.9%).
29 samples accepted (2.9%).
33 samples accepted (3.0%).
36 samples accepted (3.0%).
38 samples accepted (2.9%).
40 samples accepted (2.9%).
42 samples accepted (2.8%).
44 samples accepted (2.8%).
45 samples accepted (2.6%).
49 samples accepted (2.7%).
52 samples accepted (2.7%).
56 samples accepted (2.8%).
59 samples accepted (2.8%).
60 samples accepted (2.7%).
62 samples accepted (2.7%).
63 samples accepted (2.6%).
64 samples accepted (2.6%).
68 samples accepted (2.6%).
70 samples accepted (2.6%).
72 samples accepted (2.6%).
73 samples accepted (2.5%).
73 samples accepted (2.4%).
75 samples accepted (2.4%).
79 samples accepted (2.5%).
84 samples accepted (2.5%).
88 samples accepted (2.6%).
90 samples accepted (2.6%).
92 samples accepted (2.6%).
98 samples accepted (2.6%).
100 samples accepted (2.6%).
102 samples accepted (2.6%).
106 samples accepted (2.6%).
110 samples accepted (2.7%).
113 samples accepted (2.7%).
116 samples accepted (2.7%).
121 samples accepted (2.8%).
126 samples accepted (2.8%).
128 samples accepted (2.8%).
130 samples accepted (2.8%).
133 samples accepted (2.8%).
135 samples accepted (2.8%).
137 samples accepted (2.7%).
139 samples accepted (2.7%).
140 samples accepted (2.7%).
143 samples accepted (2.7%).
145 samples accepted (2.7%).
147 samples accepted (2.7%).
150 samples accepted (2.7%).
153 samples accepted (2.7%).
154 samples accepted (2.7%).
158 samples accepted (2.7%).
161 samples accepted (2.7%).
166 samples accepted (2.7%).
169 samples accepted (2.7%).
170 samples accepted (2.7%).
172 samples accepted (2.7%).
177 samples accepted (2.7%).
180 samples accepted (2.7%).
183 samples accepted (2.7%).
185 samples accepted (2.7%).
190 samples accepted (2.8%).
192 samples accepted (2.7%).
194 samples accepted (2.7%).
195 samples accepted (2.7%).
197 samples accepted (2.7%).
198 samples accepted (2.7%).

In [29]:

    

pmc_posterior['epsilon'][-1]


    

    
Out[29]:





0.055113344227209454

In [30]:

    

fig, axes = plt.subplots(4, 1, figsize=(10,10))

for i,ax in enumerate(axes):
    ax.hist(pmc_posterior[-1][0][i,:], histtype='step', lw=3)
    ax.axvline(theta_0[i], color='r', ls=':')


    

In [38]:

    

from scipy.stats import beta
fig, ax = plt.subplots(1,1)
eccs = np.linspace(0,1,100)
for a,b in zip(pmc_posterior[i][0][2,:], pmc_posterior[i][0][3,:]):
    ax.plot(eccs, beta(a,b).pdf(eccs), color='k', alpha=0.1)


    

In [36]:

    

%matplotlib inline
import matplotlib.pyplot as plt

i1 = 2
i2 = 3

for i in range(pmc_posterior.shape[0]):
    plt.figure()
    plt.plot(pmc_posterior[i][0][i1,:], pmc_posterior[i][0][i2,:], '.');
    plt.title('epsilon = {}'.format(pmc_posterior[i]['epsilon']))
    plt.plot(theta_0[i1], theta_0[i2], 'rx', ms=20)
    plt.xlim(*model.bounds[i1])
    plt.ylim(*model.bounds[i2])


    

In [ ]: