brent pymc test

gistfile1.py

import pymc as mc
import numpy as np
import pyfits as pf

arr = pf.getdata('stack3_gold.fits')
x=y=np.arange(0,71)
x,y=np.meshgrid(x,y)
err_map = pf.getdata('stack3wht_gold.fits')
def model((x,y),arr):
    amp = mc.Uniform('amp',lower=0,upper=np.amax(arr),doc='Amplitude')
    x0 = mc.Uniform('x0',lower=21,upper=51,doc='xo')
    y0 = mc.Uniform('y0',lower=21,upper=51,doc='yo')
    sigx = mc.Uniform('sigx',lower=0.1,upper=10,doc='Sigma in X')
    sigy = mc.Uniform('sigy',lower=0.1,upper=10,doc='Sigma in Y')
    theta = mc.Uniform('theta',lower=0,upper=2*np.pi,doc='Rotation')
    offset = mc.Uniform('c',lower=np.amin(arr),upper=np.amax(arr),doc='Vertical offset')
    
    @mc.deterministic(plot=False,trace=False)
    def gaussian((x, y), amplitude, xo, yo, sigma_x, sigma_y, theta, offset):
        xo = float(xo)
        yo = float(yo)    
        a = (mc.cos(theta)**2)/(2*sigma_x**2) + (mc.sin(theta)**2)/(2*sigma_y**2)
        b = -(mc.sin(2*theta))/(4*sigma_x**2) + (mc.sin(2*theta))/(4*sigma_y**2)
        c = (mc.sin(theta)**2)/(2*sigma_x**2) + (mc.cos(theta)**2)/(2*sigma_y**2)
        gauss = offset+amplitude*mc.exp(-1*(a*((x-xo)**2)+2*b*(x-xo)*(y-yo)+c*((y-yo)**2)))
        return np.ravel(gauss)
    flux = mc.Normal('flux',mu=gaussian,tau=err_map,value=arr,observed=True,doc='Observed Flux')
    return locals()
mdl = mc.MCMC(model((x,y),arr))
mdl.sample(iter=1e5,burn=9e4)