0.6.2-dev0
FORCES
FORtran lib for Comp. Env. Sys.
Loading...
Searching...
No Matches
mo_linfit.f90
Go to the documentation of this file.
1!> \file mo_linfit.f90
2!> \brief \copybrief mo_linfit
3!> \details \copydetails mo_linfit
4
5!> \brief Fitting a straight line.
6!> \details This module provides a routine to fit a straight line with model I or model II regression.
7!> \authors Matthias Cuntz
8!> \date Mar 2011
9!> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
10!! FORCES is released under the LGPLv3+ license \license_note
11MODULE mo_linfit
12
13 USE mo_kind, ONLY : sp, dp
14
15 Implicit NONE
16
17 PUBLIC :: linfit ! Fitting straight line (without error bars on input), Model I or Model II
18
19 ! ------------------------------------------------------------------
20
21 ! NAME
22 ! linfit
23
24 ! PURPOSE
25 !> \brief Fits a straight line to input data by minimizing chi^2.
26
27 !> \details Given a set of data points x(1:ndata), y(1:ndata),
28 !> fit them to a straight line \f$ y = a+bx \f$ by minimizing chi2.\n
29 !> Model I minimizes y vs. x while Model II takes the geometric mean of y vs. x and x vs. y.
30 !> Returned is the fitted line at x.
31 !> Optional returns are a, b and their respective probable uncertainties siga and sigb,
32 !> and the chi-square chi2.
33
34 ! CALLING SEQUENCE
35 ! out = linfit(x, y, a=a, b=b, siga=siga, sigb=sigb, chi2=chi2, model2=model2)
36
37 ! INTENT(IN)
38 !> \param[in] "real(sp/dp) :: x(:)" 1D-array with input x
39 !> \param[in] "real(sp/dp) :: y(:)" 1D-array with input y
40
41 ! INTENT(INOUT)
42 ! None
43
44 ! INTENT(OUT)
45 ! None
46
47 ! INTENT(IN), OPTIONAL
48 !> \param[in] "logical, optional :: model2" If present, use geometric mean regression
49 !> instead of ordinary least square
50
51 ! INTENT(INOUT), OPTIONAL
52 ! None
53
54 ! INTENT(OUT), OPTIONAL
55 !> \param[out] "real(sp/dp), dimension(M) :: a" intercept
56 !> \param[out] "real(sp/dp), dimension(M) :: b" slope
57 !> \param[out] "real(sp/dp), dimension(M) :: siga" error on intercept
58 !> \param[out] "real(sp/dp), dimension(M) :: sigb" error on slope
59 !> \param[out] "real(sp/dp) :: chisq" Minimum chi^2
60
61 ! RETURN
62 !> \return real(sp/dp), dimension(:), allocatable :: out — fitted values at x(:).
63
64 ! RESTRICTIONS
65 ! None
66
67 ! EXAMPLE
68 ! ytmp = linfit(x, y, a=inter, b=slope, model2=.true.)
69
70 ! LITERATURE
71 ! Model I follows closely
72 ! Press WH, Teukolsky SA, Vetterling WT, & Flannery BP - Numerical Recipes in Fortran 90 -
73 ! The Art of Parallel Scientific Computing, 2nd Edition, Volume 2 of Fortran Numerical Recipes,
74 ! Cambridge University Press, UK, 1996
75 ! For Model II (geometric mean regression)
76 ! Robert R. Sokal and F. James Rohlf, Biometry: the principle and practice of statistics
77 ! in biological research, Freeman & Co., ISBN 0-7167-2411-1
78
79 ! HISTORY
80 ! Written Matthias Cuntz, March 2011 - following closely the routine fit of Numerical Recipes
81 ! - linfit Model II: geometric mean regression
82 ! Modified Matthias Cuntz, Nov 2011 - sp, dp
83 ! - documentation
84 INTERFACE linfit
85 MODULE PROCEDURE linfit_sp, linfit_dp
86 END INTERFACE linfit
87
88 ! ------------------------------------------------------------------
89
90 PRIVATE
91
92 ! ------------------------------------------------------------------
93
94CONTAINS
95
96 ! ------------------------------------------------------------------
97
98 FUNCTION linfit_dp(x, y, a, b, siga, sigb, chi2, model2)
99
100 IMPLICIT NONE
101
102 REAL(dp), DIMENSION(:), INTENT(IN) :: x, y
103 REAL(dp), OPTIONAL, INTENT(OUT) :: a, b, siga, sigb, chi2
104 LOGICAL, OPTIONAL, INTENT(IN) :: model2
105 REAL(dp), DIMENSION(:), allocatable :: linfit_dp
106
107 REAL(dp) :: sigdat, nx, sx, sxoss, sy, st2
108 REAL(dp), DIMENSION(size(x)), TARGET :: t
109 REAL(dp) :: aa, bb, cchi2
110 LOGICAL :: mod2
111 REAL(dp) :: mx, my, sx2, sy2, sxy, sxy2, syx2, ssigb
112 !REAL(dp) :: r
113
114 if (size(x) /= size(y)) stop 'linfit_dp: size(x) /= size(y)'
115 if (.not. allocated(linfit_dp)) allocate(linfit_dp(size(x)))
116 if (present(model2)) then
117 mod2 = model2
118 else
119 mod2 = .false.
120 end if
121
122 if (mod2) then
123 nx = real(size(x), dp)
124 mx = sum(x) / nx
125 my = sum(y) / nx
126 sx2 = sum((x - mx) * (x - mx))
127 sy2 = sum((y - my) * (y - my))
128 sxy = sum((x - mx) * (y - my))
129 !r = sxy / sqrt(sx2) / sqrt(sy2)
130 bb = sign(sqrt(sy2 / sx2), sxy)
131 aa = my - bb * mx
132 linfit_dp(:) = aa + bb * x(:)
133 if (present(a)) a = aa
134 if (present(b)) b = bb
135 if (present(chi2)) then
136 t(:) = y(:) - linfit_dp(:)
137 chi2 = dot_product(t, t)
138 end if
139 if (present(siga) .or. present(sigb)) then
140 syx2 = (sy2 - sxy * sxy / sx2) / (nx - 2.0_dp)
141 sxy2 = (sx2 - sxy * sxy / sy2) / (nx - 2.0_dp)
142 ! syx2 should be >0
143 ! ssigb = sqrt(syx2/sx2)
144 ssigb = sqrt(abs(syx2) / sx2)
145 if (present(sigb)) sigb = ssigb
146 if (present(siga)) then
147 ! siga = sqrt(syx2*(1.0_dp/nX+mx*mx/sx2))
148 siga = sqrt(abs(syx2) * (1.0_dp / nx + mx * mx / sx2))
149 ! Add Extra Term for Error in xmean which is not in Sokal & Rohlf.
150 ! They take the error estimate of the chi-squared error for a.
151 ! siga = sqrt(siga*siga + bb*bb*sxy2/nx)
152 siga = sqrt(siga * siga + bb * bb * abs(sxy2) / nx)
153 end if
154 end if
155 else
156 nx = real(size(x), dp)
157 sx = sum(x)
158 sy = sum(y)
159 sxoss = sx / nx
160 t(:) = x(:) - sxoss
161 bb = dot_product(t, y)
162 st2 = dot_product(t, t)
163 bb = bb / st2
164 aa = (sy - sx * bb) / nx
165 linfit_dp(:) = aa + bb * x(:)
166 if (present(a)) a = aa
167 if (present(b)) b = bb
168 if (present(chi2) .or. present(siga) .or. present(sigb)) then
169 t(:) = y(:) - linfit_dp(:)
170 cchi2 = dot_product(t, t)
171 if (present(chi2)) chi2 = cchi2
172 if (present(siga) .or. present(sigb)) then
173 sigdat = sqrt(cchi2 / (size(x) - 2))
174 if (present(siga)) then
175 siga = sqrt((1.0_dp + sx * sx / (nx * st2)) / nx)
176 siga = siga * sigdat
177 end if
178 if (present(sigb)) then
179 sigb = sqrt(1.0_dp / st2)
180 sigb = sigb * sigdat
181 end if
182 end if
183 end if
184 end if
185
186 END FUNCTION linfit_dp
187
188
189 FUNCTION linfit_sp(x, y, a, b, siga, sigb, chi2, model2)
190
191 IMPLICIT NONE
192
193 REAL(sp), DIMENSION(:), INTENT(IN) :: x, y
194 REAL(sp), OPTIONAL, INTENT(OUT) :: a, b, siga, sigb, chi2
195 LOGICAL, OPTIONAL, INTENT(IN) :: model2
196 REAL(sp), DIMENSION(:), allocatable :: linfit_sp
197
198 REAL(sp) :: sigdat, nx, sx, sxoss, sy, st2
199 REAL(sp), DIMENSION(size(x)), TARGET :: t
200 REAL(sp) :: aa, bb, cchi2
201 LOGICAL :: mod2
202 REAL(sp) :: mx, my, sx2, sy2, sxy, sxy2, syx2, ssigb
203 !REAL(sp) :: r
204
205 if (size(x) /= size(y)) stop 'linfit_sp: size(x) /= size(y)'
206 if (.not. allocated(linfit_sp)) allocate(linfit_sp(size(x)))
207 if (present(model2)) then
208 mod2 = model2
209 else
210 mod2 = .false.
211 end if
212
213 if (mod2) then
214 nx = real(size(x), sp)
215 mx = sum(x) / nx
216 my = sum(y) / nx
217 sx2 = sum((x - mx) * (x - mx))
218 sy2 = sum((y - my) * (y - my))
219 sxy = sum((x - mx) * (y - my))
220 !r = sxy / sqrt(sx2) / sqrt(sy2)
221 bb = sign(sqrt(sy2 / sx2), sxy)
222 aa = my - bb * mx
223 linfit_sp(:) = aa + bb * x(:)
224 if (present(a)) a = aa
225 if (present(b)) b = bb
226 if (present(chi2)) then
227 t(:) = y(:) - linfit_sp(:)
228 chi2 = dot_product(t, t)
229 end if
230 if (present(siga) .or. present(sigb)) then
231 syx2 = (sy2 - sxy * sxy / sx2) / (nx - 2.0_sp)
232 sxy2 = (sx2 - sxy * sxy / sy2) / (nx - 2.0_sp)
233 ! syx2 should be >0
234 ! ssigb = sqrt(syx2/sx2)
235 ssigb = sqrt(abs(syx2) / sx2)
236 if (present(sigb)) sigb = ssigb
237 if (present(siga)) then
238 ! siga = sqrt(syx2*(1.0_sp/nX+mx*mx/sx2))
239 siga = sqrt(abs(syx2) * (1.0_sp / nx + mx * mx / sx2))
240 ! Add Extra Term for Error in xmean which is not in Sokal & Rohlf.
241 ! They take the error estimate of the chi-squared error for a.
242 ! siga = sqrt(siga*siga + bb*bb*sxy2/nx)
243 siga = sqrt(siga * siga + bb * bb * abs(sxy2) / nx)
244 end if
245 end if
246 else
247 nx = real(size(x), sp)
248 sx = sum(x)
249 sy = sum(y)
250 sxoss = sx / nx
251 t(:) = x(:) - sxoss
252 bb = dot_product(t, y)
253 st2 = dot_product(t, t)
254 bb = bb / st2
255 aa = (sy - sx * bb) / nx
256 linfit_sp(:) = aa + bb * x(:)
257 if (present(a)) a = aa
258 if (present(b)) b = bb
259 if (present(chi2) .or. present(siga) .or. present(sigb)) then
260 t(:) = y(:) - linfit_sp(:)
261 cchi2 = dot_product(t, t)
262 if (present(chi2)) chi2 = cchi2
263 if (present(siga) .or. present(sigb)) then
264 sigdat = sqrt(cchi2 / (size(x) - 2))
265 if (present(siga)) then
266 siga = sqrt((1.0_sp + sx * sx / (nx * st2)) / nx)
267 siga = siga * sigdat
268 end if
269 if (present(sigb)) then
270 sigb = sqrt(1.0_sp / st2)
271 sigb = sigb * sigdat
272 end if
273 end if
274 end if
275 end if
276
277 END FUNCTION linfit_sp
278
279 ! ------------------------------------------------------------------
280
281END MODULE mo_linfit
mo_linfit::linfit
Fits a straight line to input data by minimizing chi^2.
Definition mo_linfit.f90:84
mo_kind
Define number representations.
Definition mo_kind.F90:17
mo_kind::sp
integer, parameter sp
Single Precision Real Kind.
Definition mo_kind.F90:44
mo_kind::dp
integer, parameter dp
Double Precision Real Kind.
Definition mo_kind.F90:46
mo_linfit
Fitting a straight line.
Definition mo_linfit.f90:11