Initial commit

REAME.txt

@@ -0,0 +1,8 @@
+
+Python module to implement Bayesian inference routines.
+
+Neale Gibson
+[email protected]
+[email protected]
+
+This module has merged my MCMC routines and Gaussian processes modules. The MCMC now includes orthogonal stepping, affine invariant methods (emcee), and I've made many additions to the Gaussian Process classes to include multiplicative systematics (not fully tested), easy integration to the inference and optimisation tools, better plotting interfaces, and more kernels. I've also added various methods to get the Bayesian evidence, eg importance sampling - best initiated from MCMC samples. Please contact me should you wish to use this code and have questions. These codes or earlier variations have been used for the inference in my recent papers.

examples/Aff_test1.py

@@ -0,0 +1,74 @@
+#!/usr/bin/env python
+
+import numpy as np
+import pylab
+import os
+
+import Infer
+import MyMCMC
+from MyFunctions import Transit_aRs
+from MyFunctions import LogLikelihood_iid_mf
+
+def Posterior(p,*args):
+  if p[-1] < 0.: return -np.inf
+  if np.any(p < 0.): return -np.inf
+  return LogLikelihood_iid_mf(p,*args)
+
+#light curve parameters
+lc_pars = [.0,2.5,11.,.1,0.6,0.2,0.3,1.,0.]
+wn = 0.0003
+
+#create the data set (ie training data)
+time = np.linspace(-0.1,0.1,300)
+flux = Transit_aRs(lc_pars,time) + np.random.normal(0,wn,time.size)
+
+#guess parameter values and guess uncertainties
+guess_pars = lc_pars + [wn]
+err_pars = [0.0001,0,0.05,0.003,0.002,0.0,0.0,0.01,0.0001,0.00010]
+
+#plot the light curve + guess function
+pylab.figure(1)
+pylab.errorbar(time,flux,yerr=wn,fmt='.')
+pylab.plot(time,Transit_aRs(guess_pars[:-1],time),'r--')
+
+#first optimise the function
+guess_pars = Infer.Optimise(Posterior,guess_pars[:],(Transit_aRs,time,flux),fixed=(np.array(err_pars) == 0)*1,method='BFGS')
+
+#run a normal MCMC
+chain_len = 60000
+conv = 30000
+thin = 10
+no_ch=3
+adapt_lims = (2000,conv,10)
+glob_lims = (2000,conv,10)
+# Infer.MCMC(Posterior,guess_pars[:],(Transit_aRs,time,flux),chain_len,err_pars,n_chains=no_ch,adapt_limits=adapt_lims,glob_limits=glob_lims,thin=thin)
+# #Get parameters values/errors from chains
+# par,par_err = Infer.AnalyseChains(conv/thin,n_chains=no_ch)
+# bf_par = Infer.GetBestFit(n_chains=no_ch)
+# print "Best Fit log p =", LogLikelihood_iid_mf(bf_par,Transit_aRs,time,flux)
+# pylab.figure(3)
+# Infer.PlotCorrelations(conv/thin,n_chains=no_ch,p=np.where(np.array(par_err)>0.)[0])
+
+#run an affine inv MCMC
+n = 300
+chain_len = chain_len/n
+conv = conv/n
+no_ch=3
+Infer.AffInvMCMC(Posterior,guess_pars[:],(Transit_aRs,time,flux),n,chain_len,err_pars,n_chains=no_ch)
+#Get parameters values/errors from chains
+par,par_err = Infer.AnalyseChains(conv*n,n_chains=no_ch)
+bf_par = Infer.GetBestFit(n_chains=no_ch)
+print "Best Fit log p =", LogLikelihood_iid_mf(bf_par,Transit_aRs,time,flux)
+pylab.figure(4)
+Infer.PlotCorrelations(conv*n,n_chains=no_ch,p=np.where(np.array(par_err)>0.)[0])
+
+#plot the chains and correlations
+#lab = [r'$T_0$',r'$a/R\star$',r'$\rho$',r'$b$']
+#pylab.figure(2)
+#Infer.PlotChains(conv/thin,n_chains=no_ch,p=[0,2,3,4],labels=lab)
+
+#plot fitted function
+pylab.figure(1)
+pylab.plot(time,Transit_aRs(par[:-1],time),'g-')
+
+raw_input()

