notebook.community

Edit and run



In [32]:

    
import numpy as np
import pandas as pd
import emcee 
from sklearn.ensemble import ExtraTreesRegressor, RandomForestRegressor
from sklearn import tree
from sklearn.model_selection import cross_val_score
import re

%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
import corner



In [2]:

    
DF = pd.read_table('simulations.dat', sep='\t')
DF['age'] = DF['age'] / 10**9
DF.head()









    Out[2]:







  
    
      
      id
      M
      Y
      Z
      alpha
      diffusion
      age
      beta
      Teff
      L
      ...
      r01_43
      r13_44
      r10_43
      nu_0_45
      nu_2_44
      r01_44
      nu_3_8
      nu_2_9
      nu_3_44
      nu_2_8
    
  
  
    
      0
      20001
      0.776883
      0.285652
      0.007474
      0.720276
      1
      8.398173
      -0.025153
      5022.223586
      0.546163
      ...
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
    
    
      1
      20002
      0.811883
      0.225652
      0.015805
      2.220276
      1
      6.198423
      0.024847
      4830.176155
      0.270065
      ...
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
    
    
      2
      20003
      0.741883
      0.265652
      0.003534
      1.220276
      1
      10.597923
      -0.075153
      5530.540668
      0.533502
      ...
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
    
    
      3
      20004
      0.759383
      0.235652
      0.005139
      1.970276
      1
      5.098548
      -0.050153
      5324.783056
      0.351822
      ...
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
    
    
      4
      20007
      0.724383
      0.295652
      0.010868
      1.470276
      1
      7.298298
      0.099847
      5081.494960
      0.336307
      ...
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
      NaN
    
  

5 rows × 300 columns



In [3]:

    
DF.shape









    Out[3]:





(12720, 300)



In [4]:

    
KIC = pd.read_table('../regression/perturb/bigG/7970740_perturb.dat')
KIC = KIC.drop(['Dnu0', 'dnu02', 'r02', 'r01', 'r10'], axis=1)
KIC = KIC.drop([i for i in KIC.columns 
    if re.search('nu_|r01_', i)], axis=1)
KIC_columns = KIC.columns.values
KIC_columns[1] = 'Fe_H'
KIC.columns = KIC_columns
KIC.head()









    Out[4]:







  
    
      
      Teff
      Fe_H
      r02_17
      r02_18
      r02_19
      r02_20
      r02_21
      r02_22
      r02_23
      r02_24
      ...
      r10_19
      r10_20
      r10_21
      r10_22
      r10_23
      r10_24
      r10_25
      r10_26
      r10_27
      r10_28
    
  
  
    
      0
      5309.000000
      -0.540000
      0.063186
      0.056338
      0.054542
      0.055206
      0.049283
      0.049298
      0.052405
      0.036912
      ...
      0.015446
      0.015210
      0.014888
      0.014634
      0.014336
      0.015857
      0.014204
      0.010339
      0.007882
      0.019809
    
    
      1
      5406.247480
      -0.572623
      0.064753
      0.057509
      0.054785
      0.056270
      0.049662
      0.048558
      0.050471
      0.036498
      ...
      0.015442
      0.015770
      0.014949
      0.014335
      0.014835
      0.015585
      0.014184
      0.010739
      0.002031
      0.019091
    
    
      2
      5291.495691
      -0.513386
      0.062448
      0.055563
      0.054344
      0.056732
      0.049226
      0.049792
      0.050888
      0.037943
      ...
      0.015297
      0.015168
      0.014699
      0.014357
      0.013725
      0.016313
      0.014485
      0.009978
      0.005982
      0.018987
    
    
      3
      5326.634508
      -0.440346
      0.069188
      0.052609
      0.054532
      0.053740
      0.048868
      0.048610
      0.054871
      0.034841
      ...
      0.015781
      0.015021
      0.014749
      0.014997
      0.013261
      0.014480
      0.014027
      0.011441
      0.010553
      0.019370
    
    
      4
      5274.490584
      -0.492826
      0.063603
      0.056196
      0.056307
      0.054735
      0.048241
      0.047175
      0.052728
      0.041067
      ...
      0.015662
      0.015422
      0.015030
      0.014258
      0.013659
      0.016100
      0.015520
      0.011253
      0.012393
      0.028254
    
  

