Fitting beta_app and other stuff in pymc3

pymc_fit_G.py

import pymc3 as pm
import numpy as np


betas, betas_err, thetas = np.loadtxt("/home/ilya/github/lags_utils.py/data.txt",
                                      unpack=True)
thetas /= 180/np.pi


# Use beta_app and theta from Hovatta
with pm.Model() as model:
    G = pm.Uniform("G", 1, 100)
    beta = pm.Deterministic("beta", T.sqrt(G**2-1)/G)
    beta_ = pm.Deterministic("beta_", beta*T.sin(thetas)/(1-beta*T.cos(thetas)))
    likelihood = pm.Normal("beta_obs", mu=beta_, sd=betas_err, observed=betas)
    trace = pm.sample(njobs=4, tune=1000)
    ax1 = pm.traceplot(trace, varnames=[G])
    ax2 = pm.plot_posterior(trace, varnames=[G])


# Use only beta_app
with pm.Model() as model:
    G = pm.Uniform("G", 1, 100)
    beta = pm.Deterministic("beta", T.sqrt(G**2-1)/G)
    thetas = pm.Uniform("thetas", 0, np.pi/3, shape=len(betas))
    beta_ = pm.Deterministic("beta_", beta*T.sin(thetas)/(1-beta*T.cos(thetas)))
    likelihood = pm.Normal("beta_obs", mu=beta_, sd=betas_err, observed=betas)
    trace = pm.sample(njobs=4, tune=1000)
    ax1 = pm.traceplot(trace, varnames=["G", "thetas"])
    ax2 = pm.plot_posterior(trace, varnames=["G", "thetas"])
    samples = np.hstack((trace.get_values("G").reshape(2000, 1),
                         trace.get_values("thetas")))
    fig_corner = corner.corner(samples[::2], labels=['G']+["theta_{}".format(i) for
                                                        i in range(len(betas_err))],
                               show_titles=True,
                               quantiles=[0.16, 0.50, 0.84])