User-defined Model and Fit Statistic

An advanced demo for ADASS 2010 illustrating a script that incorporates a user-defined model and a user-defined fit statistic with a prior into an X-ray spectral fit.

#!/usr/bin/env python
import sherpa.astro.ui as sherpa
from sherpa.models import Parameter, ArithmeticModel
from sherpa.stats import truncation_value
import numpy

class MyPowerLaw(ArithmeticModel):
    """
    Class to characterize a power-law function

    f(x) = ampl (x/ref) ^ (-gamma)

    """
    def __init__(self, name='mypowerlaw'):

        self.gamma = Parameter(name, "gamma", 1.0, min=-10, max=10)
        self.ref = Parameter(name, "ref", 1.0, alwaysfrozen=True)
        self.ampl = Parameter(name, "ampl", 1.0, min=0)

        ArithmeticModel.__init__(self, name,
                                 (self.gamma, self.ref, self.ampl))

    def calc(self, p, x, xhi=None, *args, **kwargs):
        """
        Params
        `p`   list of ordered parameter values.
        `x`   ndarray of domain values, bin midpoints or lower
              bin edges.
        `xhi` ndarray of upper bin edges.

        returns ndarray of calculated function values.
        """

        if xhi is None:
            return self.point(p, x)
        return self.integrated(p, x, xhi)


    @staticmethod
    def point(p, x):
        """
                                         / x_i  \(-p[0])
         point version, f_i(x_i) =  p[2] |------|
                                         \ p[1] /

        Params
        `p`  list of ordered parameter values.
        `x`  ndarray of bin midpoints.
        
        returns ndarray of calculated function values at
                bin midpoints
        """

        if xhi is None:
            return p[2]*numpy.power(x/p[1], -p[0])


    @staticmethod
    def integrated(p, xlo, xhi):
        """
        
         integrated form from lower bin edge to upper edge 
        
         ⌠ xhi_i              p[2]    /    /xhi \ (-p[0])     /xlo \ (-p[0]) \
         |      f_i(x_i) dx = ----- * |xhi |----|       - xlo |----|         |
         ⌡ xlo_i              1-p[0]  \    \p[1]/             \p[1]/         /
        
        Params
        `p`   list of ordered parameter values.
        `xlo` ndarray of lower bin edges.
        `xhi` ndarray of upper bin edges.

        returns ndarray of integrated function values over 
                lower and upper bin edges.
        """
        if p[0] == 1.0:
            return p[2] * p[1] * (numpy.log(xhi) - numpy.log(xlo))

        p1 = numpy.power(xlo, 1.0-p[0])/(1.0-p[0])
        p2 = numpy.power(xhi, 1.0-p[0])/(1.0-p[0])
        return  p[2] / numpy.power(p[1], -p[0]) * (p2 - p1)

# <demo> stop

def logNormal(x, mu, sigma):
    """
    Compute natural log of Normal distribution PDF analytically

    Params
    `x`      x
    `mu`     mean
    `sigma`  standard deviation

    returns ln(f_X(x; mu, sigma))
    """
    sigma_sqr = sigma*sigma
    return (-numpy.log(numpy.sqrt(2*numpy.pi*sigma_sqr)) -
             numpy.square(x-mu)/2/sigma_sqr)

# <demo> stop

# stat interface instancemethod
def calc_stat(self, data, model, staterror=None,
              syserror=None, weight=None):
    """
    Calculate the fit statistic value and an array of statistic
    value contributions per bin.

    Params
    `data`   ndarray of observed data points
    `model`  ndarray of predicted data points
    `staterror`  ndarray of statistical error on observed data points
    `syserror`   ndarray of systematic error on observed data points
    `weight`     ndarray of statistic weights

    returns a tuple of the fit statistic value and array of
            stat contributions per bin.
    """
    gamma = self.gamma.val             # nH (abs1.nh)
    beta  = self.beta.val              # powerlaw index (p1.gamma)
    xi    = numpy.log(self.alpha.val)  # log powerlaw norm (p1.ampl)

    phigamma = logNormal(gamma, self.mugamma.val, self.siggamma.val)
    phibeta  = logNormal(beta, self.mubeta.val, self.sigbeta.val)
    phixi    = logNormal(xi, self.muxi.val, self.sigxi.val)

    # truncation_value used as an approximation for
    # non-positive values in model
    if truncation_value > 0:
        model[model<=0.0] = truncation_value

    # Poisson log-likelihood summation
    likelihood = sum(-model + data*numpy.log(model))

    # Add the prior to log-likelihood
    prior = (phigamma+phibeta+phixi)
    stat = likelihood + prior

    # reverse sign to minimize log-likelihood!
    return (-stat, numpy.ones_like(data))


