#!/usr/bin/env python
"""
A contrived toy model to see if I can get this thing to work

"""
import pylab
import numpy as np
import pymc
from pymc.Matplot import plot as mcplot

N = 1000
true_alpha = 2
noise = 0.001
# Simulated flux densities
s_arr = pymc.rpareto(true_alpha, 1, size=N)

alpha = pymc.Uniform('alpha', 1, 3)
flux = pymc.Pareto('flux', alpha=alpha, m=1, size=N)
flux_meas = pymc.Normal('flux_meas', mu=flux, tau=noise**-2, observed=True,
                             value=pymc.rnormal(s_arr, tau=noise**-2, size=N))



M = pymc.MCMC([alpha, flux, flux_meas])
M.sample(5000, 1000, 10)
mcplot(M)
pylab.show()