Line data Source code
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
11 : MODULE 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 :
94 : CONTAINS
95 :
96 : ! ------------------------------------------------------------------
97 :
98 4 : 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 2 : REAL(dp) :: sigdat, nx, sx, sxoss, sy, st2
108 202 : REAL(dp), DIMENSION(size(x)), TARGET :: t
109 2 : REAL(dp) :: aa, bb, cchi2
110 : LOGICAL :: mod2
111 2 : REAL(dp) :: mx, my, sx2, sy2, sxy, sxy2, syx2, ssigb
112 : !REAL(dp) :: r
113 :
114 2 : if (size(x) /= size(y)) stop 'linfit_dp: size(x) /= size(y)'
115 8 : if (.not. allocated(linfit_dp)) allocate(linfit_dp(size(x)))
116 2 : if (present(model2)) then
117 2 : mod2 = model2
118 : else
119 : mod2 = .false.
120 : end if
121 :
122 2 : if (mod2) then
123 1 : nx = real(size(x), dp)
124 101 : mx = sum(x) / nx
125 101 : my = sum(y) / nx
126 101 : sx2 = sum((x - mx) * (x - mx))
127 101 : sy2 = sum((y - my) * (y - my))
128 101 : sxy = sum((x - mx) * (y - my))
129 : !r = sxy / sqrt(sx2) / sqrt(sy2)
130 1 : bb = sign(sqrt(sy2 / sx2), sxy)
131 1 : aa = my - bb * mx
132 101 : linfit_dp(:) = aa + bb * x(:)
133 1 : if (present(a)) a = aa
134 1 : if (present(b)) b = bb
135 1 : if (present(chi2)) then
136 101 : t(:) = y(:) - linfit_dp(:)
137 101 : chi2 = dot_product(t, t)
138 : end if
139 1 : if (present(siga) .or. present(sigb)) then
140 1 : syx2 = (sy2 - sxy * sxy / sx2) / (nx - 2.0_dp)
141 1 : sxy2 = (sx2 - sxy * sxy / sy2) / (nx - 2.0_dp)
142 : ! syx2 should be >0
143 : ! ssigb = sqrt(syx2/sx2)
144 1 : ssigb = sqrt(abs(syx2) / sx2)
145 1 : if (present(sigb)) sigb = ssigb
146 1 : if (present(siga)) then
147 : ! siga = sqrt(syx2*(1.0_dp/nX+mx*mx/sx2))
148 1 : 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 1 : siga = sqrt(siga * siga + bb * bb * abs(sxy2) / nx)
153 : end if
154 : end if
155 : else
156 1 : nx = real(size(x), dp)
157 101 : sx = sum(x)
158 101 : sy = sum(y)
159 1 : sxoss = sx / nx
160 101 : t(:) = x(:) - sxoss
161 101 : bb = dot_product(t, y)
162 101 : st2 = dot_product(t, t)
163 1 : bb = bb / st2
164 1 : aa = (sy - sx * bb) / nx
165 101 : linfit_dp(:) = aa + bb * x(:)
166 1 : if (present(a)) a = aa
167 1 : if (present(b)) b = bb
168 1 : if (present(chi2) .or. present(siga) .or. present(sigb)) then
169 101 : t(:) = y(:) - linfit_dp(:)
170 101 : cchi2 = dot_product(t, t)
171 1 : if (present(chi2)) chi2 = cchi2
172 1 : if (present(siga) .or. present(sigb)) then
173 1 : sigdat = sqrt(cchi2 / (size(x) - 2))
174 1 : if (present(siga)) then
175 1 : siga = sqrt((1.0_dp + sx * sx / (nx * st2)) / nx)
176 1 : siga = siga * sigdat
177 : end if
178 1 : if (present(sigb)) then
179 1 : sigb = sqrt(1.0_dp / st2)
180 1 : sigb = sigb * sigdat
181 : end if
182 : end if
183 : end if
184 : end if
185 :
186 2 : END FUNCTION linfit_dp
187 :
188 :
189 4 : 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 2 : REAL(sp) :: sigdat, nx, sx, sxoss, sy, st2
199 202 : REAL(sp), DIMENSION(size(x)), TARGET :: t
200 2 : REAL(sp) :: aa, bb, cchi2
201 : LOGICAL :: mod2
202 2 : REAL(sp) :: mx, my, sx2, sy2, sxy, sxy2, syx2, ssigb
203 : !REAL(sp) :: r
204 :
205 2 : if (size(x) /= size(y)) stop 'linfit_sp: size(x) /= size(y)'
206 8 : if (.not. allocated(linfit_sp)) allocate(linfit_sp(size(x)))
207 2 : if (present(model2)) then
208 2 : mod2 = model2
209 : else
210 : mod2 = .false.
211 : end if
212 :
213 2 : if (mod2) then
214 1 : nx = real(size(x), sp)
215 101 : mx = sum(x) / nx
216 101 : my = sum(y) / nx
217 101 : sx2 = sum((x - mx) * (x - mx))
218 101 : sy2 = sum((y - my) * (y - my))
219 101 : sxy = sum((x - mx) * (y - my))
220 : !r = sxy / sqrt(sx2) / sqrt(sy2)
221 1 : bb = sign(sqrt(sy2 / sx2), sxy)
222 1 : aa = my - bb * mx
223 101 : linfit_sp(:) = aa + bb * x(:)
224 1 : if (present(a)) a = aa
225 1 : if (present(b)) b = bb
226 1 : if (present(chi2)) then
227 101 : t(:) = y(:) - linfit_sp(:)
228 101 : chi2 = dot_product(t, t)
229 : end if
230 1 : if (present(siga) .or. present(sigb)) then
231 1 : syx2 = (sy2 - sxy * sxy / sx2) / (nx - 2.0_sp)
232 1 : sxy2 = (sx2 - sxy * sxy / sy2) / (nx - 2.0_sp)
233 : ! syx2 should be >0
234 : ! ssigb = sqrt(syx2/sx2)
235 1 : ssigb = sqrt(abs(syx2) / sx2)
236 1 : if (present(sigb)) sigb = ssigb
237 1 : if (present(siga)) then
238 : ! siga = sqrt(syx2*(1.0_sp/nX+mx*mx/sx2))
239 1 : 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 1 : siga = sqrt(siga * siga + bb * bb * abs(sxy2) / nx)
244 : end if
245 : end if
246 : else
247 1 : nx = real(size(x), sp)
248 101 : sx = sum(x)
249 101 : sy = sum(y)
250 1 : sxoss = sx / nx
251 101 : t(:) = x(:) - sxoss
252 101 : bb = dot_product(t, y)
253 101 : st2 = dot_product(t, t)
254 1 : bb = bb / st2
255 1 : aa = (sy - sx * bb) / nx
256 101 : linfit_sp(:) = aa + bb * x(:)
257 1 : if (present(a)) a = aa
258 1 : if (present(b)) b = bb
259 1 : if (present(chi2) .or. present(siga) .or. present(sigb)) then
260 101 : t(:) = y(:) - linfit_sp(:)
261 101 : cchi2 = dot_product(t, t)
262 1 : if (present(chi2)) chi2 = cchi2
263 1 : if (present(siga) .or. present(sigb)) then
264 1 : sigdat = sqrt(cchi2 / (size(x) - 2))
265 1 : if (present(siga)) then
266 1 : siga = sqrt((1.0_sp + sx * sx / (nx * st2)) / nx)
267 1 : siga = siga * sigdat
268 : end if
269 1 : if (present(sigb)) then
270 1 : sigb = sqrt(1.0_sp / st2)
271 1 : sigb = sigb * sigdat
272 : end if
273 : end if
274 : end if
275 : end if
276 :
277 2 : END FUNCTION linfit_sp
278 :
279 : ! ------------------------------------------------------------------
280 :
281 : END MODULE mo_linfit
|
