I am trying to use pymc3 to compute the distance to stars based on their spectrum (a histogram of their light).

Basically, you compare the observed spectrum with one of a similar star whose wavelength calibration is known:

In this basic model, there is only have one parameter, the redshift (z). My approach with pymc2 was interpolating at each iteration the model flux:

```
from os import environ
environ["MKL_THREADING_LAYER"] = "GNU"
import theano
import theano.tensor as tt
import pymc3 as pm
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate.interpolate import interp1d
# Model data
waveData = np.arange(401, 550, 1)
fluxData = np.array([1.303, 1.345, 1.3 , 1.352, 1.334, 1.33 , 1.294, 1.269, 1.18 ,
0.717, 1.069, 1.22 , 1.237, 1.298, 1.276, 1.283, 1.263, 1.257,
1.25 , 1.237, 1.279, 1.273, 1.263, 1.26 , 1.212, 1.21 , 1.208,
1.23 , 1.17 , 1.113, 1.125, 1.16 , 1.055, 0.596, 1.045, 1.17 ,
1.166, 1.196, 1.143, 1.197, 1.211, 1.205, 1.203, 1.182, 1.194,
1.181, 1.147, 1.158, 1.171, 1.19 , 1.21 , 1.178, 1.142, 1.163,
1.132, 1.182, 1.148, 1.144, 1.14 , 1.156, 1.166, 1.143, 1.162,
1.141, 1.136, 1.164, 1.143, 1.123, 1.151, 1.138, 1.113, 1.136,
1.108, 1.125, 1.16 , 1.13 , 1.128, 1.133, 1.109, 1.134, 1.12 ,
1.11 , 1.1 , 1.078, 1.002, 0.672, 0.982, 1.061, 1.075, 1.104,
1.08 , 1.039, 1.12 , 1.064, 1.093, 1.084, 1.098, 1.078, 1.072,
1.041, 1.055, 1.044, 1.064, 0.993, 1.026, 1.04 , 1.053, 1.008,
1.078, 1.017, 1.023, 1.036, 1.042, 1.023, 1.031, 1.042, 0.957,
1.041, 1.039, 1.022, 1.001, 1.023, 1.031, 1.04 , 1.022, 1.032,
0.963, 1.019, 1.033, 0.993, 1.023, 1.004, 0.972, 0.988, 1.013,
1.028, 0.977, 1.007, 1.002, 0.985, 0.974, 0.999, 0.955, 1.017,
0.999, 0.982, 0.994, 0.962, 0.981])
# Generating synthetic observation
z_True = 1.225
waveObs = np.arange(420, 500, 1) * z_True
fluxObs = interp1d(waveData, fluxData, bounds_error=True)(np.arange(420, 500, 1))
sigma_cont = 0.05
# Limits to the redshift
z_min_ssp = np.around((waveObs[-1] / waveData[-1]), decimals=2)
z_max_ssp = np.around((waveObs[0] / waveData[0]), decimals=2)
# Declare data tensors
waveObs_tt = theano.shared(waveData)
waveData_tt = theano.shared(waveData)
fluxData_tt = theano.shared(fluxData)
# Pymc3 model
with pm.Model() as model:
# Prior redshift
z = pm.Uniform('z', lower=z_min_ssp, upper=z_max_ssp)
# Proposal wavelength
wave_i = waveData_tt * z
# Interpolate model flux at observed wavelength range
flux_i = interp1d(wave_i, fluxData_tt, bounds_error=True)(waveObs_tt)
# Likelihood
Y = pm.Normal('Y', mu=flux_i, sd=sigma_cont, observed=fluxObs)
for RV in model.basic_RVs:
print(RV.name, RV.logp(model.test_point))
# Launch model
trace = pm.sample(5000, tune=1000)
```

This is not working, however, as interp1d is not tensor friendly.

I wonder if anyone would share how they do their interpolations or if there is an interpolation ops available.

Thanks!