reflectometry.reduction.wsolve

1 """ 2 Solve a potentially over-determined system with uncertainty in 3 the values. 4 5 Given: A x = y +/- dy 6 Use: s = wsolve(A,y,dy) 7 8 wsolve uses the singular value decomposition for increased accuracy. 9 Estimates the uncertainty for the solution from the scatter in the data. 10 11 The returned model object s provides: 12 13 s.x solution 14 s.std uncertainty estimate assuming no correlation 15 s.rnorm residual norm 16 s.DoF degrees of freedom 17 s.cov covariance matrix 18 s.ci(p) confidence intervals at point p 19 s.pi(p) prediction intervals at point p 20 s(p) predicted value at point p 21 22 Example 23 ======= 24 25 Weighted system:: 26 27 import numpy,wsolve 28 A = numpy.matrix("1,2,3;2,1,3;1,1,1",'d').A 29 xin = numpy.array([1,2,3],'d') 30 dy = numpy.array([0.2,0.01,0.1]) 31 y = numpy.random.normal(numpy.dot(A,xin),dy) 32 print A,y,dy 33 s = wsolve.wsolve(A,y,dy) 34 print "xin,x,dx", xin, s.x, s.std 35 36 Note there is a counter-intuitive result that scaling the estimated 37 uncertainty in the data does not affect the computed uncertainty in 38 the fit. This is the correct result --- if the data were indeed 39 selected from a process with ten times the uncertainty, you would 40 expect the scatter in the data to increase by a factor of ten as 41 well. When this new data set is fitted, it will show a computed 42 uncertainty increased by the same factor. Monte carlo simulations 43 bear this out. The conclusion is that the dataset carries its own 44 information about the variance in the data, and the weight vector 45 serves only to provide relative weighting between the points. 46 """ 47 48 # FIXME: test second example 49 # 50 # Example 2: weighted overdetermined system y = x1 + 2*x2 + 3*x3 + e 51 # 52 # A = fullfact([3,3,3]); xin=[1;2;3]; 53 # y = A*xin; dy = rand(size(y))/50; y+=dy.*randn(size(y)); 54 # [x,s] = wsolve(A,y,dy); 55 # dx = s.normr*sqrt(sumsq(inv(s.R'))'/s.df); 56 # res = [xin, x, dx] 57 58 59 import numpy as N 60 61 # Grab erfc from scipy if it is available; if not then we can only 62 # calculate confidence intervals for sigma = 1 63 try: 64 from scipy import stats 65 from scipy.special import erfc 66 except: 67 # If scipy is missing we will not be able to calculate confidence 68 # intervals or prediction intervals. 69 pass 70

71 -class LinearModel(object):

72 """ 73 Model evaluator for linear solution to Ax = y. 74 75 Computes a confidence interval (range of likely values for the 76 mean at x) or a prediction interval (range of likely values 77 seen when measuring at x). The prediction interval tells 78 you the width of the distribution at x. This should be the same 79 regardless of the number of measurements you have for the value 80 at x. The confidence interval tells you how well you know the 81 mean at x. It should get smaller as you increase the number of 82 measurements. Error bars in the physical sciences usually show 83 a 1-alpha confidence value of erfc(1/sqrt(2)), representing 84 a 1 sigma standandard deviation of uncertainty in the mean. 85 86 Confidence intervals for linear system are given by:: 87 88 x' p +/- sqrt( Finv(1-a,1,df) var(x' p) ) 89 90 where for confidence intervals:: 91 92 var(x' p) = sigma^2 (x' inv(A'A) x) 93 94 and for prediction intervals:: 95 96 var(x' p) = sigma^2 (1 + x' inv(A'A) x) 97 98 99 Stored properties:: 100 101 DoF = len(y)-len(x) = degrees of freedom 102 rnorm = 2-norm of the residuals y-Ax 103 x = solution to the equation Ax = y 104 105 Computed properties:: 106 107 cov = covariance matrix [ inv(A'A); O(n^3) ] 108 var = parameter variance [ diag(cov); O(n^2)] 109 std = standard deviation of parameters [ sqrt(var); O(n^2) ] 110 p = test statistic for chisquare goodness of fit [ chi2.sf; O(1) ] 111 112 Methods:: 113 114 ci(A,sigma=1): return confidence interval evaluated at A 115 pi(A,alpha=0.05): return prediction interval evaluated at A 116 117 """

