linfitex.pro
3.57 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
;+
; NAME:
; LINFITEX
;
; AUTHOR:
; Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
; craigm@lheamail.gsfc.nasa.gov
; UPDATED VERSIONs can be found on my WEB PAGE:
; http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
; Model function for fitting line with errors in X and Y
;
; MAJOR TOPICS:
; Curve and Surface Fitting
;
; CALLING SEQUENCE:
; parms = MPFIT('LINFITEX', start_parms, $
; FUNCTARGS={X: X, Y: Y, SIGMA_X: SIGMA_X, SIGMA_Y: SIGMA_Y}, $
; ...)
;
; DESCRIPTION:
;
; LINFITEX is a model function to be used with MPFIT in order to
; fit a line to data with errors in both "X" and "Y" directions.
; LINFITEX follows the methodology of Numerical Recipes, in the
; section entitled, "Straight-Line Data with Errors in Both
; Coordinates."
;
; The user is not meant to call LINFITEX directly. Rather, the
; should pass LINFITEX as a user function to MPFIT, and MPFIT will in
; turn call LINFITEX.
;
; Each data point will have an X and Y position, as well as an error
; in X and Y, denoted SIGMA_X and SIGMA_Y. The user should pass
; these values using the FUNCTARGS convention, as shown above. I.e.
; the FUNCTARGS keyword should be set to a single structure
; containing the fields "X", "Y", "SIGMA_X" and "SIGMA_Y". Each
; field should have a vector of the same length.
;
; Upon return from MPFIT, the best fit parameters will be,
; P[0] - Y-intercept of line on the X=0 axis.
; P[1] - slope of the line
;
; NOTE that LINFITEX requires that AUTODERIVATIVE=1, i.e. MPFIT
; should compute the derivatives associated with each parameter
; numerically.
;
; INPUTS:
; P - parameters of the linear model, as described above.
;
; KEYWORD INPUTS:
; (as described above, these quantities should be placed in
; a FUNCTARGS structure)
; X - vector, X position of each data point
; Y - vector, Y position of each data point
; SIGMA_X - vector, X uncertainty of each data point
; SIGMA_Y - vector, Y uncertainty of each data point
;
; RETURNS:
; Returns a vector of residuals, of the same size as X.
;
; EXAMPLE:
;
; ; X and Y values
; XS = [2.9359964E-01,1.0125043E+00,2.5900450E+00,2.6647639E+00,3.7756164E+00,4.0297413E+00,4.9227958E+00,6.4959011E+00]
; YS = [6.0932738E-01,1.3339731E+00,1.3525699E+00,1.4060204E+00,2.8321848E+00,2.7798350E+00,2.0494456E+00,3.3113062E+00]
;
; ; X and Y errors
; XE = [1.8218818E-01,3.3440986E-01,3.7536234E-01,4.5585755E-01,7.3387712E-01,8.0054945E-01,6.2370265E-01,6.7048335E-01]
; YE = [8.9751285E-01,6.4095122E-01,1.1858428E+00,1.4673588E+00,1.0045623E+00,7.8527629E-01,1.2574003E+00,1.0080348E+00]
;
; ; Best fit line
; p = mpfit('LINFITEX', [1d, 1d], $
; FUNCTARGS={X: XS, Y: YS, SIGMA_X: XE, SIGMA_Y: YE}, $
; perror=dp, bestnorm=chi2)
; yfit = p[0] + p[1]*XS
;
;
; REFERENCES:
;
; Press, W. H. 1992, *Numerical Recipes in C*, 2nd Ed., Cambridge
; University Press
;
; MODIFICATION HISTORY:
; Written, Feb 2009
; Documented, 14 Apr 2009, CM
; $Id: linfitex.pro,v 1.3 2009/04/15 04:17:52 craigm Exp $
;
;-
; Copyright (C) 2009, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-
;
function linfitex, p, $
x=x, y=y, sigma_x=sigma_x, sigma_y=sigma_y, $
_EXTRA=extra
a = p[0] ;; Intercept
b = p[1] ;; Slope
f = a + b*x
resid = (y - f)/sqrt(sigma_y^2 + (b*sigma_x)^2)
return, resid
end