Index: /trunk/doc/pslib/psLibADD.tex =================================================================== --- /trunk/doc/pslib/psLibADD.tex (revision 2036) +++ /trunk/doc/pslib/psLibADD.tex (revision 2037) @@ -1,3 +1,3 @@ -%%% $Id: psLibADD.tex,v 1.46 2004-09-10 00:40:42 price Exp $ +%%% $Id: psLibADD.tex,v 1.47 2004-10-09 01:34:35 price Exp $ \documentclass[panstarrs]{panstarrs} @@ -27,5 +27,5 @@ % version Date Description 00 & 2004 Mar 11 & Hacking \\ \hline -01 & 2004 May 21 & Added section on 2D Chebyshev fitting. \\ \hline +01 & 2004 May 21 & Added section on 2D Chebyshev fitting, then removed. \\ \hline 02 & 2004 Jun 22 & modified stats specification \\ \hline 03--05 & ??? & ??? \\ \hline @@ -140,12 +140,15 @@ We define the following statistical terms, assuming there is a set of -data elements $x_i$. +data elements $x_i$ with (standard) errors $\sigma_i$. \subparagraph{Mean} -The mean is defined as: -\begin{equation} -\bar{x} = \frac{1}{N} \sum_{i = 1}^{N} x_i -\end{equation} +The (weighted) mean is defined as: +\begin{equation} +\bar{x} = \sum_i \frac{x_i}{\sigma_i^2} \ / \ \sum_i \frac{1}{\sigma_i^2} +\end{equation} + +In the event that all the errors are identical, this reduces to the +standard definition of the mean. \subparagraph{Median} @@ -157,5 +160,6 @@ if the number of values is even. This median should be avoided for samples which are large (e.g., $N > 10^4$ elements) as the basic -robust median is quicker and more accurate. +robust median is quicker and more accurate. Errors are ignored when +calculating the sample median. \subparagraph{Upper and Lower Quartiles} @@ -171,5 +175,6 @@ to these values. The sample quartiles should be avoided for samples which are large (e.g., $N > 10^4$ elements) as the robust quartiles -are quicker and more accurate. +are quicker and more accurate. Errors are ignored when calculating +the sample quartiles. \subparagraph{Standard Deviation} @@ -187,4 +192,11 @@ \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]} \end{equation} + +If the errors are known, then the sample standard deviation is: + +\begin{equation} +\sigma = \left( \sum_i \frac{1}{\sigma_i^2} \right) ^{-1/2} +\end{equation} + \paragraph{Clipped Statistics} @@ -220,4 +232,9 @@ \end{enumerate} +If the errors in the input values are known, then the clips are made +on the basis of the errors in the input values instead of the standard +deviation of the sample: values are excluded for which $|x_i - +\bar{x}| > k \sigma_i$. + \paragraph{Robust Statistics} @@ -266,4 +283,15 @@ bins in the vicinity of the upper and lower quartile points. +If the errors in the input values are known, then the same approach is +used, except that the histograms become probability density functions +(PDFs). In this case, the input values are spread out, so that they +do not simply contribute a single unit to the histogram, but rather +contribute a fraction of a value, equivalent to the weight. In the +interests of speed, a boxcar PDF may be used to represent each input +value (as opposed to a Gaussian), where the width is equal to the +variance (the square of the error) and each input value contributes +constant area. Then the mean, median, mode, standard deviation and +quartiles are estimated in the same manner as above. + \subsubsection{Matrix Operations} @@ -456,27 +484,28 @@ (section~\ref{LUdecomp}). -\paragraph{Non-linear Fitting: Levenberg-Marquardt Method} - -For models in which the system of equations defined by the partial -derivatives cannot be solved with the linear technique, other options -are necessary. The Levenberg-Marquardt Method (LMM; see NR \S 15.5) -may be used for these situations. In LMM, we make a guess at the -input parameters, measure the $\chi^2$, vary the parameters by a -particular choice based on the gradient, measure the $\chi^2$ again, -and adjust the parameters and the parameter varient based on the -results. - -Given a set of $N$ data values $y_i$ with errors $\sigma_i$, dependent -on values $x_i$, we would like to find the parameters $a_k$ of the -function $f(x_i; a_k)$ which minimize the $\chi^2$, defined in the -usual manner (\ref{chisq}). We start with a set of parameter guesses, -$a_k$. We calculate the gradient $\beta_k$ and the Hessian matrix -$\alpha_{j,k}$ at this parameter selection as follows: - -\begin{eqnarray} -\beta_k & = & \sum_{i=1}^{N} \frac{\partial f(x_i)}{\partial a_k} \frac{(y_i - f(x_i))}{\sigma_i^2} \\ -\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} -\end{eqnarray} -% +\subsubsection{Non-linear Minimization} + +\paragraph{Levenberg-Marquardt Method} + +In the Levenberg-Marquardt Method (LMM; see NR \S 15.5), we make a +guess at the input parameters, evaluate the function of interest, vary +the parameters by a particular choice based on the gradient, evaluate +the function again, and adjust the parameters and the parameter +varient based on the results. + +The LMM only works if the second derivative of the function can be +considered negligible, as in the case of minimizing $\chi^2$. + +Given some ordinates, $x_i$, we would like to find the parameters, +$a_k$, of the function $f(x_i; a_k)$ which minimize $\chi^2$ for some +measurements, $y_i$ and associated errors, $\sigma_i$. We start with +a set of parameter guesses, $a_k$. We calculate the gradient +$\beta_k$ and the Hessian matrix $\alpha_{j,k}$ at this parameter +selection as follows: +\begin{eqnarray} +\beta_k & = & \frac{\partial \chi^2}{\partial a_k} \\ +\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} +\end{eqnarray} + We now define the new parameter guess for $a_k$ based on the gradient and Hessian by defining $A_{j,k}$ as a variant on $\alpha_{j,k}$ as @@ -493,19 +522,22 @@ \end{equation} % -where $a^\prime_k$ represents our new attempt at a parameter -guess. We use this parameter set to calculate $\chi^2$. If the new -value of $\chi^2$ is lower than the previous guess, we accept this new -set of parameters and decrease $\lambda$ by a factor of 10, otherwise -we keep the old set, and increase the value of $\lambda$ by a factor -of 10. We repeat this process until the value of the reduced $\chi^2$ -changes by much less than 1.0. The resulting values of $a_k$ are the -best-fit parameters for the system. If the errors are normally -distributed, the formal errors on the parameters are then calculated -by setting $\lambda = 0$ and calculating the covarience matrix -$C_{i,j}$, the inverse of the matrix $\alpha_{j,k}$. -% -The covariance matrix allows simple calculation of the confidence -limits of the parameters. - +where $a^\prime_k$ represents our new attempt at a parameter guess. We +use this parameter set to evaluate the function. If the new value of +the function is lower than the previous guess, we accept this new set +of parameters and decrease $\lambda$ by a factor of 10, otherwise we +keep the old set, and increase the value of $\lambda$ by a factor of +10. We repeat this process until the value of the function changes by +much less than the tolerance. The resulting values of $a_k$ are the +best-fit parameters for the system. + +The covariance matrix, $C_{i,j}$, which is the inverse of the matrix +$\alpha_{j,k}$ allows simple calculation of the confidence limits of +the parameters. + + +%If the errors are normally distributed, the formal errors on the +%parameters are then calculated by setting $\lambda = 0$ and +%calculating the covarience matrix $C_{i,j}$, the inverse of the matrix +%$\alpha_{j,k}$. %The independent 68.3\% confidence limit on parameter $a_k$ is then %$\sqrt{C_{k,k}}$. Confidence contours for sets of parameters may be @@ -533,5 +565,5 @@ -\paragraph{Non-linear fitting: Powell's method} +\paragraph{Powell's method} Powell's method is a type of ``Direction Set'' methods in