In [ ]:

    

import healpy as hp


    

In [ ]:

    

import pymc as mc
import pandas as pd


    

In [ ]:

    

import patsy


    

In [ ]:

    

nside = 8
npix = hp.nside2npix(nside)


    

In [ ]:

    

vectors = hp.pix2vec(nside, np.arange(npix))


    

In [ ]:

    

true_monopole = 1.
true_dipole = np.array([.5, .5, 0])


    

In [ ]:

    

m = true_monopole + np.dot(true_dipole, vectors)
m += np.random.normal(0, .1, len(m))


    

In [ ]:

    

hp.mollview(m)


    

In [ ]:

    

hp.fit_dipole(m)


    

In [ ]:

    

labels = ["dx", "dy", "dz"]
data = pd.DataFrame({labels[i]:vectors[i] for i in range(3)})
data.index.name = "pix"
data["m"] = m


    

In [ ]:

    

outcome, predictors = patsy.dmatrices("m ~ dx + dy + dz + 1", data)
betas = np.linalg.lstsq(predictors, outcome)[0].ravel()
for name, beta in zip(predictors.design_info.column_names, betas):
      print("%s: %s" % (name, beta))


    

In [ ]:

    

import statsmodels.api as sm


    

In [ ]:

    

sm.OLS?


    

In [ ]:

    

sm.WLS?


    

In [ ]:

    

smresult = sm.WLS(outcome, predictors, weights=10).fit()
print smresult.summary()


    

In [ ]:

    

import statsmodels.api as sm


    

In [ ]:

    

import statsmodels.formula.api as smf


    

In [ ]:

    

data.head()


    

In [ ]:

    

results = smf.ols("m ~ dx + dy + dz + 1", data=data).fit()
print results.summary()


    

In [ ]:

    

with mc.Model() as model:
    mc.glm.glm("m ~ dx + dy + dz + 1", data)
    step = mc.NUTS() # Instantiate MCMC sampling algorithm
    trace = mc.sample(2000, step) # draw


    

In [ ]:

    

from pymc.plots import MultiTrace, kdeplot_op
vars = trace.varnames

if isinstance(trace, MultiTrace):
    trace = trace.combined()

n = len(vars)
f, ax = subplots(n, 2, squeeze=False, figsize=(15,15))

for i, v in enumerate(vars):
    d = np.squeeze(trace[v])

    if trace[v].dtype.kind == 'i':
        ax[i, 0].hist(d, bins=sqrt(d.size))
    else:
        kdeplot_op(ax[i, 0], d)
    ax[i, 0].set_title(str(v))
    ax[i, 1].plot(d, alpha=.35)

    ax[i, 0].set_ylabel("frequency")
    ax[i, 1].set_ylabel("sample value")

plt.tight_layout()


    

In [ ]: