Fitting

Currently multicomponent Gaussian function is available. This is done by creating a fitting object, setting up the fit and actually fitting the data. Fitting can either be done on a single scantable selection or on an entire scantable using the auto_fit function. If single value fitting is used, and the current selection includes multiple spectra (beams, IFs, scans etc) then the first spectrum in the scantable will be used for fitting.

 ASAP>f = fitter()
 ASAP>f.set_function(gauss=2) # Fit two Gaussians
 ASAP>f.set_scan(scans)
 ASAP>selection = selector()
 ASAP>selection.set_polarisations(1) # Fit the second polarisation
 ASAP>scans.set_selection(selection)
 ASAP>scans.set_unit('km/s')  # Make fit in velocity units
 ASAP>f.fit(1)                # Run the fit on the second row in the table
 ASAP>f.plot()                # Show fit in a plot window
 ASAP>f.get_parameters()      # Return the fit paramaters

This auto-guesses the initial values of the fit and works well for data without extra confusing features. Note that the fit is performed in whatever unit the abscissa is set to.

If you want to confine the fitting to a smaller range (e.g. to avoid band edge effects or RFI you must set a mask.

  ASAP>f = fitter()
  ASAP>f.set_function(gauss=2)
  ASAP>scans.set_unit('km/s')  # Set the mask in channel units
  ASAP>msk = s.create_mask([1800,2200])
  ASAP>scans.set_unit('km/s')  # Make fit in velocity units
  ASAP>f.set_scan(s,msk)
  ASAP>f.fit()
  ASAP>f.plot()
  ASAP>f.get_parameters()

If you wish, the initial parameter guesses can be specified and specific parameters can be fixed:

  ASAP>f = fitter()
  ASAP>f.set_function(gauss=2)
  ASAP>f.set_scan(s,msk)
  ASAP>f.fit() # Fit using auto-estimates
  # Set Peak, centre and fwhm for the second gaussian.
  # Force the centre to be fixed
  ASAP>f.set_gauss_parameters(0.4,450,150,0,1,0,component=1)
  ASAP>f.fit() # Re-run the fit

The fitter plot function has a number of options to either view the fit residuals or the individual components (by default it plots the sum of the model components).

Examples:

  # Plot the residual
  ASAP>f.plot(residual=True)

  # Plot the first 2 componentsa
  ASAP>f.plot(components=[0,1])

  # Plot the first and third component plus the model sum
  ASAP>f.plot(components=[-1,0,2])  # -1 means the compoment sum