Editing Gaussian function (section)

== Estimation of parameters ==
{{see also|Normal distribution#Estimation of parameters}}

A number of fields such as [[Photometry (astronomy)|stellar photometry]], [[Gaussian beam]] characterization, and [[emission spectrum#Emission spectroscopy|emission/absorption line spectroscopy]] work with sampled Gaussian functions and need to accurately estimate the height, position, and width parameters of the function. There are three unknown parameters for a 1D Gaussian function (''a'', ''b'', ''c'') and five for a 2D Gaussian function <math>(A; x_0,y_0; \sigma_X,\sigma_Y)</math>. 

The most common method for estimating the Gaussian parameters is to take the logarithm of the data and [[polynomial fitting|fit a parabola]] to the resulting data set.<ref name="Caruana Searle Heller Shupack 1986 pp. 1162–1167">{{cite journal | last1=Caruana | first1=Richard A. | last2=Searle | first2=Roger B. | last3=Heller | first3=Thomas. | last4=Shupack | first4=Saul I. | title=Fast algorithm for the resolution of spectra | journal=Analytical Chemistry | publisher=American Chemical Society (ACS) | volume=58 | issue=6 | year=1986 | issn=0003-2700 | doi=10.1021/ac00297a041 | pages=1162–1167}}</ref><ref name="Guo">[https://dx.doi.org/10.1109/MSP.2011.941846 Hongwei Guo, "A simple algorithm for fitting a Gaussian function," IEEE Sign. Proc. Mag. 28(9): 134-137 (2011).]</ref> While this provides a simple [[curve fitting]] procedure, the resulting algorithm may be biased by excessively weighting small data values, which can produce large errors in the profile estimate. One can partially compensate for this problem through [[weighted least squares]] estimation, reducing the weight of small data values, but this too can be biased by allowing the tail of the Gaussian to dominate the fit. In order to remove the bias, one can instead use an [[iteratively reweighted least squares]] procedure, in which the weights are updated at each iteration.<ref name="Guo" />
It is also possible to perform [[non-linear regression]] directly on the data, without involving the [[logarithmic data transformation]]; for more options, see [[probability distribution fitting]].

=== Parameter precision ===

Once one has an algorithm for estimating the Gaussian function parameters, it is also important to know how [[Accuracy and precision|precise]] those estimates are. Any [[least squares]] estimation algorithm can provide numerical estimates for the variance of each parameter (i.e., the variance of the estimated height, position, and width of the function). One can also use [[Cramér–Rao bound]] theory to obtain an analytical expression for the lower bound on the parameter variances, given certain assumptions about the data.<ref name="Hagen1">[https://dx.doi.org/10.1364/AO.46.005374 N. Hagen, M. Kupinski, and E. L. Dereniak, "Gaussian profile estimation in one dimension," Appl. Opt. 46:5374–5383 (2007)]</ref><ref name="Hagen2">[https://dx.doi.org/10.1364/AO.47.006842 N. Hagen and E. L. Dereniak, "Gaussian profile estimation in two dimensions," Appl. Opt. 47:6842–6851 (2008)]</ref>
# The noise in the measured profile is either [[Independent and identically-distributed random variables|i.i.d.]] Gaussian, or the noise is [[Poisson distribution|Poisson-distributed]].
# The spacing between each sampling (i.e. the distance between pixels measuring the data) is uniform.
# The peak is "well-sampled", so that less than 10% of the area or volume under the peak (area if a 1D Gaussian, volume if a 2D Gaussian) lies outside the measurement region.
# The width of the peak is much larger than the distance between sample locations (i.e. the detector pixels must be at least 5 times smaller than the Gaussian FWHM).
When these assumptions are satisfied, the following [[covariance matrix]] '''K''' applies for the 1D profile parameters <math>a</math>, <math>b</math>, and <math>c</math> under i.i.d. Gaussian noise and under Poisson noise:<ref name="Hagen1" />
<math display="block"> \mathbf{K}_{\text{Gauss}} = \frac{\sigma^2}{\sqrt{\pi} \delta_X Q^2} \begin{pmatrix} \frac{3}{2c} &0 &\frac{-1}{a} \\ 0 &\frac{2c}{a^2} &0 \\ \frac{-1}{a} &0 &\frac{2c}{a^2} \end{pmatrix} \ ,
\qquad
\mathbf{K}_\text{Poiss} = \frac{1}{\sqrt{2 \pi}} \begin{pmatrix} \frac{3a}{2c} &0 &-\frac{1}{2} \\ 0 &\frac{c}{a} &0 \\ -\frac{1}{2} &0 &\frac{c}{2a} \end{pmatrix} \ ,</math>
where <math>\delta_X</math> is the width of the pixels used to sample the function, <math>Q</math> is the quantum efficiency of the detector, and <math>\sigma</math> indicates the standard deviation of the measurement noise. Thus, the individual variances for the parameters are, in the Gaussian noise case,
<math display="block">\begin{align}
\operatorname{var} (a) &= \frac{3 \sigma^2}{2 \sqrt{\pi} \, \delta_X Q^2 c} \\
\operatorname{var} (b) &= \frac{2 \sigma^2 c}{\delta_X \sqrt{\pi} \, Q^2 a^2} \\
\operatorname{var} (c) &= \frac{2 \sigma^2 c}{\delta_X \sqrt{\pi} \, Q^2 a^2}
\end{align}</math>

and in the Poisson noise case,
<math display="block">\begin{align}
\operatorname{var} (a) &= \frac{3a}{2 \sqrt{2 \pi} \, c} \\
\operatorname{var} (b) &= \frac{c}{\sqrt{2 \pi} \, a} \\
\operatorname{var} (c) &= \frac{c}{2 \sqrt{2 \pi} \, a}.
\end{align} </math>

For the 2D profile parameters giving the amplitude <math>A</math>, position <math>(x_0,y_0)</math>, and width <math>(\sigma_X,\sigma_Y)</math> of the profile, the following covariance matrices apply:<ref name="Hagen2" />

<math display="block">\begin{align}
\mathbf{K}_\text{Gauss} = \frac{\sigma^2}{\pi \delta_X \delta_Y Q^2} & \begin{pmatrix}
\frac{2}{\sigma_X \sigma_Y} &0 &0 &\frac{-1}{A \sigma_Y} &\frac{-1}{A \sigma_X} \\
 0 &\frac{2 \sigma_X}{A^2 \sigma_Y} &0 &0 &0 \\
 0 &0 &\frac{2 \sigma_Y}{A^2 \sigma_X} &0 &0 \\
 \frac{-1}{A \sigma_y} &0 &0 &\frac{2 \sigma_X}{A^2 \sigma_y} &0 \\
 \frac{-1}{A \sigma_X} &0 &0 &0 &\frac{2 \sigma_Y}{A^2 \sigma_X}
\end{pmatrix} \\[6pt]
\mathbf{K}_{\operatorname{Poisson}} = \frac{1}{2 \pi} & \begin{pmatrix}
 \frac{3A}{\sigma_X \sigma_Y} &0 &0 &\frac{-1}{\sigma_Y} &\frac{-1}{\sigma_X} \\
 0 & \frac{\sigma_X}{A \sigma_Y} &0 &0 &0 \\ 0 &0 &\frac{\sigma_Y}{A \sigma_X} &0 &0 \\
 \frac{-1}{\sigma_Y} &0 &0 &\frac{2 \sigma_X}{3A \sigma_Y} &\frac{1}{3A} \\
 \frac{-1}{\sigma_X} &0 &0 &\frac{1}{3A} &\frac{2 \sigma_Y}{3A \sigma_X}
 \end{pmatrix}.
\end{align}</math>
where the individual parameter variances are given by the diagonal elements of the covariance matrix.