5 rows × 28 columns



In [5]:

    
y_star = KIC.iloc[0]
y_star









    Out[5]:





Teff      5309.000000
Fe_H        -0.540000
r02_17       0.063186
r02_18       0.056338
r02_19       0.054542
r02_20       0.055206
r02_21       0.049283
r02_22       0.049298
r02_23       0.052405
r02_24       0.036912
r02_25       0.035124
r02_26       0.028170
r02_27       0.014598
r02_28       0.034549
r02_29       0.030393
r10_16       0.020654
r10_17       0.019698
r10_18       0.018950
r10_19       0.015446
r10_20       0.015210
r10_21       0.014888
r10_22       0.014634
r10_23       0.014336
r10_24       0.015857
r10_25       0.014204
r10_26       0.010339
r10_27       0.007882
r10_28       0.019809
Name: 0, dtype: float64



In [6]:

    
Sigma = np.cov(KIC.iloc[:,2:].T)
Sigma = np.vstack([np.zeros(Sigma.shape[0]), np.zeros(Sigma.shape[0]), Sigma])
Sigma = np.insert(Sigma, 0, np.zeros(Sigma.shape[0]), axis=1)
Sigma = np.insert(Sigma, 0, np.zeros(Sigma.shape[0]), axis=1)
Sigma[0,0] = 77**2
Sigma[1,1] = 0.1**2

np.set_printoptions(precision=10)
print(np.sqrt(np.diag(Sigma)))

plt.imshow(Sigma[2:,2:], interpolation='nearest')
plt.show()









    



[7.7000000000e+01 1.0000000000e-01 3.9029146453e-03 2.1725166110e-03
 1.4763450735e-03 1.0333801995e-03 1.1141153559e-03 9.9657954298e-04
 1.5794116586e-03 3.0378212737e-03 4.6829660573e-03 6.9970515064e-03
 5.1684757992e-03 7.7142395021e-03 8.8689859465e-03 8.2823371490e-04
 5.6268399444e-04 4.4734607701e-04 4.0459075105e-04 4.5440933645e-04
 5.3669735521e-04 5.7708904568e-04 8.9493725612e-04 1.2624834117e-03
 1.7150952707e-03 2.5820145538e-03 3.9392679640e-03 5.1059887744e-03]



In [7]:

    
Sigma_inv = np.linalg.inv(Sigma)
Sigma_inv

plt.imshow(Sigma_inv, interpolation='nearest')
plt.show()



In [8]:

    
Sigma_inv.shape









    Out[8]:





(28, 28)



In [9]:

    
# train a forest to learn Teff, FeH, nu_max, Dnu0, dnu02 = f(age, M, Y, Z, alpha, overshoot, undershoot, diffusion)
y_columns = KIC.columns.tolist() #['Teff', 'Fe_H', 'nu_max', 'Dnu0', 'dnu02', 'r02', 'r01', 'r10']
X_columns = ['age', 'M', 'Y', 'Z', 'alpha', 'beta']
X_labels = [r'$\tau$', "$M$", "$Y$", "$Z$", r'$\alpha$', r'$\beta$']
#data_ = DF.loc[:,list(set(X_columns+y_columns))].dropna()
data_ = DF.loc[:,X_columns+y_columns].dropna()
data_.shape









    Out[9]:





(9206, 34)



In [10]:

    
X = data_.loc[:, X_columns]
X.head()