118 - def __init__(self, x=None, DoF=None, SVinv=None, rnorm=None):

119 """ 120 121 """ 122 # V,S where USV' = A 123 self.x = x 124 self.DoF = DoF 125 self.rnorm = rnorm 126 self._SVinv = SVinv

127 128 # covariance matrix invC = A'A = (USV')'USV' = VSU'USV' = VSSV' 129 # C = inv(A'A) = inv(VSSV') = inv(V')inv(SS)inv(V) = Vinv(SS)V' 130 # diag(inv(A'A)) is sum of the squares of the columns inv(S) V' 131 # and is also the sum of the squares of the rows of V inv(S)

132 - def _cov(self):

133 # FIXME: don't know if we need to scale by C, but it will 134 # at least make things consistent 135 C = self.rnorm**2/self.DoF if self.DoF>0 else 1 136 return C * N.dot(self._SVinv,self._SVinv.T)

137 - def _var(self):

138 C = self.rnorm**2/self.DoF if self.DoF>0 else 1 139 return C * N.sum( self._SVinv**2, axis=1)

140 - def _std(self):

141 return N.sqrt(self._var())

142 - def _p(self):

143 from scipy import stats 144 return stats.chi2.sf(self.rnorm**2,self.DoF)

145 146 cov = property(_cov,doc="covariance matrix") 147 var = property(_var,doc="result variance") 148 std = property(_std,doc="result standard deviation") 149 p = property(_p,doc="probability of rejection") 150

151 - def _interval(self,X,alpha,pred):

152 """ 153 Helper for computing prediction/confidence intervals. 154 """ 155 156 # Comments from QR decomposition solution to Ax = y: 157 # 158 # Rather than A'A we have R from the QR decomposition of A, but 159 # R'R equals A'A. Note that R is not upper triangular since we 160 # have already multiplied it by the permutation matrix, but it 161 # is invertible. Rather than forming the product R'R which is 162 # ill-conditioned, we can rewrite x' inv(A'A) x as the equivalent 163 # x' inv(R) inv(R') x = t t', for t = x' inv(R) 164 # 165 # We have since switched to an SVD solver, which gives us 166 # 167 # invC = A'A = (USV')'USV' = VSU'USV' = VSSV' 168 # C = inv(A'A) = inv(VSSV') = inv(V')inv(SS)inv(V) 169 # = Vinv(SS)V' = Vinv(S) inv(S)V' 170 # 171 # Substituting, we get 172 # 173 # x' inv(A'A) x = t t', for t = x' Vinv(S) 174 # 175 # Since x is a vector, t t' is the inner product sum(t**2). 176 # Note that LAPACK allows us to do this simultaneously for many 177 # different x using sqrt(sum(T**2,axis=1)), with T = X' Vinv(S). 178 # 179 # Note: sqrt(F(1-a;1,df)) = T(1-a/2;df) 180 # 181 y = N.dot(X,self.x).ravel() 182 s = stats.t.ppf(1-alpha/2,self.DoF)*self.rnorm/N.sqrt(self.DoF) 183 t = N.dot(X,self._SVinv) 184 dy = s*N.sqrt(pred + N.sum( t**2, axis=1)) 185 return y,dy

186

187 - def __call__(self, A):

188 """ 189 Return the prediction for a linear system at points in the 190 rows of A. 191 """ 192 return N.dot(N.asarray(A),self.x)

193

194 - def ci(self, A, sigma=1):

195 """ 196 Compute the calculated values and the confidence intervals 197 for the linear model evaluated at A. 198 199 sigma=1 corresponds to a 1-sigma confidence interval 200 201 Confidence intervals are sometimes expressed as 1-alpha values, 202 where alpha = erfc(sigma/sqrt(2)). 203 """ 204 alpha = erfc(sigma/N.sqrt(2)) 205 return self._interval(N.asarray(A),alpha,0)

206

207 - def pi(self, A, p=0.05):

208 """ 209 Compute the calculated values and the prediction intervals 210 for the linear model evaluated at A. 211 212 p = 1-alpha = 0.05 corresponds to 95% prediction interval 213 """ 214 return self._interval(N.asarray(A),p,1)

215

216 -def wsolve(A,y,dy=1,rcond=1e-12):

