Friday, November 5, 2010

Simple Simulation

A simple demo for ADASS 2010 to fake data with Poisson noise and fit.

#!/usr/bin/env python

import sherpa.astro.ui as sherpa

# create simple grid [0.1, 10] with a bin size of 0.01
sherpa.dataspace1d(0.1,10,0.01, dstype=sherpa.Data1D)

# <demo> stop

# shorthand for normgauss1d model
gauss = sherpa.normgauss1d

# define reference model as sum of two normalized gaussians
sherpa.set_model(gauss.ref1+gauss.ref2)

# initialize parameter values
ref1.pos = 2
ref1.fwhm = 2.5

# define area under curve ref1
ref1.ampl = 75

ref2.pos = 8

ref2.fwhm = 5

# define area under curve ref2
ref2.ampl = 125

sherpa.plot_model()

# <demo> stop

# evaluate the model, add poisson noise, and populate the dataspace
sherpa.fake()

# plot faked data as 1D line plot with errorbars
sherpa.plot_data()

# <demo> stop

# fit the faked data with new instances of normgauss1d

sherpa.set_model(gauss.g1+gauss.g2)

# narrow the parameter space boundaries according to data
sherpa.guess(g1)

g1.pos = 2

sherpa.freeze(g1)

sherpa.guess(g2)
g2.pos = 8

sherpa.fit()

sherpa.plot_fit()

# <demo> stop

sherpa.thaw(g1)

sherpa.fit()

sherpa.plot_fit()

# <demo> stop

# Run confidence to estimate limits of best-fit parameter values
sherpa.conf()

# Overplot reference model to show goodness-of-fit
# The reference model is show in orange compared to fitted model in red.
sherpa.set_model(ref1+ref2)

sherpa.get_model_plot_prefs()['linecolor']='aquamarine'
sherpa.plot_model(overplot=True)

# <demo> stop

No comments:

Post a Comment