Merge branch 'master', remote branch 'origin' into gradientBranch

Author: jsalvatier

docs/jss/jss_article_static.tex

@@ -42,7 +42,7 @@ def manual2article(text):
 
     # Proglang and pkg
     text = re.sub(r"\\pdfbookmark.*", r"", text)
-    for pkgname in ['PyMC','NumPy','SciPy','PyTables','Matplotlib','Pylab','Pyrex']:
+    for pkgname in ['PyMC','NumPy','SciPy','PyTables','Matplotlib','Pylab','Pyrex','WinBUGS','JAGS','Hierarchical Bayes Compiler','pyTables','pydot','IPython','nose','gfortran','gcc','Enthought Python Distribution','Gnu Compiler Collection','GCC','g77','subversion']:
         text = re.sub(pkgname, r'\\pkg{%s}'%pkgname,text)
     for proglang in ['Python','Fortran']:
         text = re.sub(' ' + proglang, r' \\proglang{%s}'%proglang, text)

docs/jss/pymc.bib

@@ -267,7 +267,7 @@ \section{Variables' values and log-probabilities}
 
 \section{Fitting the model with MCMC}
 
-PyMC provides several objects that fit probability models (linked collections of variables) like ours. The primary such object, \code{MCMC}, fits models with the Markov chain Monte Carlo algorithm \cite{Gamerman:1997tb}. To create an \code{MCMC} object to handle our model, import \module{DisasterModel.py} and use it as an argument for \code{MCMC}:
+PyMC provides several objects that fit probability models (linked collections of variables) like ours. The primary such object, \code{MCMC}, fits models with a Markov chain Monte Carlo algorithm \cite{Gamerman:1997tb}. To create an \code{MCMC} object to handle our model, import \module{DisasterModel.py} and use it as an argument for \code{MCMC}:
 \begin{verbatim}
    >>> from pymc.examples import DisasterModel
    >>> from pymc import MCMC

docs/manual2article.py

@@ -604,7 +604,11 @@ def remember(self, chain=-1, trace_index = None):
 
         for variable in self._variables_to_tally:
             if isinstance(variable, Stochastic):
-                variable.value = self.trace(variable.__name__, chain=chain)[trace_index]
+                try:
+                    variable.value = self.trace(variable.__name__, chain=chain)[trace_index]
+                except:
+                    cls, inst, tb = sys.exc_info()
+                    warnings.warn('Unable to remember value of variable %s. Original error: \n\n%s: %s'%(variable,cls.__name__,inst.message))
 
     def trace(self, name, chain=-1):
         """Return the trace of a tallyable object stored in the database.

docs/tutorial.tex

@@ -1146,7 +1146,7 @@ def gamma_expval(alpha, beta):
 
     Expected value of gamma distribution.
     """
-    return asarray(alpha) / beta
+    return 1. * asarray(alpha) / beta
 
 def gamma_like(x, alpha, beta):
     R"""
@@ -1410,7 +1410,7 @@ def inverse_gamma_expval(alpha, beta):
 
     Expected value of inverse gamma distribution.
     """
-    return 1. / (asarray(beta) * (alpha-1.))
+    return 1. * asarray(beta) / (alpha-1.)
 
 def inverse_gamma_like(x, alpha, beta):
     R"""

pymc/Model.py

@@ -7,8 +7,7 @@
 # want to change the 'm' parameters
 # below.
 
-# The MCMC takes several hours on my machine. To make it run faster, thin the
-# dataset in getdata.py
+# The MCMC takes several hours on my machine. To make it run faster, thin the dataset.
 
 
 
@@ -19,35 +18,41 @@
 from pylab import *
 from numpy import *
 import model
+import pymc as pm
 
 
 data = csv2rec('duffy-jittered.csv')
-DuffySampler = MCMC(model.make_model(data), db='hdf5', dbcomplevel=1, dbcomplib='zlib')
+# Convert from record array to dictionary
+data = dict([(k,data[k]) for k in data.dtype.names])
+# Create model. Use the HDF5 database backend, because the trace will take up a lot of memory.
+DuffySampler = pm.MCMC(model.make_model(**data), db='hdf5', dbcomplevel=1, dbcomplib='zlib', dbname='duffydb.hdf5')
 
 # ====================
 # = Do the inference =
 # ====================
-# Use the HDF5 database backend, because the trace will take up a lot of memory.
-DuffySampler.use_step_method(GPEvaluationGibbs, walker_v, V, d)
+# Use GPEvaluationGibbs step methods.
+DuffySampler.use_step_method(pm.gp.GPEvaluationGibbs, DuffySampler.sp_sub_b, DuffySampler.V_b, DuffySampler.tilde_fb)
+DuffySampler.use_step_method(pm.gp.GPEvaluationGibbs, DuffySampler.sp_sub_s, DuffySampler.V_s, DuffySampler.tilde_fs)
+# Run the MCMC.
 DuffySampler.isample(50000,10000,100)
 
-n = len(DuffySampler.trace('V')[:])
-
+n = len(DuffySampler.trace('V_b')[:])
 
 # ==========================
 # = Mean and variance maps =
 # ==========================
 
 import tables
 covariate_raster = tables.openFile('africa.hdf5')
-xplot = covariate_raster.lon[:]*pi/180.
-yplot = covariate_raster.lat[:]*pi/180.
+xplot = covariate_raster.root.lon[:]*pi/180.
+yplot = covariate_raster.root.lat[:]*pi/180.
+
+data = ma.masked_array(covariate_raster.root.data[:], mask = covariate_raster.root.mask[:])
 
-data = ma.masked_array(covariate_raster.data[:], mask = covariate_raster.mask[:])
 where_unmasked = np.where(True-data.mask)
-dpred = dstack(meshgrid(xplot,yplot))[where_unmasked]
+dpred = dstack(meshgrid(xplot,yplot))[::-1][where_unmasked]
+africa = covariate_raster.root.data[:][where_unmasked]
 
-# This computation is O(m^2)
 Msurf = zeros(data.shape)
 E2surf = zeros(data.shape)
 
@@ -56,64 +61,46 @@
     # Reset all variables to their values at frame i of the trace
     DuffySampler.remember(0,i)
     # Evaluate the observed mean
-    Msurf_b, Vsurf_b = point_eval(DuffySampler.sp_sub_b.M_obs.value, DuffySampler.sp_sub_b.C_obs.value, dpred)
-    Msurf_0, Vsurf_0 = point_eval(DuffySampler.sp_sub_0.M_obs.value, DuffySampler.sp_sub_0.C_obs.value, dpred)
+    store_africa_val(DuffySampler.sp_sub_b.M_obs.value, dpred, africa)
+    Msurf_b, Vsurf_b = pm.gp.point_eval(DuffySampler.sp_sub_b.M_obs.value, DuffySampler.sp_sub_b.C_obs.value, dpred)
+    Msurf_s, Vsurf_s = pm.gp.point_eval(DuffySampler.sp_sub_s.M_obs.value, DuffySampler.sp_sub_s.C_obs.value, dpred)
+    Vsurf_b += DuffySampler.V_b.value
+    Vsurf_s += DuffySampler.V_s.value
 
-    freq_b = pm.invlogit(pm.rnormal(Msurf_b, 1/Vsurf_b))
-    freq_0 = pm.invlogit(pm.rnormal(Msurf_0, 1/Vsurf_0))
+    freq_b = pm.invlogit(Msurf_b +pm.rnormal(0,1)*np.sqrt(Vsurf_b))
+    freq_s = pm.invlogit(Msurf_s +pm.rnormal(0,1)*np.sqrt(Vsurf_s))
 
-    samp_i = (freq_b*freq_0+(1-freq_b)*DuffySampler.p1.value)**2
+    samp_i = (freq_b*freq_s+(1-freq_b)*DuffySampler.p1.value)**2
 
-    Msurf[where_unmasked] += samp_i/n
+    Msurf[where_unmasked] += samp_i/float(n)
     # Evaluate the observed covariance with one argument
-    E2surf[where_unmasked] += samp_i**2/n
+    E2surf[where_unmasked] += samp_i**2/float(n)
 
 # Get the posterior variance and standard deviation
 Vsurf = E2surf - Msurf**2
 SDsurf = sqrt(Vsurf)
 
-
+Msurf = ma.masked_array(Msurf, mask=covariate_raster.root.mask[:])
+SDsurf = ma.masked_array(SDsurf, mask=covariate_raster.root.mask[:])
+covariate_raster.close()
 
 # Plot mean and standard deviation surfaces
 close('all')
-imshow(Msurf, extent=[x.min(),x.max(),y.min(),y.max()],interpolation='nearest')
-plot(x,y,'r.',markersize=4)
-axis([x.min(),x.max(),y.min(),y.max()])
+imshow(Msurf[::-1,:], extent=np.array([xplot.min(),xplot.max(),yplot.min(),yplot.max()])*180./pi,interpolation='nearest')
+plot(DuffySampler.lon,DuffySampler.lat,'r.',markersize=4)
+axis(np.array([xplot.min(),xplot.max(),yplot.min(),yplot.max()])*180./pi)
 title('Posterior predictive mean surface')
+xlabel('Degrees longitude')
+ylabel('Degrees latitude')
 colorbar()
 savefig('duffymean.pdf')
 
 figure()
-imshow(SDsurf, extent=[x.min(),x.max(),y.min(),y.max()],interpolation='nearest')
-plot(x,y,'r.',markersize=4)
-axis([x.min(),x.max(),y.min(),y.max()])
+imshow(SDsurf[::-1,:], extent=np.array([xplot.min(),xplot.max(),yplot.min(),yplot.max()])*180./pi,interpolation='nearest')
+plot(DuffySampler.lon,DuffySampler.lat,'r.',markersize=4)
+axis(np.array([xplot.min(),xplot.max(),yplot.min(),yplot.max()])*180./pi)
 title('Posterior predictive standard deviation surface')
+xlabel('Degrees longitude')
+ylabel('Degrees latitude')
 colorbar()
-savefig('duffyvar.pdf')
-
-
-# # ====================
-# # = Realization maps =
-# # ====================
-# 
-# # Use thinner input arrays, this computation is O(m^6)!!
-# m = 101
-# xplot = linspace(x.min(),x.max(),m)
-# yplot = linspace(y.min(),y.max(),m)
-# dpred = dstack(meshgrid(yplot,xplot))
-# 
-# indices = random_integers(n,size=2)
-# for j,i in enumerate(indices):
-#     # Reset all variables to their values at frame i of the trace
-#     DuffySampler.remember(0,i)
-#     # Evaluate the Gaussian process realisation
-#     R = DuffySampler.walker_v.f.value(dpred)
-# 
-#     # Plot the realization
-#     figure()
-#     imshow(R,extent=[x.min(),x.max(),y.min(),y.max()],interpolation='nearest')
-#     plot(x,y,'r.',markersize=4)
-#     axis([x.min(),x.max(),y.min(),y.max()])
-#     title('Realization from the posterior predictive distribution')
-#     colorbar()
-#     savefig('elevdraw%i.pdf'%j)
+savefig('duffyvar.pdf')