217 """ 218 Given a linear system y = A*x + e(dy), estimates x,dx 219 220 A is an n x m array 221 y is an n x k array or vector of length n 222 dy is a scalar or an n x 1 array 223 x is a m x k array 224 """ 225 # The ugliness v[:,N.newaxis] transposes a vector 226 # The ugliness N.dot(a,b) is a*b for a,b matrices 227 # The ugliness vh.T.conj() is the hermitian transpose 228 229 # Make sure inputs are arrays 230 A,y,dy = N.asarray(A),N.asarray(y),N.asarray(dy) 231 result_dims = y.ndim 232 if dy.ndim == 1: dy = dy[:,N.newaxis] 233 if y.ndim == 1: y = y[:,N.newaxis] 234 235 # Apply weighting if dy is not a scalar 236 # If dy is a scalar, it cancels out of both sides of the equation 237 # Note: with A,dy arrays instead of matrices, A/dy operates element-wise 238 # Since dy is a row vector, this divides each row of A by the corresponding 239 # element of dy. 240 if dy.ndim == 2: A,y = A/dy,y/dy 241 242 # Singular value decomposition: A = U S V.H 243 # Since A is an array, U, S, VH are also arrays 244 # The zero indicates an economy decomposition, with u nxm rathern than nxn 245 u,s,vh = N.linalg.svd(A,0) 246 247 # FIXME what to do with ill-conditioned systems? 248 #if s[-1]<rcond*s[0]: raise ValueError, "matrix is singular" 249 #s[s<rcond*s[0]] = 0. # Can't do this because 1/s below will fail 250 251 # Solve: x = V inv(S) U.H y 252 # S diagonal elements => 1/S is inv(S) 253 # A*D, D diagonal multiplies each column of A by the corresponding diagonal 254 # D*A, D diagonal multiplies each row of A by the corresponding diagonal 255 # Computing V*inv(S) is slightly faster than inv(S)*U.H since V is smaller 256 # than U.H. Similarly, U.H*y is somewhat faster than V*U.H 257 SVinv = vh.T.conj()/s 258 Uy = N.dot(u.T.conj(), y) 259 x = N.dot(SVinv, Uy) 260 261 DoF = y.shape[0] - x.shape[0] 262 rnorm = N.linalg.norm(y - N.dot(A,x)) 263 264 return LinearModel(x=x, DoF=DoF, SVinv=SVinv, rnorm=rnorm)

265

266 -def _poly_matrix(x,degree,origin=False):

267 """ 268 Generate the matrix A used to fit a polynomial using a linear solver. 269 """ 270 if origin: 271 n = N.array(range(degree,0,-1)) 272 else: 273 n = N.array(range(degree,-1,-1)) 274 return N.asarray(x)[:,None]**n[None,:]

275

276 -class PolynomialModel(object):

277 """ 278 Model evaluator for best fit polynomial p(x) = y. 279 280 Stored properties:: 281 282 DoF = len(y)-len(x) = degrees of freedom 283 rnorm = 2-norm of the residuals y-Ax 284 coeff = coefficients 285 degree = polynomial degree 286 287 Computed properties:: 288 289 cov = covariance matrix [ inv(A'A); O(n^3) ] 290 var = coefficient variance [ diag(cov); O(n^2)] 291 std = standard deviation of coefficients [ sqrt(var); O(n^2) ] 292 p = test statistic for chisquare goodness of fit [ chi2.sf; O(1) ] 293 294 Methods:: 295 296 __call__(x): return the polynomial evaluated at x 297 ci(x,sigma=1): return confidence interval evaluated at x 298 pi(x,alpha=0.05): return prediction interval evaluated at x 299 300 Note that the covariance matrix will not include the ones column if 301 the polynomial goes through the origin. 302 """

303 - def __init__(self, s, origin=False):

304 self.origin = origin 305 self.coeff = N.ravel(s.x) 306 if origin: self.coeff = N.hstack((self.coeff,0)) 307 self.degree = len(self.coeff)-1 308 self.DoF = s.DoF 309 self.rnorm = s.rnorm 310 self._conf = s

311 - def _cov(self):

312 return self._conf.cov

313 - def _std(self):

314 return N.sqrt(self._var())

315 - def _var(self):

316 var = N.ravel(self._conf.var) 317 if self.origin: var = N.hstack((var,0)) 318 return var

319 - def _p(self):

320 return self._conf.p

