Interpolation of a 1d tensor

Questions
1

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:
shift_representation
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!

2

There is a interpolate distribution class you can try. For an example see http://docs.pymc.io/notebooks/updating_priors.html

1 Like
3

Hello @junpenglao! Thanks a lot for your reply.

I have been reading the interpolated example you posted but I do not understand how I could apply it here:

My 1d tensor is not a distribution I want to sample but rather a curve (which I have shifted according to the prior) and now I want to trim according to the observed data so I can compare it with.

Is there something I am missing?

Thanks again!

4

I see. In that case, you can either wrap the interp1d as a theano function using the theano.tensor.as_op decorator, or try to implement a custom interpolation function.

2 Likes
5

Thanks again, I understand writing a custom theano ops is the way to go.