examples/Aff_test2.py

@@ -0,0 +1,71 @@
+#!/usr/bin/env python
+
+import numpy as np
+import pylab
+import os
+
+import Infer
+import MyMCMC
+from MyFunctions import Transit_aRs
+from MyFunctions import LogLikelihood_iid_mf
+
+#define a test function: 2D Rosenbrock function (inverted)
+def InvRosenbrock(p):
+  R = (1.-p[0])**2. + 100.*(p[1]-p[0]**2.)**2.  
+  return -R
+
+#multivariate - sum of p/2 coupled Rosenbrock functions
+def InvNRosenbrock(p):
+  R = 0.
+  for i in range(len(p)/2):
+    R += 100.*(p[2*i]**2-p[2*i+1])**2. + (p[2*i]-1.)**2.
+  return -0.4*R
+
+gp = [1.,1.,1.,1.,1.,1.,1.,1.]
+ep = [0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1]
+
+gp = np.array(gp) + np.random.normal(0,1,len(gp))
+
+#first optimise the function
+gp = Infer.Optimise(InvNRosenbrock,gp[:],(),fixed=(np.array(ep) == 0)*1)
+
+#run a normal MCMC
+chain_len = 50000
+conv = 2000
+thin = 1
+no_ch=3
+adapt_lims = (200,conv,10)
+glob_lims = (200,conv,10)
+Infer.MCMC(InvNRosenbrock,gp[:],(),chain_len,ep,n_chains=no_ch,adapt_limits=adapt_lims,glob_limits=glob_lims,thin=thin)
+#Get parameters values/errors from chains
+par,par_err = Infer.AnalyseChains(conv/thin,n_chains=no_ch)
+bf_par = Infer.GetBestFit(n_chains=no_ch)
+print "Best Fit log p =", InvNRosenbrock(bf_par,)
+pylab.figure(3)
+Infer.PlotCorrelations(conv/thin,n_chains=no_ch,p=np.where(np.array(par_err)>0.)[0])
+
+#run an affine inv MCMC
+n = 200
+chain_len = chain_len/n
+conv = conv
+no_ch=3
+Infer.AffInvMCMC(InvNRosenbrock,gp[:],(),n,chain_len,ep,n_chains=no_ch)
+#Get parameters values/errors from chains
+par,par_err = Infer.AnalyseChains(conv,n_chains=no_ch)
+bf_par = Infer.GetBestFit(n_chains=no_ch)
+print "Best Fit log p =", InvNRosenbrock(bf_par,)
+pylab.figure(4)
+Infer.PlotCorrelations(conv*n,n_chains=no_ch,p=np.where(np.array(par_err)>0.)[0])
+
+# 
+# 
+# #plot the chains and correlations
+# #lab = [r'$T_0$',r'$a/R\star$',r'$\rho$',r'$b$']
+# #pylab.figure(2)
+# #Infer.PlotChains(conv/thin,n_chains=no_ch,p=[0,2,3,4],labels=lab)
+# 
+# #plot fitted function
+# pylab.figure(1)
+# pylab.plot(time,Transit_aRs(par[:-1],time),'g-')
+
+raw_input()

examples/BayesModelSelection.py

