working on sampling script

models/galaxy_model/run_sampler.py

@@ -6,6 +6,11 @@ import numpy as np
 from scipy.stats import gaussian_kde as kde
 import matplotlib.pyplot as plt
 
+import mergerratemodel as mr
+import prior_check 
+
+from PTMCMCSampler.PTMCMCSampler import PTSampler as ptmcmc
+
 
 def readChains(nFreqBins):
 
@@ -27,7 +32,7 @@ def makeKDEs(samples,nFreqBins):
     make the KDEs to use in the likelihood
     '''
    
-    KDE = {i: kde(samples[:,i]) for i in range(nFrequencies)}
+    KDE = {i: kde(samples[:,i]) for i in range(nFreqBins)}
 
     return KDE
 
@@ -38,7 +43,24 @@ def in_bounds(par):
 
 def log_prior(par):
     if in_bounds(par):
-        param = dict(Phi0 = par[0], PhiI = par[1], M0 = par[2], alpha = par[3], alphaI = par[4], f0 = par[5], beta = par[6], gamma = par[7], delta = par[8], t0 = par[9], epsilon = par[10], zeta = par[11], eta = par[12], Ms = par[13], theta = par[14], sigma = par[15], e0 = par[16], rho = par[17])        
+        param = dict(Phi0 = par[0], 
+                     PhiI = par[1], 
+                     M0 = par[2], 
+                     alpha = par[3], 
+                     alphaI = par[4], 
+                     f0 = par[5], 
+                     beta = par[6], 
+                     gamma = par[7], 
+                     delta = par[8], 
+                     t0 = par[9], 
+                     epsilon = par[10], 
+                     zeta = par[11], 
+                     eta = par[12], 
+                     Ms = par[13], 
+                     theta = par[14], 
+                     sigma = par[15], 
+                     e0 = par[16], 
+                     rho = par[17])        
         return prior_check.check_p(param)
     else:
         return -np.inf
@@ -49,7 +71,7 @@ def log_likelihood(par):
     M1 = np.linspace(9,12,25)
     q = np.linspace(0.25,1,10)
     z = np.linspace(0.,1.5,5)
-
+    f = frequencies
     initpar = dict(Phi0 = par[0], 
                    PhiI = par[1], 
                    M0 = par[2], 
@@ -71,20 +93,48 @@ def log_likelihood(par):
     model = mr.mergerrate(M1,q,z,f,**initpar).hmodelt(fbin=None)[0]
 
     logL = 0.
-    for i in range(nFrequencies):
-        logL += KDE[i].logpdf(model)
+    for i in range(numberOfFrequencyBinsToUse):
+        logL += KDEs[i].logpdf(model[i])
     return logL
 
 
 
-
+# prior bounds
+bounds=[(-3.4,-2.4),
+        (-0.6,0.2),
+        (11,11.5),
+        (-1.5,-1.),
+        (-0.2,0.2),
+        (0.01,0.05),
+        (0.,2.),
+        (-0.5,0.5),
+        (-0.2,0.2),
+        (0.1,10.),
+        (-0.5,0.5),
+        (-3.,1.),
+        (-0.2,0.2),
+        (7.75,8.75),
+        (0.9,1.1),
+        (0.2,0.5),
+        (0.01,0.99),
+        (-2.,2.)]
+
+# read in the data
 numberOfFrequencyBinsToUse = 5
 frequencies, samples = readChains(numberOfFrequencyBinsToUse)
+
+# make the kdes
 KDEs = makeKDEs(samples, numberOfFrequencyBinsToUse)
 
+# ptmcmc things
+x0 = np.array([-2.8,-0.2,11.25,-1.25,0.,0.025,0.8,0.,0.,1.,0.,-0.5,0.,8.25,1.,0.4,0.5,0.])
+ndim = len(x0)
+cov = np.diag(np.ones(ndim) * 0.01**2)
+N = 1000 # 1000000
 
+sampler = ptmcmc(ndim, log_likelihood, log_prior, cov, outDir='./output', resume=False)
+sampler.sample(x0, N, SCAMweight=30, AMweight=15, DEweight=50)
 
-to do - add sampler