IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 2037


Ignore:
Timestamp:
Oct 8, 2004, 3:34:35 PM (22 years ago)
Author:
Paul Price
Message:

Added stuff on what to do with errors in statistics.
Fixed up LMM.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/doc/pslib/psLibADD.tex

    r1771 r2037  
    1 %%% $Id: psLibADD.tex,v 1.46 2004-09-10 00:40:42 price Exp $
     1%%% $Id: psLibADD.tex,v 1.47 2004-10-09 01:34:35 price Exp $
    22\documentclass[panstarrs]{panstarrs}
    33
     
    2727% version     Date         Description
    282800 & 2004 Mar 11 & Hacking \\ \hline
    29 01 & 2004 May 21 & Added section on 2D Chebyshev fitting. \\ \hline
     2901 & 2004 May 21 & Added section on 2D Chebyshev fitting, then removed. \\ \hline
    303002 & 2004 Jun 22 & modified stats specification \\ \hline
    313103--05 & ??? & ??? \\ \hline
     
    140140
    141141We define the following statistical terms, assuming there is a set of
    142 data elements $x_i$.
     142data elements $x_i$ with (standard) errors $\sigma_i$.
    143143
    144144\subparagraph{Mean}
    145145
    146 The mean is defined as:
    147 \begin{equation}
    148 \bar{x} = \frac{1}{N} \sum_{i = 1}^{N} x_i
    149 \end{equation}
     146The (weighted) mean is defined as:
     147\begin{equation}
     148\bar{x} = \sum_i \frac{x_i}{\sigma_i^2} \ / \ \sum_i \frac{1}{\sigma_i^2}
     149\end{equation}
     150
     151In the event that all the errors are identical, this reduces to the
     152standard definition of the mean.
    150153
    151154\subparagraph{Median}
     
    157160if the number of values is even.  This median should be avoided for
    158161samples which are large (e.g., $N > 10^4$ elements) as the basic
    159 robust median is quicker and more accurate. 
     162robust median is quicker and more accurate.  Errors are ignored when
     163calculating the sample median.
    160164
    161165\subparagraph{Upper and Lower Quartiles}
     
    171175to these values.  The sample quartiles should be avoided for samples
    172176which are large (e.g., $N > 10^4$ elements) as the robust quartiles
    173 are quicker and more accurate.
     177are quicker and more accurate.  Errors are ignored when calculating
     178the sample quartiles.
    174179
    175180\subparagraph{Standard Deviation}
     
    187192\sigma = \sqrt{\frac{1}{N - 1} \left[ \sum_{i = 1}^{N} (x_i - \bar{x})^2 - \frac{1}{N} \left(\sum_{i = 1}^{N} (x_i - \bar{x})\right)^2 \ \right]}
    188193\end{equation}
     194
     195If the errors are known, then the sample standard deviation is:
     196
     197\begin{equation}
     198\sigma = \left( \sum_i \frac{1}{\sigma_i^2} \right) ^{-1/2}
     199\end{equation}
     200
    189201
    190202\paragraph{Clipped Statistics}
     
    220232\end{enumerate}
    221233
     234If the errors in the input values are known, then the clips are made
     235on the basis of the errors in the input values instead of the standard
     236deviation of the sample: values are excluded for which $|x_i -
     237\bar{x}| > k \sigma_i$.
     238
    222239\paragraph{Robust Statistics}
    223240
     
    266283bins in the vicinity of the upper and lower quartile points.
    267284
     285If the errors in the input values are known, then the same approach is
     286used, except that the histograms become probability density functions
     287(PDFs).  In this case, the input values are spread out, so that they
     288do not simply contribute a single unit to the histogram, but rather
     289contribute a fraction of a value, equivalent to the weight.  In the
     290interests of speed, a boxcar PDF may be used to represent each input
     291value (as opposed to a Gaussian), where the width is equal to the
     292variance (the square of the error) and each input value contributes
     293constant area.  Then the mean, median, mode, standard deviation and
     294quartiles are estimated in the same manner as above.
     295
    268296
    269297\subsubsection{Matrix Operations}
     
    456484(section~\ref{LUdecomp}).
    457485
    458 \paragraph{Non-linear Fitting: Levenberg-Marquardt Method}
    459 
    460 For models in which the system of equations defined by the partial
    461 derivatives cannot be solved with the linear technique, other options
    462 are necessary.  The Levenberg-Marquardt Method (LMM; see NR \S 15.5)
    463 may be used for these situations.  In LMM, we make a guess at the
    464 input parameters, measure the $\chi^2$, vary the parameters by a
    465 particular choice based on the gradient, measure the $\chi^2$ again,
    466 and adjust the parameters and the parameter varient based on the
    467 results.
    468 
    469 Given a set of $N$ data values $y_i$ with errors $\sigma_i$, dependent
    470 on values $x_i$, we would like to find the parameters $a_k$ of the
    471 function $f(x_i; a_k)$ which minimize the $\chi^2$, defined in the
    472 usual manner (\ref{chisq}).  We start with a set of parameter guesses,
    473 $a_k$.  We calculate the gradient $\beta_k$ and the Hessian matrix
    474 $\alpha_{j,k}$ at this parameter selection as follows:
    475 
    476 \begin{eqnarray}
    477 \beta_k & = & \sum_{i=1}^{N} \frac{\partial f(x_i)}{\partial a_k} \frac{(y_i - f(x_i))}{\sigma_i^2} \\
    478 \alpha_{j,k} & = & \sum_{i=1}^{N} \frac{\partial f(x_i)}{\partial a_k} \frac{\partial f(x_i)}{\partial a_j} \frac{1}{\sigma_i^2}
    479 \end{eqnarray}
    480 %
     486\subsubsection{Non-linear Minimization}
     487
     488\paragraph{Levenberg-Marquardt Method}
     489
     490In the Levenberg-Marquardt Method (LMM; see NR \S 15.5), we make a
     491guess at the input parameters, evaluate the function of interest, vary
     492the parameters by a particular choice based on the gradient, evaluate
     493the function again, and adjust the parameters and the parameter
     494varient based on the results.
     495
     496The LMM only works if the second derivative of the function can be
     497considered negligible, as in the case of minimizing $\chi^2$.
     498
     499Given some ordinates, $x_i$, we would like to find the parameters,
     500$a_k$, of the function $f(x_i; a_k)$ which minimize $\chi^2$ for some
     501measurements, $y_i$ and associated errors, $\sigma_i$.  We start with
     502a set of parameter guesses, $a_k$.  We calculate the gradient
     503$\beta_k$ and the Hessian matrix $\alpha_{j,k}$ at this parameter
     504selection as follows:
     505\begin{eqnarray}
     506\beta_k & = & \frac{\partial \chi^2}{\partial a_k} \\
     507\alpha_{j,k} & = & \sum_i \frac{1}{\sigma_i^2} \frac{\partial f(c_i)}{\partial a_k} \frac{\partial f(c_i)}{\partial a_j}
     508\end{eqnarray}
     509
    481510We now define the new parameter guess for $a_k$ based on the gradient
    482511and Hessian by defining $A_{j,k}$ as a variant on $\alpha_{j,k}$ as
     
    493522\end{equation}
    494523%
    495 where $a^\prime_k$ represents our new attempt at a parameter
    496 guess. We use this parameter set to calculate $\chi^2$.  If the new
    497 value of $\chi^2$ is lower than the previous guess, we accept this new
    498 set of parameters and decrease $\lambda$ by a factor of 10, otherwise
    499 we keep the old set, and increase the value of $\lambda$ by a factor
    500 of 10.  We repeat this process until the value of the reduced $\chi^2$
    501 changes by much less than 1.0.  The resulting values of $a_k$ are the
    502 best-fit parameters for the system.  If the errors are normally
    503 distributed, the formal errors on the parameters are then calculated
    504 by setting $\lambda = 0$ and calculating the covarience matrix
    505 $C_{i,j}$, the inverse of the matrix $\alpha_{j,k}$.
    506 %
    507 The covariance matrix allows simple calculation of the confidence
    508 limits of the parameters.
    509 
     524where $a^\prime_k$ represents our new attempt at a parameter guess. We
     525use this parameter set to evaluate the function.  If the new value of
     526the function is lower than the previous guess, we accept this new set
     527of parameters and decrease $\lambda$ by a factor of 10, otherwise we
     528keep the old set, and increase the value of $\lambda$ by a factor of
     52910.  We repeat this process until the value of the function changes by
     530much less than the tolerance.  The resulting values of $a_k$ are the
     531best-fit parameters for the system.
     532
     533The covariance matrix, $C_{i,j}$, which is the inverse of the matrix
     534$\alpha_{j,k}$ allows simple calculation of the confidence limits of
     535the parameters.
     536
     537
     538%If the errors are normally distributed, the formal errors on the
     539%parameters are then calculated by setting $\lambda = 0$ and
     540%calculating the covarience matrix $C_{i,j}$, the inverse of the matrix
     541%$\alpha_{j,k}$.
    510542%The independent 68.3\% confidence limit on parameter $a_k$ is then
    511543%$\sqrt{C_{k,k}}$.  Confidence contours for sets of parameters may be
     
    533565
    534566
    535 \paragraph{Non-linear fitting: Powell's method}
     567\paragraph{Powell's method}
    536568
    537569Powell's method is a type of ``Direction Set'' methods in
Note: See TracChangeset for help on using the changeset viewer.