@@ -0,0 +1,79 @@
+#!/usr/bin/env python
+
+import numpy as np
+import pylab
+import os
+
+import Infer
+import MyMCMC
+from MyFunctions import LogLikelihood_iid_mf
+
+#define some polynomials
+def Model1(p,x):
+  return p[0]+p[1]*x
+def Model2(p,x):
+  return p[0]+p[1]*x+p[2]*x**2
+def Model3(p,x):
+  return p[0]+p[1]*x+p[2]*x**2+p[3]*x**3
+
+#generate some noisy data
+par1 = [0.,0.3,0.01]
+par2 = [0.2,0.2,-0.002,0.01]
+par3 = [0.2,0.2,-0.002,0.002,0.01]
+x = np.linspace(-5,5,100)
+y1 = Model1(par1,x) + np.random.normal(0,par1[-1],x.size)
+y2 = Model2(par2,x) + np.random.normal(0,par2[-1],x.size)
+y3 = Model3(par3,x) + np.random.normal(0,par3[-1],x.size)
+y = y3
+
+#plot the data
+pylab.errorbar(x,y1,yerr=par1[-1],ls='.')
+pylab.errorbar(x,y2,yerr=par2[-1],ls='.')
+pylab.errorbar(x,y3,yerr=par3[-1],ls='.')
+pylab.plot(x,Model1(par1[:-1],x),'r--')
+pylab.plot(x,Model2(par2[:-1],x),'r--')
+pylab.plot(x,Model3(par3[:-1],x),'r--')
+
+#MCMC params
+chain_len = 25000
+conv = 10000
+thin=10
+no_ch=2
+limits = (conv/4,conv,3)
+
+#first try the linear model
+gp = par1[:]
+ep = [0.001,0.001,0.001]
+Infer.MCMC(LogLikelihood_iid_mf,gp,(Model1,x,y),chain_len,ep,n_chains=no_ch,adapt_limits=limits,glob_limits=limits,thin=thin)
+par,par_err = Infer.AnalyseChains(conv/thin,n_chains=no_ch)
+pylab.plot(x,Model1(par[:-1],x),'b')
+m,K = Infer.NormalFromMCMC(conv/thin,n_chains=no_ch)
+E1,par,par_err = Infer.ImportanceSamp(LogLikelihood_iid_mf,(Model1,x,y),m,K,chain_len)
+E1,par,par_err = Infer.ImportanceSamp(LogLikelihood_iid_mf,(Model1,x,y),m,3.*K,chain_len)
+
+#now try the quadratic model
+gp = par2[:]
+ep = [0.001,0.001,0.001,0.001]
+Infer.MCMC(LogLikelihood_iid_mf,gp,(Model2,x,y),chain_len,ep,n_chains=no_ch,adapt_limits=limits,glob_limits=limits,thin=thin)
+par,par_err = Infer.AnalyseChains(conv/thin,n_chains=no_ch)
+pylab.plot(x,Model2(par[:-1],x),'g')
+m,K = Infer.NormalFromMCMC(conv/thin,n_chains=no_ch)
+E2,par,par_err = Infer.ImportanceSamp(LogLikelihood_iid_mf,(Model2,x,y),m,K,chain_len)
+E2,par,par_err = Infer.ImportanceSamp(LogLikelihood_iid_mf,(Model2,x,y),m,3.*K,chain_len)
+
+#3rd order...
+gp = par3[:]
+ep = [0.001,0.001,0.001,0.001,0.001]
+Infer.MCMC(LogLikelihood_iid_mf,gp,(Model3,x,y),chain_len,ep,n_chains=no_ch,adapt_limits=limits,glob_limits=limits,thin=thin)
+par,par_err = Infer.AnalyseChains(conv/thin,n_chains=no_ch)
+pylab.plot(x,Model2(par[:-1],x),'g')
+m,K = Infer.NormalFromMCMC(conv/thin,n_chains=no_ch)
+E3,par,par_err = Infer.ImportanceSamp(LogLikelihood_iid_mf,(Model3,x,y),m,K,chain_len)
+E3,par,par_err = Infer.ImportanceSamp(LogLikelihood_iid_mf,(Model3,x,y),m,3.*K,chain_len)
+
+print E1, E2, E3
+Infer.BayesFactor(E1,E2)
+Infer.BayesFactor(E1,E3,H1='H1',H2='H3')
+Infer.BayesFactor(E2,E3,H1='H2',H2='H3')
+
+raw_input()