# <demo> stop

sherpa.load_pha("3c273.pi")
sherpa.notice(0, 6.0)

sherpa.add_model(MyPowerLaw)
sherpa.set_source(sherpa.xsphabs.abs1 * mypowerlaw.p1)

sherpa.load_user_stat("stat", calc_stat,
               priors=dict(mugamma=0.017, # hyperparameters
                           mubeta=1.26,
                           muxi=-8.5,
                           siggamma=1,
                           sigbeta=1,
                           sigxi=1,
                           gamma=abs1.nh, # fit parameters
                           beta=p1.gamma,
                           alpha=p1.ampl))

print stat

sherpa.set_stat(stat)
sherpa.set_method("neldermead")

sherpa.fit()

sherpa.plot_fit_delchi()

# <demo> stop

Estimate Flux Errors

A demo for ADASS 2010 to estimate the integrated flux from an X-ray spectral fit using simulations.

#!/usr/bin/env python
import sherpa.astro.ui as sherpa
import numpy

# load the data and fit with an absorbed power-law
sherpa.load_pha("3c273.pi")
sherpa.notice(0.5,7.0)
sherpa.subtract()

sherpa.set_stat("chi2datavar")
sherpa.set_model(sherpa.xsphabs.abs1*sherpa.xspowerlaw.p1)
abs1.nh=0.07
sherpa.freeze(abs1.nh)
sherpa.fit()

# <demo> stop

# simulate fluxes by sampling model parameters for 100 iterations.

flux100=sherpa.sample_energy_flux(0.5,7.,num=100)

# slice the flux column from simulation table.

fluxes = flux100[:,0]

print "Mean:", numpy.mean(fluxes)

print "Median:", numpy.median(fluxes)

print "Std Dev:", numpy.std(fluxes)

# determine 95% and 50% quantiles using numpy array sorting

sf=numpy.sort(fluxes)

print "95% quantile:", sf[0.95*(len(fluxes)-1)]
print "50% quantile:", sf[0.5*(len(fluxes)-1)]

# <demo> stop

# run 7.5e3 flux simulations and bin a histogram using the result

hist=sherpa.get_energy_flux_hist(0.5,7.,num=7500)
hist.histo_prefs['linecolor']=None
sherpa.plot_energy_flux(recalc=False)

# <demo> stop

# fit a Gaussian to the histogram
sherpa.load_arrays(2,hist.xlo,hist.xhi,hist.y,sherpa.Data1DInt)

sherpa.set_model(2,sherpa.gauss1d.g0)
g0.integrate=False
sherpa.guess(2,g0)
g0.ampl.min=0
g0.ampl.max=100
g0.ampl=1

sherpa.set_stat("leastsq")
sherpa.fit(2)

sherpa.plot_model(2, overplot=True)

# <demo> stop

# calculate sigma instead of FWHM
# FWHM = sigma * sqrt(8 ln 2)

fac0=numpy.sqrt(8*numpy.log(2))
sig=g0.fwhm.val/fac0

print "Sigma: ",sig,"Position: ", g0.pos.val

# write binned histogram arrays to file

sherpa.save_arrays("hist.dat",[hist.xlo,hist.xhi,hist.y],
                   fields=["xlo","xhi","y"], clobber=True)

# <demo> stop

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

Sherpa Mac Disk Image

