Let's see what's possible with astropy.models ...
import numpy as np
from astropy.models import models, fitting, parameters
# We use iminuit to check the results
# http://iminuit.github.io/iminuit/
import iminuit
# Let's use the example from here:
# http://nbviewer.ipython.org/5030045
# Data
xdata = np.array([0, 1, 2, 3, 4, 5])
ydata = np.array([1, 1, 5, 7, 8, 12])
sigma = np.array([1, 2, 1, 2, 1, 2], dtype=float)
# Model
def f(x, slope, intercept):
"""Linear function (two fit parameters)"""
return slope * x + intercept
# Compute parameters and errors with iminuit
def chi2(slope, intercept):
"""Define fit statistic, interpreting sigma as errors on ydata"""
chi = (ydata - f(xdata, slope, intercept)) / sigma
return np.sum(chi ** 2)
minuit = iminuit.Minuit(chi2, pedantic=False, print_level=0)
minuit.migrad()
minuit.hesse()
print 'values:', minuit.values
print 'errors:', minuit.errors
values: {'slope': 1.952830188677304, 'intercept': 0.7704402515790711} errors: {'slope': 0.3071475582636637, 'intercept': 0.8504530789683833}
# Compute parameters and errors with astropy.models
model = models.Poly1DModel(1)
# TODO: Is it possible to rename parameters to what I want?
#model.parnames = ['intercept', 'slope']
fitter = fitting.LinearLSQFitter(model)
fitter(xdata, ydata, weights=sigma ** (-1))
print dict(slope=model.c1, intercept=model.c0)
# TODO: Currently no way to compute parameter errors with astropy.models?
{'slope': [1.9528301886792441], 'intercept': [0.77044025157232865]}
import astropy.models.parameters as p
p._Parameter
astropy.models.parameters._Parameter
# Let's try implementing this as a user-defined model in astropy
from astropy.models.models import ParametricModel
from astropy.models.parameters import _Parameter
class Line(ParametricModel):
parnames = ['slope', 'intercept']
def __init__(self, slope, intercept, paramdim=1):
self.linear = True
self.ndim = 1
self.outdim = 1
self._slope = _Parameter(name='slope', val=slope, mclass=self, paramdim=paramdim)
self._intercept = _Parameter(name='intercept', val=intercept, mclass=self, paramdim=paramdim)
ParametricModel.__init__(self, self.parnames, paramdim=paramdim)
def eval(self, x, params):
return params[0] * x + params[1]
def __call__(self, x):
x, format = _convert_input(x, self.paramdim)
result = self.eval(x, self.psets)
return _convert_output(result, format)
model = Line(slope=1, intercept=1)
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) <ipython-input-7-20332935293d> in <module>() 19 return _convert_output(result, format) 20 ---> 21 model = Line(slope=1, intercept=1) <ipython-input-7-20332935293d> in __init__(self, slope, intercept, paramdim) 7 def __init__(self, slope, intercept, paramdim=1): 8 self.linear = True ----> 9 self.ndim = 1 10 self.outdim = 1 11 self._slope = _Parameter(name='slope', val=slope, mclass=self, paramdim=paramdim) AttributeError: can't set attribute
ERROR: AttributeError: can't set attribute [IPython.core.interactiveshell]
# Let's try implementing a this as a user-defined fit statistic in astropy
from astropy.models.fitting import Fitter
MAXITER = 100
EPS = 1e-10
class SLSQPFitter(Fitter):
def __init__(self, model, fixed=None, tied=None, bounds=None,
eqcons=None, ineqcons=None):
Fitter.__init__(self, model, fixed=fixed, tied=tied, bounds=bounds,
eqcons=eqcons, ineqcons=ineqcons)
if self.model.linear:
raise ModelLinearityException('Model is linear in parameters, '
'non-linear fitting methods should not be used.')
def errorfunc(self, fps, *args):
meas = args[0]
self.fitpars = fps
res = self.model(*args[1:]) - meas
return np.sum(res**2)
def __call__(self, x, y , maxiter=MAXITER, epsilon=EPS):
self.fitpars = optimize.fmin_slsqp(self.errorfunc, p0=self.model.parameters[:], args=(y, x),
bounds=self.constraints._bounds, eqcons=self.constraints.eqcons,
ieqcons=self.constraints.ineqcons)
model = models.Poly1DModel(1)
fitter = SLSQPFitter(model)
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-8-025b390ac3a8> in <module>() 26 27 model = models.Poly1DModel(1) ---> 28 fitter = SLSQPFitter(model) <ipython-input-8-025b390ac3a8> in __init__(self, model, fixed, tied, bounds, eqcons, ineqcons) 9 eqcons=None, ineqcons=None): 10 Fitter.__init__(self, model, fixed=fixed, tied=tied, bounds=bounds, ---> 11 eqcons=eqcons, ineqcons=ineqcons) 12 if self.model.linear: 13 raise ModelLinearityException('Model is linear in parameters, ' TypeError: __init__() got an unexpected keyword argument 'ineqcons'
ERROR: TypeError: __init__() got an unexpected keyword argument 'ineqcons' [IPython.core.interactiveshell]
Once the chi^2 examples work, let's see if we can do a Poisson likelihood fit (Cash statistic) with astropy.models
def cash(D, M, safe=False):
"""cash fit statistic where D = data, M = model
http://cxc.cfa.harvard.edu/sherpa/statistics/#cash
This is simply 2 x the negative log Poisson likelihood:
P = (M ** D) * exp(-M) / (D!)
log(P) = D * log(M) - M - log(D!)
cash = - 2 * log(P) (model-independent term dropped)
"""
stat = 2 * (M - D * log(M))
if safe:
return np.where(M > 0, stat, 1)
else:
return stat