examples/BlockedGibbs_test1.py

@@ -0,0 +1,70 @@
+#!/usr/bin/env python
+
+import numpy as np
+import pylab
+import os
+
+import Infer
+import MyMCMC
+from MyFunctions import Transit_aRs
+from MyFunctions import LogLikelihood_iid_mf
+
+#light curve parameters
+lc_pars = [.0,2.5,11.,.1,0.6,0.2,0.3,1.,0.]
+wn = 0.0003
+
+#create the data set (ie training data)
+time = np.arange(-0.1,0.1,0.001)
+flux = Transit_aRs(lc_pars,time) + np.random.normal(0,wn,time.size)
+
+#guess parameter values and guess uncertainties
+guess_pars = lc_pars + [wn]
+err_pars = [0.00001,0,0.2,0.0003,0.02,0.0,0.0,0.001,0.0001,0.0001]
+
+#plot the light curve + guess function
+pylab.figure(1)
+pylab.errorbar(time,flux,yerr=wn,fmt='.')
+pylab.plot(time,Transit_aRs(guess_pars[:-1],time),'r--')
+
+#first optimise the function
+guess_pars = Infer.Optimise(LogLikelihood_iid_mf,guess_pars[:],(Transit_aRs,time,flux),fixed=(np.array(err_pars) == 0)*1)
+
+#run a standard MCMC
+chain_len = 40000
+conv = 10000
+thin = 10
+no_ch=2
+adapt_lims = (2000,conv,5)
+glob_lims = (2000,conv,5)
+glob_lims = (0,0,0)
+Infer.MCMC(LogLikelihood_iid_mf,guess_pars[:],(Transit_aRs,time,flux),chain_len,err_pars,n_chains=no_ch,adapt_limits=adapt_lims,glob_limits=glob_lims,thin=thin)
+#Get parameters values/errors from chains
+par,par_err = Infer.AnalyseChains(conv/thin,n_chains=no_ch)
+bf_par = Infer.GetBestFit(n_chains=no_ch)
+print "Best Fit log p =", LogLikelihood_iid_mf(bf_par,Transit_aRs,time,flux)
+pylab.figure(2)
+Infer.PlotCorrelations(conv/thin,n_chains=no_ch,p=np.where(np.array(par_err)>0.)[0])
+
+#run a Blocked Gibbs MCMC
+gibbs_in = [1,0,2,2,2,0,0,1,1,1]
+no_steps = max(gibbs_in)
+chain_len = 20000
+conv = 10000/2
+thin = 10
+no_ch=2
+adapt_lims = (1000,conv,5)
+glob_lims = (1000,conv,5)
+glob_lims = (0,0,0)
+Infer.BGMCMC(LogLikelihood_iid_mf,guess_pars[:],(Transit_aRs,time,flux),chain_len,err_pars,gibbs_in,n_chains=no_ch,adapt_limits=adapt_lims,glob_limits=glob_lims,thin=thin)
+#Get parameters values/errors from chains
+par,par_err = Infer.AnalyseChains(conv/thin,n_chains=no_ch)
+bf_par = Infer.GetBestFit(n_chains=no_ch)
+print "Best Fit log p =", LogLikelihood_iid_mf(bf_par,Transit_aRs,time,flux)
+pylab.figure(3)
+Infer.PlotCorrelations(conv/no_steps/thin,n_chains=no_ch,p=np.where(np.array(par_err)>0.)[0])
+
+#plot fitted function
+pylab.figure(1)
+pylab.plot(time,Transit_aRs(par[:-1],time),'g-')
+
+raw_input()