with pm.Model() as gw_model:
neta = pm.Normal(‘neta’ , mu = 0.1 , sd = 0.2)
mass = pm.Normal(‘mass’ , mu = 25 , sd = 100)
amplitude = 10**-3
Mass = 30
distance = 1
Neta = 0.25
tc=0
phic =0
with gw_model :
phase = 2np.pifrequenciestc + phic + 3/(128 * neta * (np.pimass)(5/3))frequencies(-4/3)
exponent = 1jphase + 1j(np.pi/4)
strain = (amplitude * mass*(5/6))/(distancenp.pi**(2/3)) ((5*neta)/(24))**(-7/6)*np.exp(exponent)
diff = np.linalg.norm(waveform - strain)
likelihood = pm.Normal(‘obs’ , mu = diff
, sd =1 , observed = np.random.normal( 0 , 1, len(waveform)))
The model above is a very simplified model of a gravitational wave signal, I’m looking at using probabilistic programming for parameter estimation. The likelihood function at the most basic revolves around taking the norm of the difference between singals. The example above shows what I’m trying to do however when i evaluate the log probability at any point it returns a complex value.
i.e.
gw_model.logp({‘mass’: 100 , ‘neta’:0.1})
returns
array((-10766.848421086901-16136.450316121165j))
This (I think is leading to me being unable to evaluate the gradients so can’t use VI or sampling) ). I’m new to PyMC3 so may have made a mistake somewhere
the entire code can be found at
https://github.com/rgreen1995/Probabilistic-Programming/blob/master/PyMC3/PyMC_make_model.ipynb
Thanks