In [11]:

    
ys = data_.loc[:, [i for i in y_columns if i not in X_columns]]
ys.head()









    Out[11]:







  
    
      
      Teff
      Fe_H
      r02_17
      r02_18
      r02_19
      r02_20
      r02_21
      r02_22
      r02_23
      r02_24
      ...
      r10_19
      r10_20
      r10_21
      r10_22
      r10_23
      r10_24
      r10_25
      r10_26
      r10_27
      r10_28
    
  
  
    
      0
      5022.223586
      -0.490474
      0.078714
      0.076353
      0.073776
      0.071417
      0.068856
      0.066595
      0.064218
      0.062057
      ...
      0.027489
      0.025164
      0.024356
      0.022492
      0.021573
      0.020378
      0.019242
      0.018492
      0.017306
      0.016905
    
    
      2
      5530.540668
      -0.796330
      0.069752
      0.066596
      0.063632
      0.060687
      0.058197
      0.055656
      0.053445
      0.051366
      ...
      0.019706
      0.019338
      0.017466
      0.016923
      0.016058
      0.015088
      0.015102
      0.014103
      0.014194
      0.013571
    
    
      4
      5081.494960
      -0.234200
      0.068950
      0.066202
      0.063436
      0.061174
      0.058782
      0.056601
      0.054814
      0.052832
      ...
      0.020414
      0.018690
      0.017891
      0.018260
      0.016495
      0.016631
      0.016484
      0.014804
      0.015465
      0.014555
    
    
      5
      4966.956239
      -0.741142
      0.069361
      0.066217
      0.063178
      0.060358
      0.057819
      0.055288
      0.053120
      0.050975
      ...
      0.019748
      0.018927
      0.017172
      0.016857
      0.015553
      0.015067
      0.014783
      0.013872
      0.014243
      0.013244
    
    
      7
      5409.063014
      -0.406669
      0.051582
      0.047466
      0.043100
      0.039070
      0.035389
      0.031672
      0.028385
      0.025322
      ...
      0.022608
      0.020189
      0.015852
      0.014096
      0.011631
      0.008566
      0.007720
      0.005432
      0.003902
      0.003571
    
  

5 rows × 28 columns



In [44]:

    
forest = RandomForestRegressor(n_estimators=8)
print(cross_val_score(forest, X, ys, cv=3))
forest.fit(X, ys)
#forest.oob_score_









    



[0.9723 0.9736 0.9747]






    Out[44]:





RandomForestRegressor(bootstrap=True, criterion='mse', max_depth=None,
           max_features='auto', max_leaf_nodes=None,
           min_impurity_decrease=0.0, min_impurity_split=None,
           min_samples_leaf=1, min_samples_split=2,
           min_weight_fraction_leaf=0.0, n_estimators=8, n_jobs=None,
           oob_score=False, random_state=None, verbose=0, warm_start=False)



In [ ]:

    
#forest = tree.DecisionTreeRegressor(
#    n_estimators=8,
#    n_jobs=1,
#    oob_score=1, bootstrap=1)
#forest.fit(X, ys)
#forest.oob_score_



In [39]:

    
#forest = RandomForestRegressor(n_estimators=16)
#chi2s = [-np.log10(np.dot(np.dot(y-y_star, Sigma_inv), y-y_star)/2) for idx, y in ys.iterrows()]
#chi2s = [-np.log(np.dot(np.dot(y-y_star, Sigma_inv), y-y_star)/2) for idx, y in ys.iterrows()]
#print(cross_val_score(forest, X, ys, cv=3))
#forest.fit(X, chi2s)
#forest.oob_score_









    



[0.9758 0.9776 0.9786]



In [45]:

    
X_max = X.max()
X_min = X.min()

def lnprior(theta):
    if any(theta > X_max) or any(theta < X_min): #X.max()
        return -np.inf
    return 0.0

def lnlike(theta):
    y = forest.predict([theta])[0]
    resid = y-y_star
    chi2 = np.log(np.dot(np.dot(resid, Sigma_inv), resid))
    return(-chi2/2)
    #return forest.predict([theta])[0]

def lnprob(theta):
    lp = lnprior(theta)
    if not np.isfinite(lp):
        return -np.inf
    return lp + lnlike(theta)



In [46]:

    
theta = [9.9, 0.75408, 0.258478, 0.005, 1.85, -0.0017]

lnprob(theta)









    Out[46]:





-2.7240912276571105



In [49]:

    
ndim, nwalkers = 6, 20
pos = [[10., 0.77, 0.26, 0.005, 1.8, 0.] + 1e-3*np.random.randn(ndim) 
    for i in range(nwalkers)]

np.set_printoptions(precision=4, linewidth=100)
pos









    Out[49]:





[array([ 9.9995e+00,  7.7107e-01,  2.6011e-01,  6.5588e-03,  1.8011e+00, -1.6479e-04]),
 array([ 9.9987e+00,  7.6927e-01,  2.5912e-01,  6.4745e-03,  1.7995e+00, -1.8838e-03]),
 array([ 1.0000e+01,  7.7142e-01,  2.5814e-01,  5.8477e-03,  1.8003e+00, -4.8730e-04]),
 array([1.0002e+01, 7.6943e-01, 2.6018e-01, 4.4436e-03, 1.7993e+00, 6.1778e-04]),
 array([9.9993e+00, 7.6885e-01, 2.5918e-01, 5.4226e-03, 1.8014e+00, 2.8782e-03]),
 array([9.9990e+00, 7.7056e-01, 2.6008e-01, 4.3258e-03, 1.7985e+00, 6.7903e-04]),
 array([ 9.9988e+00,  7.6951e-01,  2.6018e-01,  3.6671e-03,  1.7994e+00, -1.2159e-03]),
 array([ 9.9993e+00,  7.6977e-01,  2.6090e-01,  5.1890e-03,  1.8000e+00, -8.5293e-04]),
 array([9.9987e+00, 7.7096e-01, 2.5979e-01, 5.4510e-03, 1.7991e+00, 1.9193e-04]),
 array([1.0000e+01, 7.6955e-01, 2.5952e-01, 5.3777e-03, 1.7997e+00, 7.6795e-04]),
 array([1.0000e+01, 7.6854e-01, 2.5905e-01, 4.2826e-03, 1.7993e+00, 5.9654e-05]),
 array([ 1.0003e+01,  7.7110e-01,  2.6168e-01,  7.0142e-03,  1.7983e+00, -9.0540e-04]),
 array([ 1.0000e+01,  7.7066e-01,  2.6163e-01,  5.5054e-03,  1.7973e+00, -5.8777e-04]),
 array([1.0002e+01, 7.6825e-01, 2.5891e-01, 4.8120e-03, 1.8001e+00, 1.9693e-03]),
 array([ 1.0001e+01,  7.7005e-01,  2.5991e-01,  5.4813e-03,  1.7998e+00, -3.7214e-04]),
 array([ 9.9989e+00,  7.7038e-01,  2.6336e-01,  5.2743e-03,  1.8002e+00, -5.5185e-04]),
 array([ 1.0001e+01,  7.7164e-01,  2.5927e-01,  5.8817e-03,  1.8005e+00, -8.6105e-04]),
 array([ 9.9999e+00,  7.7113e-01,  2.6094e-01,  5.2491e-03,  1.8007e+00, -1.6372e-03]),
 array([ 9.9985e+00,  7.7002e-01,  2.5923e-01,  3.8643e-03,  1.7996e+00, -6.8963e-04]),
 array([ 1.0001e+01,  7.6910e-01,  2.5995e-01,  5.4106e-03,  1.8003e+00, -9.1566e-04])]



In [54]:

    
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob)
sampler.run_mcmc(pos, 5000)









    Out[54]:





(array([[ 7.0529e+00,  7.7028e-01,  2.4726e-01,  1.3589e-02,  1.5115e+00, -9.9424e-03],
        [ 6.9510e+00,  8.0012e-01,  2.8381e-01,  5.5025e-03,  1.8209e+00,  8.0135e-02],
        [ 7.6644e+00,  7.0210e-01,  2.7763e-01,  1.6693e-02,  1.4561e+00, -8.6081e-02],
        [ 5.3017e+00,  8.1491e-01,  2.8350e-01,  1.9214e-02,  2.3769e+00,  5.2905e-02],
        [ 1.1439e+01,  8.2289e-01,  2.7491e-01,  1.3977e-02,  2.1475e+00,  6.9682e-02],
        [ 7.8509e+00,  7.4709e-01,  2.8434e-01,  8.0869e-03,  1.8550e+00, -4.8550e-02],
        [ 9.4714e+00,  8.1814e-01,  2.2582e-01,  6.0953e-03,  2.1017e+00,  1.5761e-02],
        [ 1.3121e+01,  7.9407e-01,  2.5964e-01,  1.2943e-02,  1.7088e+00, -4.1017e-02],
        [ 1.2065e+01,  7.2055e-01,  2.6183e-01,  1.3980e-02,  2.2677e+00, -9.5346e-02],
        [ 1.2744e+01,  7.2584e-01,  2.8385e-01,  7.3300e-03,  1.4274e+00, -8.3929e-02],
        [ 1.3254e+01,  7.3514e-01,  2.8038e-01,  2.7304e-03,  2.1079e+00,  1.7154e-02],
        [ 1.2153e+01,  7.5131e-01,  2.3756e-01,  1.7001e-02,  1.5791e+00, -9.4401e-02],
        [ 9.1163e+00,  7.0101e-01,  2.3986e-01,  5.9385e-03,  1.9669e+00, -1.2298e-02],
        [ 1.1653e+01,  7.3532e-01,  2.7866e-01,  1.7501e-02,  2.3424e+00,  3.9572e-02],
        [ 6.3977e+00,  7.1625e-01,  2.8269e-01,  1.7570e-02,  2.4988e+00, -3.7899e-02],
        [ 1.0889e+01,  7.5802e-01,  2.9897e-01,  1.7149e-02,  8.8151e-01, -2.3228e-02],
        [ 1.0184e+01,  7.7356e-01,  2.2769e-01,  1.0570e-02,  1.3523e+00, -4.5247e-02],
        [ 1.1108e+01,  8.0445e-01,  2.8348e-01,  1.3175e-02,  1.5186e+00,  5.2871e-05],
        [ 8.0044e+00,  8.3535e-01,  2.8986e-01,  1.2907e-02,  2.2849e+00,  6.8093e-02],
        [ 8.3073e+00,  7.1284e-01,  2.9724e-01,  1.8768e-02,  9.5310e-01,  4.7513e-02]]),
 array([-2.4821, -3.4169, -2.5492, -2.5202, -3.0175, -2.5611, -2.8924, -2.4405, -3.1312, -2.964 ,
        -3.4978, -2.7008, -2.4936, -3.4679, -3.1688, -3.1036, -2.8441, -2.9872, -2.7408, -3.3687]),
 ('MT19937',
  array([2610935452, 1910583459, 4271360869, 2076275873, 1954976454, 1548584202, 1502570744,
          473532763, 1817730251, 1470782210, 1167316908, 4126075200, 2943902468, 4073404919,
           49727157, 1269502416,   18021379,  782575297, 3866575697, 3966901029,  869883991,
         1311213367,   25348983, 2634864539,  222190038, 2303035182, 3352207642, 1301117068,
         3774438861, 1243067828, 3571927162,  862809450,  873549925, 2554296548, 3150529068,
           67684396, 1084437297, 4163012200,  749310372, 3648035334,  299228144,  637906932,
          175048994,  592483393, 3171433140, 1460050393, 1277927670, 2257019689, 3670090600,
          338917966, 4083654068, 2620061764,  683424177, 3838037200, 1969507765, 1821917278,
         1251347103, 3380376332, 3874410771, 1536839600, 1062155831, 2372723738, 1110798607,
           49755852, 1136450178, 1738957999,  935532952,  888302506, 1074445999, 1808300512,
         3326007625, 1257115254, 1726651039, 1276538812, 2746736594, 1549021448, 2486937579,
         1413597135, 2095649471, 3809658198, 3094216100, 1750695071, 1149237860,  342053660,
         2160813682,  807709573, 2782764469,  523898277, 2827059685, 1381076377, 2038677466,
           89649236, 1128054915, 1589969158, 3364912494,  671202571,  511788423, 2972235247,
         3706450898, 1433860410, 2227207292,  477303241, 4003172497,  803511402, 2298673679,
         1132759884, 3489882459, 2452156097,  292605914, 3486507923, 1066880046, 3218845662,
          350249474, 2674322878, 4095623145, 3456986158, 3035201532,  944863439, 2080949159,
         1499436501, 1819405493,  368049202, 1669608210, 3549965963, 1274844594, 3917375582,
          183939899, 2947564999,   85886564, 2677972035,  578294624, 1741126030, 2051625773,
         3394499362,  167948159, 1505462422,  410034818, 1384430734, 4125235928, 1298115391,
         3912782638,  911649870, 4143105791, 2029040158,  265947074, 3130378965, 1196904836,
          805907446,  427805389, 2496037424, 1076960348, 2752067122, 3889874272, 3404659022,
         1856543246, 3714538117, 2246066462,  477119905,  899822943,  768219927,  546735334,
         4062907292, 3139551432, 1743629436,   33989388, 3701754005,  998768265, 1634311174,
         3243624803, 4210539573, 4133530876, 2044916770, 3250597166, 2260308501, 1757981303,
         3522030351, 2159575393, 1547201716, 1133024116,  514514747, 2059311210, 1774632042,
         2487140997,  435361845,  792963047, 2163203760, 3796974831, 4017513457,  277659844,
         3128957337, 1679359694,  403481009, 1723440482, 4144540147, 2245257862, 2590958118,
         1500066265, 1752036958, 3436246022,  434073346, 2053206975, 2924891929, 1731528759,
         2852682296, 3929596978, 3330197373, 1879318627, 3328769665, 3311574334, 4137320990,
         1161323367, 1154374774,  888735115, 3873908869, 3670846772, 2921570578, 3497102775,
          429848046,  621957283,  375192169, 2838734178, 2756204571, 3497029304, 2783075915,
         1160370487,  551454837, 4049438609, 1555272616, 2158150096,  771586153, 1069214016,
          441770037, 2390357198, 1887360126, 3882811137, 2793379239,   16734769, 2986321487,
          825273959, 1616631115, 2897932156, 3580441777, 2979142439, 1807617332, 2726519776,
         2384935448,  139681711, 1391516931, 4283265868, 2988818846, 3292770713, 3143660581,
         2145722267, 4138640187, 1914510175, 3605920814, 3109442920,  909777137, 2395593644,
         3421574355, 2256144228,  256244776, 3334345627, 2765498582,  242383193, 1113079150,
         2806252010, 4070742897, 1335726469, 2529661568, 3856345562, 1355908116,  877031621,
          728905860, 1934763300,  187726603, 3405733148, 3254978368,   74629631, 3801574889,
         3879623092, 2105474350, 2573303119, 2860492838, 3806097181, 1709531039, 3024630331,
         3098747125, 2639848162,  935023572, 1247393473, 4089488857, 1020472565, 2192692703,
         2650171548, 3909078113, 1444385444,  562200360, 2065908642, 3091231131,  576290545,
         3284338441, 1589923610, 3800511656, 1243088650, 1597897427, 1623641247, 1506244999,
          394174234, 3354170668, 1864860880, 1165639809, 3049756739, 1111198373, 1802797635,
          786835340,  125583303,  892690282, 1304495784,  300346054, 4018290157, 2153620450,
          302809858, 3275078418, 2631167561, 2970547582, 1365891914, 2287699077, 2134645251,
         1802975981, 2086011797, 1657793138,  861811487, 2970498865, 1483840131, 3657096761,
         4270718873, 4115049655, 3200021204, 1802412310, 2903487552, 2202572046, 4253391718,
         1317726670, 1255442647, 3441784630, 2564873840, 3576816952, 4293470698, 1356524226,
          546974670,  915334200, 2797555021, 3036738296,  129245447,  376764907, 2165419543,
         3764407751, 3860061533, 1850102221, 2870945150, 2286669520,  934992034, 4153442551,
         2443786624, 1420529451, 3468968345, 1431214427, 4054552918, 2261085105, 4021883509,
         2691786076, 3659339002, 2297382186,  927571360, 2428058077, 2679732548,  795905350,
         2748745057, 1121267393,  293101296,  362074778, 2733866053, 1522044943, 3202295123,
          356068938, 3030054666,  212734273, 4100664465, 3893577462, 2240855289, 1467712794,
          950172830, 1350231732, 2236372138, 1387491206, 1632943121, 1922122862, 3682520919,
         2738240165, 1025043826, 3418330049, 3988889897, 1592151443, 3077703002, 1449550609,
         3686196297, 1621504853, 2667575274, 2308093113,  885975401, 3888912252, 1762480590,
          810993515,  892849123, 1281927368, 1507715747,  864894476, 3013991588, 3883103927,
         3894851201, 3945460875, 1679780287, 2175863919, 2091406691,  917582552, 1514104827,
         1429025632, 1586382954, 2224237547, 4224918627, 4115541085,  767994204,  456304029,
         2705708920, 3219419772,  528766665,  743873048, 3388230198, 4053021625, 2693330589,
         2767848252, 1983149104, 2804700135, 3053656578, 3512396556, 1360325010, 2800901450,
         2442735218,  816937037, 1113852980,  370266618, 1761637500, 2001477124, 1318282552,
          913373217, 2802612408, 1965033824, 1858653854,   79252024, 1507182783, 3287484502,
          756455868, 3423424163,  208989050, 2931292386, 3876312884, 2630701613, 2123095812,
         2873026534, 4165767903, 3208370400,   83547191, 2916284346, 4050841486, 2303508063,
          928934701,  792125598,  866482632,   41023483, 3645281231, 3555046603,  578912021,
         2845369066,  681700085, 3336529158, 2110986177, 3915139299, 2870662277, 1729797879,
         1693169602,  339701923, 3191987117, 2384321002, 3069771509, 1556935283, 2471333871,
         4161537204, 1688119399, 1278885885, 3527632774, 1139867913, 2517832617, 2316631005,
         3566635273,  966487573, 3116555354, 1703982660,  181572616, 2814879310, 3618538923,
         3189703127, 2617528917,  534199988, 2858385507,  503300948, 2094558387, 3865485863,
          989706907, 1206753212, 2686063040, 2626064309, 2225138634, 2538868949, 2198513654,
         4209772453, 1874215163, 3734880776,  167241343, 1652036652, 1137744460,   98944193,
         2042534442, 2785821151, 2187032153, 1813256806,  675115240, 1649898846, 1980461952,
         1965543096, 2372937861, 3711812832, 4089093142, 3419872436, 1195369905, 3004220952,
         2094829955, 1541054688, 1564291503, 1398102383, 2514594820, 1872937961, 1204997903,
         2694418507, 2737963955, 2001570442, 4082736508, 2103106672,  164291855, 1643686788,
         2989569287,   53723015,  152547637, 2877148283,  224913893,  842589204, 3458814019,
         2311561545, 1821792427, 4112419001, 3201434177, 1219231100, 2357809396, 4021614955,
          956935366, 3299750452, 2998328112,  880087209, 2141473450, 4084356849, 4091418007,
         3058012003, 1425004079,  939067842, 1071124662, 3099924445,  792948297, 3134256742,
         2585731683, 2941597337,  523862506,  314429881, 4059380751,  410141839, 1533732045,
          482125507,  794197431, 1671533722,  147564756, 4242184193, 3011511146, 2650942629,
         4283340331, 3340193083,  350351991, 3708407955, 3400981769, 1630725221, 3918001174,
         1914302661, 2160522965, 3783517780,  609426896, 1794361627, 1629695840,  443739409,
         3590114635, 1738918797, 2875891530,  389879194, 2892557446,  998976603, 2046994627,
          748382047], dtype=uint32),
  246,
  0,
  0.0))



In [55]:

    
fig, axes = plt.subplots(ndim, figsize=(10, 7), sharex=True)
samples = sampler.chain
labels = X_labels
for i in range(ndim):
    ax = axes[i]
    ax.plot(samples[:, :, i], "k", alpha=0.3)
    ax.set_xlim(0, len(samples))
    ax.set_ylabel(labels[i])
    ax.yaxis.set_label_coords(-0.1, 0.5)

axes[-1].set_xlabel("step number");



In [52]:

    
samples = sampler.chain[:, 1000:, :].reshape((-1, ndim))



In [53]:

    
fig = corner.corner(samples, labels=X_labels)
fig.show()









    



---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-53-5439358bed5e> in <module>()
----> 1 fig = corner.corner(samples, labels=X_labels)
      2 fig.show()

~/.local/lib/python3.7/site-packages/corner/corner.py in corner(xs, bins, range, weights, color, smooth, smooth1d, labels, label_kwargs, show_titles, title_fmt, title_kwargs, truths, truth_color, scale_hist, quantiles, verbose, fig, max_n_ticks, top_ticks, use_math_text, hist_kwargs, **hist2d_kwargs)
    145         assert len(xs.shape) == 2, "The input sample array must be 1- or 2-D."
    146         xs = xs.T
--> 147     assert xs.shape[0] <= xs.shape[1], "I don't believe that you want more " \
    148                                        "dimensions than samples!"
    149 

