Hello everyone! I am trying to update my code from pymc2 to pymc3 and I wonder if anyone had the same type of issues (or actually, if it can be done). I have this observational data (a histogram with the light from a galaxy, the stellar component is a flat continuum while the gas produces narrow lines):
As annotated in the plot I want to apply two operations: First, fit the continuum with a physical model and then fit the area of the narrow line with another physical model.
To do the second step in a pymc3 model I need to remove the fitted continuum before computed the observed area of the line. My strategy was using the ā.tag.test_valueā attribute to do this. I am not getting errors with this dummy example:
from os import environ
environ["MKL_THREADING_LAYER"] = "GNU"
import theano.tensor as tt
import pymc3 as pm
import numpy as np
from scipy.integrate import simps
import matplotlib.pyplot as plt
def planck_radiation(wav, T, h, c, k):
a = 2.0*h*c**2
b = h*c/(wav*k*T)
intensity = a / ((np.power(wav,5)) * (np.exp(b) - 1.0))
return intensity
def planck_radiation_tt(wav, T, h, c, k):
a = 2.0*h*c**2
b = h*c/(wav*k*T)
intensity = a / ((tt.power(wav,5)) * (tt.exp(b) - 1.0))
return intensity
def gaussian(wave, A, mu, sig):
return A * np.exp(-np.power(wave - mu, 2.) / (2 * np.power(sig, 2.)))
# Physical constants
h = 6.626e-34
c = 3.0e+8
k = 1.38e-23
# Observed wavelength range
obs_wave = np.arange(400,500)
# Compute stellar continuum
T_true = 10000.0
continuum_true = planck_radiation(obs_wave * 1e-9, T_true, h, c, k)
# Calculate line true shape and area
A_true, mu_true, sig_true = 1e15, 460, 1
sqrt2pi = np.sqrt((2 * np.pi))
emline_true = gaussian(obs_wave, A_true, mu_true, sig_true)
emLineArea_true = A_true*sig_true*sqrt2pi
emLineArea_err = emLineArea_true * 0.05
# Mask for the line region
idcs_mask = (obs_wave > mu_true-5) & (obs_wave < mu_true+5)
mask = np.ones(obs_wave.size)
mask[idcs_mask] = 0
# Observed flux with noise
err_continuum = np.mean(continuum_true) * 0.05
obs_flux = continuum_true + emline_true + np.random.normal(0, err_continuum, obs_wave.size)
# Pymc3 model
with pm.Model() as model:
temp = pm.Normal('temp', mu=5000.0, sd=1000.0)
A_norm = pm.HalfNormal('A_norm', sd=5)
sig = pm.HalfNormal('sigma', sd=5)
# Model continuum
continuum_flux = planck_radiation_tt(obs_wave * 1e-9, temp, h, c, k) * mask
# Likelihood model continuum (masking the line)
Y_continuum = pm.Normal('Y_continuum', mu=continuum_flux, sd=err_continuum, observed=obs_flux*mask)
# Remove from the observation the proposed continuum
emission_obs = obs_flux - continuum_flux.tag.test_value
# Integrate the line area
line_area_obs = simps(emission_obs[idcs_mask], obs_wave[idcs_mask]) / 1e15
# Model line area
line_area = A_norm*sig*sqrt2pi
# Global multivariable likelihood for all lines
Y_line = pm.Normal('Y_line', mu=line_area, sd=emLineArea_err, observed=line_area_obs)
for RV in model.basic_RVs:
print(RV.name, RV.logp(model.test_point))
# Launch model
trace = pm.sample(8000, tune=2000)
# Output trace data
print pm.summary(trace)
pm.traceplot(trace)
plt.show()
But while the continuum fitting gives a good result line area priors are not being sampled:
I wonder if anyone would share with me their approach in those cases where the observed input data needs to be treated/corrected with the model variables.
Thanks for any advice.
:
