<jupyter_start><jupyter_text>Fitting routinesFor all methods, the underlying package LMFIT handles the parabolic error calculations and passes along the optimization assignment to SciPy. Chisquare methodWhen fitting with the chisquare method, the costfunction to be minimized is$$\chi^2 = \left(\frac{y_i-HFS(x_i)}{\sigma_i}\right)^2$$with the subscript *i* referring to the datapoint, and *HFS* to the response of the :class:`.HFSModel` class. The algorithm used is the Levenberg-Marquardt algorithm, which gives quick results in a reliable fashion. The :func:`.fitting.chisquare_model` function creates the costfunction used for counting data. The function :func:`.fitting.chisquare_fit` performs the actual fit, while :func:`.fitting.chisquare_spectroscopic_fit` calculates the uncertainty on the data by taking the square root of the number of counts. One of the options for :func:`.fitting.chisquare_model` is the argument *func*. If the square root function is given, the uncertainty is calculated by taking the square root of the *fitvalue*, not the datapoint. This approach mimics the use of the Poisson distribution instead of the Gaussian distribution for the uncertainty calculation. In order to demonstrate the fit functions, toy data is needed. This is created by assuming some parameters for the basemodel, calculating the response, and then adding random noise.<jupyter_code>%matplotlib inline

import satlas as s
import numpy as np
np.random.seed(0)  #Ensure the same random numbers each time

I = 1.0
J = [1.0, 2.0]

ABC = [100, 200, 100, 200, 0, 0]
fwhm = [10, 10]
centroid = 500
scale = 100

basemodel = s.HFSModel(I, J, ABC, centroid, fwhm=fwhm, scale=scale, background_params=[10], use_saturation=False, use_racah=True)

frequency_range = (min(basemodel.locations) - 100, max(basemodel.locations) + 100)
frequency_range = np.linspace(frequency_range[0], frequency_range[1], 200)

data = np.floor(basemodel(frequency_range) + basemodel(frequency_range)**0.5 * np.random.randn(len(frequency_range)))

success, message = s.chisquare_spectroscopic_fit(basemodel, frequency_range, data)
print(success)
print(message)
print('')
basemodel.display_chisquare_fit(show_correl=False)
basemodel.plot(x=frequency_range, y=data)<jupyter_output>C:\Anaconda3\lib\site-packages\IPython\html.py:14: ShimWarning: The `IPython.html` package has been deprecated. You should import from `notebook` instead. `IPython.html.widgets` has moved to `ipywidgets`.
  "`IPython.html.widgets` has moved to `ipywidgets`.", ShimWarning)
C:\Anaconda3\lib\site-packages\matplotlib\__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.
  warnings.warn(self.msg_depr % (key, alt_key))
Chisquare fitting done: 96it [00:00, 711.07it/s]<jupyter_text>If *func=np.sqrt* is supplied, the result is<jupyter_code>basemodel = s.HFSModel(I, J, ABC, centroid, fwhm=fwhm, scale=scale, background_params=[10], use_saturation=False, use_racah=True)
success, message = s.chisquare_spectroscopic_fit(basemodel, frequency_range, data, func=np.sqrt)
print(success)
print(message)
print('')
basemodel.display_chisquare_fit(show_correl=False)
basemodel.plot(x=frequency_range, y=data)<jupyter_output>Chisquare fitting done: 126it [00:00, 818.13it/s]<jupyter_text>which is slightly different. Maximum Likelihood EstimationThe Maximum Likelihood Estimation (MLE) can be used to derive the chisquare method in the case of Gaussian uncertainties. When this is not the case (as it is for counting data, which has a Poisson distribution), a less simplified method has to be used.The MLE method works by minimizing the negative loglikelihood. This is calculated as$$-\mathcal{L}\left(\vec{\theta}\middle|x, y, \sigma\right) = \sum_i logp\left(\vec{\theta}\middle|x_i, y_i, \sigma_i\right)$$For the function *logp*, the standard choice is the loglikelihood derived from the Poisson distribution.To use this method, the fitting routine code has to be changed to<jupyter_code>basemodel = s.HFSModel(I, J, ABC, centroid, fwhm=fwhm, scale=scale, background_params=[10], use_saturation=False, use_racah=True)
success, message = s.likelihood_fit(basemodel, frequency_range, data)
print(success, message)
basemodel.display_mle_fit(show_correl=False)
basemodel.plot(x=frequency_range, y=data)<jupyter_output>Likelihood fitting done: 434it [00:00, 850.93it/s]
StarcoderdataJupyterStructuredCleanDedup