AssertionError: I don't believe that you want more dimensions than samples!



In [ ]:

    
print(samples)
list(map(lambda v: (v[1], v[2]-v[1], v[1]-v[0]),
    zip(*np.percentile(samples, [16, 50, 84], axis=0))))



In [ ]:

	id	M	Y	Z	alpha	diffusion	age	beta	Teff	L	...	r01_43	r13_44	r10_43	nu_0_45	nu_2_44	r01_44	nu_3_8	nu_2_9	nu_3_44	nu_2_8
0	20001	0.776883	0.285652	0.007474	0.720276	1	8.398173	-0.025153	5022.223586	0.546163	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
1	20002	0.811883	0.225652	0.015805	2.220276	1	6.198423	0.024847	4830.176155	0.270065	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
2	20003	0.741883	0.265652	0.003534	1.220276	1	10.597923	-0.075153	5530.540668	0.533502	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
3	20004	0.759383	0.235652	0.005139	1.970276	1	5.098548	-0.050153	5324.783056	0.351822	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
4	20007	0.724383	0.295652	0.010868	1.470276	1	7.298298	0.099847	5081.494960	0.336307	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

	Teff	Fe_H	r02_17	r02_18	r02_19	r02_20	r02_21	r02_22	r02_23	r02_24	...	r10_19	r10_20	r10_21	r10_22	r10_23	r10_24	r10_25	r10_26	r10_27	r10_28
0	5309.000000	-0.540000	0.063186	0.056338	0.054542	0.055206	0.049283	0.049298	0.052405	0.036912	...	0.015446	0.015210	0.014888	0.014634	0.014336	0.015857	0.014204	0.010339	0.007882	0.019809
1	5406.247480	-0.572623	0.064753	0.057509	0.054785	0.056270	0.049662	0.048558	0.050471	0.036498	...	0.015442	0.015770	0.014949	0.014335	0.014835	0.015585	0.014184	0.010739	0.002031	0.019091
2	5291.495691	-0.513386	0.062448	0.055563	0.054344	0.056732	0.049226	0.049792	0.050888	0.037943	...	0.015297	0.015168	0.014699	0.014357	0.013725	0.016313	0.014485	0.009978	0.005982	0.018987
3	5326.634508	-0.440346	0.069188	0.052609	0.054532	0.053740	0.048868	0.048610	0.054871	0.034841	...	0.015781	0.015021	0.014749	0.014997	0.013261	0.014480	0.014027	0.011441	0.010553	0.019370
4	5274.490584	-0.492826	0.063603	0.056196	0.056307	0.054735	0.048241	0.047175	0.052728	0.041067	...	0.015662	0.015422	0.015030	0.014258	0.013659	0.016100	0.015520	0.011253	0.012393	0.028254

	Teff	Fe_H	r02_17	r02_18	r02_19	r02_20	r02_21	r02_22	r02_23	r02_24	...	r10_19	r10_20	r10_21	r10_22	r10_23	r10_24	r10_25	r10_26	r10_27	r10_28
0	5022.223586	-0.490474	0.078714	0.076353	0.073776	0.071417	0.068856	0.066595	0.064218	0.062057	...	0.027489	0.025164	0.024356	0.022492	0.021573	0.020378	0.019242	0.018492	0.017306	0.016905
2	5530.540668	-0.796330	0.069752	0.066596	0.063632	0.060687	0.058197	0.055656	0.053445	0.051366	...	0.019706	0.019338	0.017466	0.016923	0.016058	0.015088	0.015102	0.014103	0.014194	0.013571
4	5081.494960	-0.234200	0.068950	0.066202	0.063436	0.061174	0.058782	0.056601	0.054814	0.052832	...	0.020414	0.018690	0.017891	0.018260	0.016495	0.016631	0.016484	0.014804	0.015465	0.014555
5	4966.956239	-0.741142	0.069361	0.066217	0.063178	0.060358	0.057819	0.055288	0.053120	0.050975	...	0.019748	0.018927	0.017172	0.016857	0.015553	0.015067	0.014783	0.013872	0.014243	0.013244
7	5409.063014	-0.406669	0.051582	0.047466	0.043100	0.039070	0.035389	0.031672	0.028385	0.025322	...	0.022608	0.020189	0.015852	0.014096	0.011631	0.008566	0.007720	0.005432	0.003902	0.003571