321 cov = property(_cov,doc="covariance matrix") 322 var = property(_var,doc="result variance") 323 std = property(_std,doc="result standard deviation") 324 p = property(_p,doc="probability of rejection") 325 326

327 - def __call__(self, x):

328 """ 329 Evaluate the polynomial at x. 330 """ 331 return N.polyval(self.coeff,x)

332

333 - def ci(self, x, sigma=1):

334 """ 335 Evaluate the polynomial and the confidence intervals at x. 336 337 sigma=1 corresponds to a 1-sigma confidence interval 338 """ 339 A = _poly_matrix(x,self.degree,self.origin) 340 return self._conf.ci(A,sigma)

341

342 - def pi(self, x, p=0.05):

343 """ 344 Evaluate the polynomial and the prediction intervals at x. 345 346 p = 1-alpha = 0.05 corresponds to 95% prediction interval 347 """ 348 A = _poly_matrix(x,self.degree,self.origin) 349 return self._conf.pi(A,p)

350

351 - def __str__(self):

352 # TODO: better polynomial pretty printing using formatnum 353 return "Polynomial(%s)"%self.coeff

354

355 -def wpolyfit(x,y,dy=1,degree=None,origin=False):

356 """ 357 Return the polynomial of degree n that 358 minimizes sum( (p(x_i) - y_i)**2/dy_i**2). 359 360 if origin is True, the fit should go through the origin. 361 """ 362 assert degree != None, "Missing degree argument to wpolyfit" 363 364 A = _poly_matrix(x,degree,origin) 365 s = wsolve(A,y,dy) 366 return PolynomialModel(s,origin=origin)

367 368

369 -def demo():

370 import pylab 371 372 # Make fake data 373 x = N.linspace(-15,5,15) 374 th = N.polyval([.2,3,1,5],x) # polynomial 375 dy = N.sqrt(N.abs(th)) # poisson uncertainty estimate 376 y = N.random.normal(th,dy) # but normal generator 377 378 # Fit to a polynomial 379 poly = wpolyfit(x,y,dy=dy,degree=3) 380 381 # Plot the result 382 pylab.errorbar(x,y,yerr=dy,fmt='x') 383 pylab.hold(True) 384 px=N.linspace(x[0],x[-1],200) 385 py,pdy = poly.pi(px) 386 cy,cdy = poly.ci(px) 387 pylab.plot(px,py,'g-', 388 px,py+pdy,'g-.',px,py-pdy,'g-.', 389 px,cy+cdy,'r-.',px,cy-cdy,'r-.') 390 pylab.show()

391

392 -def test():

393 """ 394 smoke test...make sure the function continues to return the same 395 result for a particular system. 396 """ 397 x = N.array([0,1,2,3,4],'d') 398 y = N.array([ 2.5, 7.9, 13.9, 21.1, 44.4],'d') 399 dy = N.array([ 1.7, 2.4, 3.6, 4.8, 6.2],'d') 400 poly = wpolyfit(x,y,dy,1) 401 px = N.array([1.5],'d') 402 py,pi = poly.pi(px) 403 py,ci = poly.ci(px) 404 405 ## Uncomment these to show target values 406 #print "Tp = [%.16g, %.16g]"%(p[0],p[1]) 407 #print "Tdp = [%.16g, %.16g]"%(dp[0],dp[1]) 408 #print "Tpi,Tci = %.16g, %.16g"%(pi,ci) 409 Tp = N.array([7.787249069840737, 1.503992847461524]) 410 Tdp = N.array([1.522338103010216, 2.117633626902384]) 411 Tpi,Tci = 7.611128464981324, 2.342860389884832 412 413 perr = N.max(N.abs(poly.coeff-Tp)) 414 dperr = N.max(N.abs(poly.std-Tdp)) 415 cierr = N.abs(ci-Tci) 416 pierr = N.abs(pi-Tpi) 417 assert perr < 1e-15,"||p-Tp||=%g"%perr 418 assert dperr < 1e-15,"||dp-Tdp||=%g"%dperr 419 assert cierr < 1e-15,"||ci-Tci||=%g"%cierr 420 assert pierr < 1e-15,"||pi-Tpi||=%g"%pierr 421 assert py == poly(px),"direct call to poly function fails"

422 423 if __name__ == "__main__": 424 test() 425 # demo() 426

Source Code for Module reflectometry.reduction.wsolve