I'm happy to announce a new deployment option for Sherpa. The latest version of Sherpa, 4.2.2, built as a Mac disk image for OSX 10.5 Leopard using native (meaning Non-CIAO) python.org 2.6 and NumPy installations. This experimental build is an effort to minimize the difficulties in building Sherpa's dependencies. This build includes XSPEC 12.6.0 models (why it's so large), CIAO region file support (including FITS), CIAO dynamic grouping, and WCS support using WCSSUBS. The user is still responsible for optional packages such as DS9, XPA, matplotlib, and pyFITS. The goal of this project is to utilize the Mac simple GUI installer and install Sherpa as easily as the SciPy and NumPy disk images.

We are currently working on a Snow Leopard/universal version with upgrades to Python 2.7 and NumPy 1.5.0.

http://cxc.harvard.edu/contrib/sherpa/downloads/sherpa-4.2.3-py2.6-python.org-macosx10.5-i386.dmg

New Sherpa Released!

New version of Sherpa was released today!

The standard CIAO version is available on Sherpa CIAO web page . The standalone version for Python users is available on Sherpa v2 page
Many enhancements include user interface to models including the support for combination of convolved and non-convolved model expressions, advanced user models, working with PSFs and table models.

Vinay's Cheat Sheet

Vinay has been writing a Sherpa Cheat Sheet on AstroStat Slog. This is a good list of basic "how to" tutorials for Sherpa.

User-defined Sherpa model types using Python Classes

Here is a quick example of a native Sherpa model (logparlaw1d) defined as a Python class. This class can then be registered as Sherpa model type using the add_model() function (found in Sherpa 4.2.2). The class attribute, name, can then be used as a model type within a standard Sherpa model expression.


To python users, instantiation of this class is obscured within the Sherpa model expression syntax. The add_model() function simply adds this class to Sherpa's function factory for creating model instances. The variable, mylogparlaw, in the local namespace is actually the Sherpa function factory. This variable's __getattr__() method has been overloaded to create instances of MyLogParLaw using the name of the attribute requested. Think of the following declarations as comparable.

lp1 = mylogparlaw.lp1

versus

lp1 = MyLogParLaw('lp1')

See the script below for more details

#!/usr/bin/env python

from sherpa.astro.ui import *
from sherpa.models import ArithmeticModel, Parameter
import numpy as np


class MyLogParaLaw(ArithmeticModel):
   """


               /  x \ - p[1] - p[2] * log10(x/p[0])
   f(x) = p[3] |----|
               \p[0]/

   """
   def __init__(self, name='mylogparlaw'):
       # p[0]
       self.ref = Parameter(name, 'ref', 1, alwaysfrozen=True)

       # p[1]
       self.c1 = Parameter(name, 'c1', 1)

       # p[2]
       self.c2 = Parameter(name, 'c2', 1)

       # p[3]
       self.ampl = Parameter(name, 'ampl', 1, 0)

       ArithmeticModel.__init__(self, name, (self.ref,self.c1,
                                             self.c2,self.ampl))


   def calc(self, p, xlo, xhi=None, *args, **kwargs):
       frac = (xlo/p[0])
       exponent = -p[1] - p[2] * np.log10(frac)
       return p[3] * np.power(frac, exponent)

if __name__ == '__main__':

   # register the class MyLogParLaw as a Sherpa model type
   add_model(MyLogParLaw)

   # Set the Sherpa model expression as the sum of MyLogParLaw
   # and Const1D
   set_model(mylogparlaw.lp1+const1d.c1)

   # Initialize the parameter values
   c1.c0 = 10
   lp1.ampl = 75
   lp1.c1 = 0.5
   lp1.c2 = 1.75
   lp1.ref = 2.1
   print lp1

   # Define a data space grid to be from [0.1, 10] in intervals
   # of 0.01
   dataspace1d(0.1,10,0.01, dstype=Data1D)

   # Evaluate the model over set of initialized parameters, add
   # poisson noise, and set as data space values.
   fake()

   # Use Chi2 with data variance as fit statistic
   set_stat('chi2datavar')

   # Set the statistical error of data space as 10% of faked
   # data values
   set_staterror(0.1, fractional=True)

   # Plot faked data space values with statistical errorbars
   plot_data()

The plot of faked data should look something like the image above.