Hi,
I’m trying to use pymc3 to do voigt fitting. I tried using the scipy wofz function but it throws errors which I assume are because theano tensors don’t work with scipy. So I define a couple functions to calculate the profile. I get the following errors using the following pymc3 code (trimmed):
def faddeeva(z):
m = tt.exp(-tt.square(z))
return (m - m * tt.erf(1.j * z))
def voigt_profile(x, gamma):
z = (x + 1.j * gamma)
return faddeeva(z).real
def calculate_profile(lambda_0, f_osc, gamma, b, N, v0, lambda_bins=None):
c = speed_of_light_cgs.d[()]
lam1 = v0
nudop = 1e8 * b / lam1
tau_X = tau_factor * N * f_osc / b
tau0 = tau_X * lambda_0 * 1e-8
u = c / b * (lam1 / lambda_bins - 1.0)
a = gamma / (4.0 * np.pi * nudop)
phi = voigt_profile(u, a)
tauphi = tau0 * phi
return np.exp(-tauphi)
for ii in nvars:
model = mc.Model()
with model:
for ii in range(1, n + 1):
vars_dict[ii] = {}
vars_key.append(["v0_" + str(ii), "N_" + str(ii), "b_" + str(ii)])
vars_dict[ii]["v0"] = mc.Uniform(
"v0_" + str(ii), minima_vel[ii - 1] - 10, minima_vel[ii - 1] + 10
)
vars_dict[ii]["b"] = mc.Uniform("b_" + str(ii), 2, 40)
vars_dict[ii]["N"] = mc.Uniform("N_" + str(ii), 10, 17)
sd = 1.0 / mc.Uniform("sd", 0, 1) ** 2
for jj in range(1, n + 1):
vv0 = vars_dict[jj]["v0"]
nn = vars_dict[jj]["N"]
bb = vars_dict[jj]["b"]
flux_total *= mc.Deterministic(
"flux_obs" + str(jj), calculate_profile(lambda_0, f_osc, gamma, bb, nn, vv0, wav)
)
observation = mc.Normal("obs", mu=flux_total, tau=sd, observed=flux_temp)
with model:
trace = mc.sample(
5000,
tune=5000,
chains=2,
discard_tuned_samples=True,
target_accept=0.99,
return_inferencedata=True,
)
I get the following error :
/u/mrmgehlm/.local/lib/python3.7/site-packages/theano/tensor/elemwise.py:826: RuntimeWarning: overflow encountered in exp
variables = ufunc(*ufunc_args, **ufunc_kwargs)
Elemwise{mul,no_inplace}.0
/u/mrmgehlm/.local/lib/python3.7/site-packages/theano/tensor/elemwise.py:826: RuntimeWarning: overflow encountered in exp
variables = ufunc(*ufunc_args, **ufunc_kwargs)
Elemwise{mul,no_inplace}.0
ERROR (theano.graph.opt): SeqOptimizer apply <theano.tensor.opt.FusionOptimizer object at 0x1520f5743950>
ERROR (theano.graph.opt): Traceback:
ERROR (theano.graph.opt): Traceback (most recent call last):
File "/u/mrmgehlm/.local/lib/python3.7/site-packages/theano/scalar/basic.py", line 366, in filter
"Value cannot accurately be converted to dtype"
TypeError: Value cannot accurately be converted to dtype (complex128) and allow_downcast is not True
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "/u/mrmgehlm/.local/lib/python3.7/site-packages/theano/graph/opt.py", line 245, in apply
sub_prof = optimizer.optimize(fgraph)
File "/u/mrmgehlm/.local/lib/python3.7/site-packages/theano/graph/opt.py", line 83, in optimize
ret = self.apply(fgraph, *args, **kwargs)
File "/u/mrmgehlm/.local/lib/python3.7/site-packages/theano/tensor/opt.py", line 7727, in apply
new_outputs = self.optimizer(fgraph, node)
File "/u/mrmgehlm/.local/lib/python3.7/site-packages/theano/tensor/opt.py", line 7671, in local_fuse
ret = local_fuse(fgraph, new_node)
File "/u/mrmgehlm/.local/lib/python3.7/site-packages/theano/tensor/opt.py", line 7671, in local_fuse
ret = local_fuse(fgraph, new_node)
File "/u/mrmgehlm/.local/lib/python3.7/site-packages/theano/tensor/opt.py", line 7671, in local_fuse
ret = local_fuse(fgraph, new_node)
[Previous line repeated 7 more times]
File "/u/mrmgehlm/.local/lib/python3.7/site-packages/theano/tensor/opt.py", line 7544, in local_fuse
tmp.tag.test_value = tv.flatten()[0]
File "/u/mrmgehlm/.local/lib/python3.7/site-packages/theano/graph/utils.py", line 266, in __setattr__
obj = self.attr_filter(obj)
File "/u/mrmgehlm/.local/lib/python3.7/site-packages/theano/scalar/basic.py", line 371, in filter
f"Could not convert {type(data)} (value={data}) to {self.dtype}",
TypeError: Could not convert <class 'numpy.complex128'> (value=(nan+nanj)) to complex128
/u/mrmgehlm/.local/lib/python3.7/site-packages/theano/tensor/elemwise.py:826: RuntimeWarning: overflow encountered in exp
variables = ufunc(*ufunc_args, **ufunc_kwargs)
Traceback (most recent call last):
File "test_fit.py", line 370, in <module>
linefit()
File "/u/mrmgehlm/.local/lib/python3.7/site-packages/click/core.py", line 829, in __call__
return self.main(*args, **kwargs)
File "/u/mrmgehlm/.local/lib/python3.7/site-packages/click/core.py", line 782, in main
rv = self.invoke(ctx)
File "/u/mrmgehlm/.local/lib/python3.7/site-packages/click/core.py", line 1066, in invoke
return ctx.invoke(self.callback, **ctx.params)
File "/u/mrmgehlm/.local/lib/python3.7/site-packages/click/core.py", line 610, in invoke
return callback(*args, **kwargs)
File "test_fit.py", line 323, in linefit
return_inferencedata=True,
File "/u/mrmgehlm/.local/lib/python3.7/site-packages/pymc3/sampling.py", line 428, in sample
check_start_vals(model.test_point, model)
File "/u/mrmgehlm/.local/lib/python3.7/site-packages/pymc3/util.py", line 240, in check_start_vals
"Initial evaluation results:\n{}".format(elem, str(initial_eval))
pymc3.exceptions.SamplingError: Initial evaluation of model at starting point failed!
Starting values:
{'v0_1_interval__': array(0.), 'b_1_interval__': array(0.), 'N_1_interval__': array(0.), 'v0_2_interval__': array(0.), 'b_2_interval__': array(0.), 'N_2_interval__': array(0.), 'sd_interval__': array(0.)}
Initial evaluation results:
v0_1_interval__ -1.39
b_1_interval__ -1.39
N_1_interval__ -1.39
v0_2_interval__ -1.39
b_2_interval__ -1.39
N_2_interval__ -1.39
sd_interval__ -1.39
obs NaN
Name: Log-probability of test_point, dtype: float64
I’m not sure what’s going on here. This exact implementation works just fine if I use a simple gaussian instead of a voigt. Any ideas/suggestions are much appreciated!
Thanks
Aniket