agn model with cpnest

models/agnostic_model/run_agn_cpnest.py

+import numpy as np
+import cpnest.model
+from scipy.stats import gaussian_kde as kde
+import matplotlib.pyplot as plt
+import hmodel
+
+
+def readChain():
+    '''
+    read in the chains
+    '''
+    samples = np.genfromtxt('../../chains/dr2new_crn_fixed/chains_h.txt')
+    return samples
+
+
+def makeKDE():
+    '''
+    make the KDEs to use in the likelihood
+    '''
+    samples = readChain()
+    KDE = kde(samples)
+    return KDE
+
+
+def priorBounds():
+	
+    logno = [-20.0,3.0]
+    beta  = [-2.0,7.0]
+    gamma = [0.2,5.0]
+    alpha = [-3.0,3.0]
+    delta = [-1.0,2.0]
+
+    priors = [logno,beta,gamma,alpha,delta]
+
+    return priors
+
+
+class kdeModel(cpnest.model.Model):
+
+
+
+
+    # parameter names and get flat prior bounds
+    names=['logn','beta','gamma','alpha','delta']
+    bounds=priorBounds()
+    
+    
+    def log_likelihood(self,param):
+    
+        # integration limits
+        MlowIntLimit  = 10.0**6.0	
+        MhighIntLimit = 10.0**11.0
+        zlowIntLimit  = 0.0 
+        zhighIntLimit = 5.0
+
+        # calculate h
+        theta = param['logn'],param['beta'],param['gamma'],param['alpha'],param['delta']
+        hc = hmodel.hmodel(theta,\
+                           MlowIntLimit,MhighIntLimit,\
+                           zlowIntLimit,zhighIntLimit)
+                                         
+        log10hModel = np.log10(float(hc))
+        
+        log_like = - kdeLike.logpdf(log10hModel)        
+        
+        return log_like
+        
+
+
+        
+kdeLike = makeKDE()
+mymodel = kdeModel()
+nest = cpnest.CPNest(mymodel,maxmcmc=1000,nlive=10000,verbose=3)
+nest.run()
+cpnest.CPNest.get_posterior_samples(nest)       
+