The n-dimensional fitting tutorial¶
The funcFit package supports n-dimensional fitting, which means that both the domain and the range may be multidimensional.
An example for a multidimensional range would be the fitting of a circle or orbit, when time (1d) is mapped to image coordinates (2d). Fitting, for example, the structure of a point spread function (PSF) would be an example for a multidimensional domain, viz., the image coordinates, which are mapped onto a one-dimensional range (the flux or intensity).
Fitting a circular orbit¶
This example demonstrates how to fit a 2d model (location in plane) depending on a single variable (time), so that there is a mapping of the form
where the \(R^m\) denotes the parameter vector.
In particular, we assume x,y-position measurements at a number of times. Furthermore, all x and y measurements have an error.
import numpy as np
import matplotlib.pylab as plt
from PyAstronomy import funcFit as fuf
# Get the circular model and assign
# parameter values
c = fuf.Circle2d()
c["r"] = 1.0
c["t0"] = 0.0
c["per"] = 3.0
# Evaluate the model at a number of
# time stamps
t = np.linspace(0.0, 10.0, 20)
pos = c.evaluate(t)
# Add some error to the "measurement"
pos += np.reshape(np.random.normal(0.0, 0.2, pos.size), pos.shape)
err = np.reshape(np.ones(pos.size), pos.shape) * 0.2
# Define free parameters and fit the model
c.thaw(["r", "t0", "per"])
c.fit(t, pos, yerr=err)
c.parameterSummary()
# Evaluate the model at a larger number of
# points for plotting
tt = np.linspace(0.0, 10.0, 200)
model = c.evaluate(tt)
# Plot the result
plt.errorbar(pos[::, 0], pos[::, 1], yerr=err[::, 1],
xerr=err[::, 0], fmt='bp')
plt.plot(model[::, 0], model[::, 1], 'r--')
plt.show()
Note
The only difference between the “normal” one-dimensional case and that exemplified here, is the dimension of the y and error arrays.
Creating “coordinate arrays”¶
If one has multiple coordinate axes, such as in the case of an image, it is often convenient to represent the coordinates using an appropriate array. funcFit provides the coordinateGrid function to help constructing an appropriate array. In particular, for n coordinate axes, the function will construct an array with n+1 dimensions, which gives the physical coordinates for every array index.
Note
The user is free to use any format to give the coordinates. For example, the meshgrid in numpy could be applied to obtain a similar result. The Fitting image data with a 2d-Gaussian model example demonstrates both use cases.
from __future__ import print_function, division
from PyAstronomy import funcFit as fuf
import numpy as np
# Constructing the two individual coordinate axes
x = np.linspace(-2., 2., 50)
y = np.linspace(-2., 2., 50)
# Applying funcFit's "coordinateGrid" helper function
# to built appropriate array-index -> coordinate mapping
# needed for nD fitting.
g = fuf.coordinateGrid(x, y)
print("(x, y) coordinates at index (11, 28): ", g[11, 28])
Fitting image data with a 2d-Gaussian model¶
The following two examples demonstrate how to fit 2d image data with a 2d Gaussian model using funcFit, i.e., we have the mapping
where, again, \(R^m\) denotes the parameter vector, \(R^2\) is the image coordinate, and the result is the “level”, “flux”, or whatever the image shows.
In this case, we need to take care of the parameter representation. In principle, we may choose the coordinate format ad libitum, however, the evaluate method of the fitting objects has to understand the format. Below, we show two conceivable approaches to specify the parameters.
Let us say, we are given both an x- and a y-axis
specifying the coordinates in the image. As demonstrated in Creating “coordinate arrays”,
the function coordinateGrid() can be used to
construct an appropriate coordinate array mapping array index to
physical coordinates. The class GaussFit2d expects such a
coordinate array, and the example below demonstrates how to use it.
Note
Below the example, the evaluate method of GaussFit2d is shown to illustrate the use of the coordinate array.
from PyAstronomy import funcFit as fuf
import numpy as np
import matplotlib.pylab as plt
# Constructing the individual coordinate axes
x = np.linspace(-2., 2., 50)
y = np.linspace(-2., 2., 50)
# Applying funcFit's "coordinateGrid" helper function
# to built appropriate array-index -> coordinate mapping
# needed for nD fitting.
g = fuf.coordinateGrid(x, y)
# Create the 2d-Gaussian model and assign
# some model parameters.
gf = fuf.GaussFit2d()
gf["sigx"] = 0.75
gf["sigy"] = 0.4
gf["A"] = 1.0
gf["rho"] = 0.4
# Get the "data" by evaluating the model
# and adding some noise. Note that the coordinate
# mapping (array g) is passed to evaluate here.
im = gf.evaluate(g)
im += np.reshape(np.random.normal(0.0, 0.1, 2500), (50, 50))
err = np.ones((50, 50))*0.1
# Thaw parameters and fit
gf.thaw(["A", "rho"])
gf.fit(g, im, yerr=err)
# Show the resulting parameter values ...
gf.parameterSummary()
# ... and plot the result.
plt.title("Image data")
plt.imshow(np.transpose(im), origin="lower")
plt.show()
plt.title("Residuals")
plt.imshow(np.transpose(im - gf.evaluate(g)), origin="lower")
plt.show()
The evaluate method of the GaussFit2d class has the following from:
# The "evaluate" method of class GaussFit2d
# Note how the coordinate grid (co argument) is used in the evaluation
# of the model. Basically co[::,::,0] gives the x-coordinate for every
# image pixel. Likewise, co[::,::,1] defines the y-coordinate.
def evaluate(self, co):
"""
Evaluates the model for current parameter values.
Parameters
----------
co : array
Specifies the points at which to evaluate the model.
"""
if (self["sigx"] <= 0.0) or (self["sigy"] <= 0.0):
raise(PE.PyAValError("Width(s) of Gaussian must be larger than zero.", \
solution="Change width ('sigx/y')."))
if self["rho"] > 1.0:
raise(PE.PyAValError("The correlation coefficient must be 0 <= rho <= 1.", \
solution="Change width ('sigx/y')."))
result = 1.0/(2.*pi*self["sigx"]*self["sigy"]*sqrt(1.-self["rho"]**2)) * \
exp( ((co[::,::,0]-self["mux"])**2/self["sigx"]**2 + (co[::,::,1]-self["muy"])**2/self["sigy"]**2 - \
2.*self["rho"]*(co[::,::,0]-self["mux"])*(co[::,::,1]-self["muy"])/(self["sigx"]*self["sigy"])) / \
(-2.*(1.-self["rho"]**2)) )
return result
Alternatively, you may find it more hand the coordinate axes directly to the
evaluate function. This is possible and demonstrated in the following
example, which relies on the GaussFit2dTuple class, whose
evaluate method expects a tuple of coordinate axes.
Note the difference in the evaluate method. In this case, we
use numpy’s meshgrid function to create coordinate mapping similar
to the one used before.
from PyAstronomy import funcFit as fuf
import numpy as np
import matplotlib.pylab as plt
# Constructing the individual coordinate axes
x = np.linspace(-2., 2., 50)
y = np.linspace(-2., 2., 50)
# Create the 2d-Gaussian model and assign
# some model parameters.
gf = fuf.GaussFit2dTuple()
gf["sigx"] = 0.75
gf["sigy"] = 0.4
gf["A"] = 1.0
gf["rho"] = 0.4
# Get the "data" by evaluating the model
# and adding some noise. Note that the coordinate
# mapping (array g) is passed to evaluate here.
im = gf.evaluate((x, y))
im += np.reshape(np.random.normal(0.0, 0.1, 2500), (50, 50))
err = np.ones((50, 50))*0.1
# Thaw parameters and fit
gf.thaw(["A", "rho"])
gf.fit((x, y), im, yerr=err)
# Show the resulting parameter values ...
gf.parameterSummary()
# ... and plot the result.
plt.title("Image data")
plt.imshow(np.transpose(im), origin="lower")
plt.show()
plt.title("Residuals")
plt.imshow(np.transpose(im - gf.evaluate((x, y))), origin="lower")
plt.show()
The associated evaluate method looks like this:
# The "evaluate" method of class GaussFit2dTuple
# A coordinate grid is constructed using numpy's meshgrid
# function. However, you may do whatever you prefer to
# map the input to the model coordinates.
def evaluate(self, co):
"""
Evaluates the model for current parameter values.
Parameters
----------
co : array
Specifies the points at which to evaluate the model.
"""
if (self["sigx"] <= 0.0) or (self["sigy"] <= 0.0):
raise(PE.PyAValError("Width(s) of Gaussian must be larger than zero.", \
solution="Change width ('sigx/y')."))
if self["rho"] > 1.0:
raise(PE.PyAValError("The correlation coefficient must be 0 <= rho <= 1.", \
solution="Change width ('sigx/y')."))
xx, yy = meshgrid(co[0], co[1])
result = 1.0/(2.*pi*self["sigx"]*self["sigy"]*sqrt(1.-self["rho"]**2)) * \
exp( ((xx-self["mux"])**2/self["sigx"]**2 + (yy-self["muy"])**2/self["sigy"]**2 - \
2.*self["rho"]*(xx-self["mux"])*(yy-self["muy"])/(self["sigx"]*self["sigy"])) / \
(-2.*(1.-self["rho"]**2)) )
return result
