Firstly, let me thank Chris, Colin and Eric for their help during the Cleveland pycon. Also thanks to Junpeng Lao for this very detailed solution which I have tried to apply to my problem:
I work in astronomy and most of our fittings depend on external data grids. I am trying to update some of my indexing and interpolating operations from pymc2 to pymc3 but I am having troubles with it.
Let’s hope this example helps: let’s say I want to use the ideal gas law (PV=nRT) to measure the specific gas constant (Rs) from some local temperature measurements. Under these conditions the equation can be rearranged to T = P / (D * Rs), where (T, P, D) are the air temperature, pressure and density. Le’ts assume the pressure remains pretty constant and that the local density magnitude can be calculated from a 2 dimensional grid which depends on two variables (x and y).
Looking at this question solution I tried to do a marginalized approach:
import os
os.environ[“MKL_THREADING_LAYER”] = “GNU”
import numpy as np
import theano
import theano.tensor as tt
import pymc3 as pm
import matplotlib.pyplot as pltConstants
R_air_true = 286.986 # J / (kg * K)
#Street measurements
P_true = 101325 # kg * m * s**-2
n_obs = 10
T_obs = np.random.uniform(low=10, high=12, size=n_obs) + 273
T_sigma = 1 # KModel ingredients data:
axis_length = 10
X_coords = np.linspace(5, 40, axis_length)
Y_coords = np.linspace(10, 20, axis_length)
density_grid = np.linspace(0.5, 2, 100).reshape((axis_length, axis_length))Converting to theano variables
X_coords_t = theano.shared(X_coords)
Y_coords_t = theano.shared(Y_coords)
density_grid_t = theano.shared(density_grid)Control the sparsity?
Dirprior = 1./axis_length
Model using street measurements
with pm.Model() as outside_temperature:
# Priors for the air gas constant R_air = pm.Normal('R_air', mu=300, sd=200) # Priors for the grid indeces beta1_x = pm.Gamma('beta1_x', alpha=Dirprior, beta=1., shape=axis_length) beta2_y = pm.Gamma('beta2_y', alpha=Dirprior, beta=1., shape=axis_length) # Variables for marginalized vectors x_idx_vector = pm.Deterministic('x', beta1_x/tt.sum(beta1_x)) y_idx_vector = pm.Deterministic('y', beta2_y/tt.sum(beta2_y)) # Values for X and Y coordinates (not needed in this example) x = x_idx_vector.dot(X_coords_t) y = y_idx_vector.dot(Y_coords_t) # Computing interpolated value from grid grid_column = density_grid_t.dot(x_idx_vector) rho_air = grid_column.dot(y_idx_vector) # Simulated expectation value T_synth = P_true / (R_air * rho_air) # Declare the model likelihood: Y_obs = pm.Normal('Y_obs', mu=T_synth, sd=T_sigma, observed=T_obs) # Run NUTS sampler trace = pm.sample(4000, tune=2000)Display summary
print pm.summary(trace)
Display trace
pm.traceplot(trace)
plt.show()
However, I am not getting a good solution and the sampling seems very slow:
Then tried generate two priors for the x and y variables indexes from which I do three interpolations: 2 for the x and y intervals and one bilinear interpolation for the grid:
import os
import numpy as np
os.environ[“MKL_THREADING_LAYER”] = “GNU”
import theano
import theano.tensor as tt
import pymc3 as pm
import matplotlib.pyplot as pltConstants
R_air_true = 286.986 # J / (kg * K)
#Street measurements
P_true = 101325 # kg * m * s**-2
n_obs = 10
T_obs = np.random.uniform(low=10, high=12, size=n_obs) + 273
T_sigma = 1 # KModel ingredients data:
axis_length = 10
X_coords = np.linspace(5, 40, axis_length)
Y_coords = np.linspace(10, 20, axis_length)
density_grid = np.linspace(0.5, 2, 100).reshape((axis_length, axis_length))Converting to theano variables
X_coords_t = theano.shared(X_coords)
Y_coords_t = theano.shared(Y_coords)
density_grid_t = theano.shared(density_grid)Model using street measurements
with pm.Model() as outside_temperature:
# Priors for the air gas constant R_air = pm.Normal('R_air', mu=300, sd=100) # Priors for the grid indeces x_idx = pm.Uniform('x_idx0', lower=0, upper=axis_length) y_idx = pm.Uniform('y_idx0', lower=0, upper=axis_length) # Use floor operaton to get the closest points in the grid x_idx1, y_idx1 = tt.cast(pm.math.floor(x_idx), 'int64'), tt.cast(pm.math.floor(y_idx), 'int64') x_idx2, y_idx2 = x_idx1 + 1, y_idx1 + 1 # Interpolate the values for X and Y coordinates x_in = X_coords_t[x_idx1] + (x_idx - x_idx1) * (X_coords_t[x_idx2] - X_coords_t[x_idx1])/(x_idx2 - x_idx1) y_in = Y_coords_t[y_idx1] + (y_idx - y_idx1) * (Y_coords_t[y_idx2] - Y_coords_t[y_idx1])/(y_idx2 - y_idx1) # Neighbouring points for the grid for X and Y coordinates: x1, x2 = X_coords_t[x_idx1], X_coords_t[x_idx2] y1, y2 = Y_coords_t[y_idx1], Y_coords_t[y_idx2] z11, z12 = density_grid_t[x_idx1, y_idx1], density_grid_t[x_idx1, y_idx2] z21, z22 = density_grid_t[x_idx2, y_idx1], density_grid_t[x_idx2, y_idx2] #Bilinear interpolation of the grid rho_air = (z11 * (x2 - x_in) * (y2 - y_in) + z21 * (x_in - x1) * (y2 - y_in) + z12 * (x2 - x_in) * (y_in - y1) + z22 * (x_in - x1) * (y_in - y1)) / ((x2 - x1) * (y2 - y1)) # Simulated expectation value T_synth = P_true / (R_air * rho_air) # Declare the model likelihood: Y_obs = pm.Normal('Y_obs', mu=T_synth, sd=T_sigma, observed=T_obs) # Run NUTS sampler trace = pm.sample(4000, tune=2000)Display summary
print pm.summary(trace)
Display trace
pm.traceplot(trace)
plt.show()
However, this code crashes very fast with this error:
Backtrace when the node is created(use Theano flag traceback.limit=N to make it longer):
File “/home/vital/PycharmProjects/thesis_pipeline/spectrum_fitting/pymc3_testing.py”, line 211, in
y1, y2 = Y_coords_t[y_idx1], Y_coords_t[y_idx2]HINT: Use the Theano flag ‘exception_verbosity=high’ for a debugprint and storage map footprint of this apply node.
I wonder if anyone would please give me any advice on the right way to explore external data grids in a pymc3 model so I can get started in the right track. Also, I have read one of the Junpeng Lao solutions advising to use the Interpolated distribution. I have not found any example of this distribution: How could it be applied here?
Thanks