001/* 002 * ModelData -- A utility class used to model experimental data. 003 * 004 * Copyright (C) 2002-2006 by Joseph A. Huwaldt 005 * All rights reserved. 006 * 007 * This library is free software; you can redistribute it and/or 008 * modify it under the terms of the GNU Lesser General Public 009 * License as published by the Free Software Foundation; either 010 * version 2 of the License, or (at your option) any later version. 011 * 012 * This library is distributed in the hope that it will be useful, 013 * but WITHOUT ANY WARRANTY; without even the implied warranty of 014 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 015 * Lesser General Public License for more details. 016 * 017 * You should have received a copy of the GNU Lesser General Public License 018 * along with this program; if not, write to the Free Software 019 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 020 * Or visit: http://www.gnu.org/licenses/lgpl.html 021 */ 022package jahuwaldt.tools.math; 023 024import java.util.Arrays; 025 026/** 027 * Methods in this class are used to create functions that model experimental data. 028 * 029 * <p> Modified by: Joseph A. Huwaldt </p> 030 * 031 * @author Joseph A. Huwaldt Date: October 21, 2002 032 * @version February 13, 2014 033 */ 034public class ModelData { 035 036 // Tolerance used in single-value-decomposition to remove singular values. 037 private static final double TOL = 1E-13; 038 039 /** 040 * Method that returns the coefficients of a polynomial of the specified degree that 041 * best models or "fits" the supplied data values in a minimization of chi-squared 042 * sense. This method assumes that each data point is normally distributed and has a 043 * standard deviation of 1.0. 044 * 045 * @param xarr Array of independent parameter data values to fit. 046 * @param yarr Array of dependent data values (associated with each X value) to be modeled. 047 * @param coef An existing array that will be filled in with the coefficients of the 048 * polynomial, in decreasing power that best fits the sample data. Example: p(x) = A + 049 * B*x + C*x^2 + D*x^3 + E*x^4 corresponds to coef[] = {A, B, C, D, E}. The number of 050 * elements in the coef array determines the order of the polynomial created. Number 051 * of elements be greater than 1. 052 * @return The chi-squared value of the polynomial returned. A value near zero 053 * indicates a perfect fit. 054 */ 055 public static double polynomial(double[] xarr, double[] yarr, double[] coef) throws RootException { 056 if (xarr == null) 057 throw new NullPointerException("xarr == null"); 058 double[] sig = new double[xarr.length]; 059 Arrays.fill(sig, 1.0); 060 BasisFunction bf = new PolynomialFit(); 061 return fit(xarr, yarr, sig, bf, coef); 062 } 063 064 /** 065 * Method that returns the coefficients of a polynomial of the specified degree that 066 * best models or "fits" the supplied data values in a minimization of chi-squared 067 * sense. 068 * 069 * @param xarr Array of independent parameter data values to fit. 070 * @param yarr Array of dependent data values (associated with each X value) to be modeled. 071 * @param sig Array of individual standard deviations for each data point. 072 * @param coef An existing array that will be filled in with the coefficients of the 073 * polynomial, in decreasing power that best fits the sample data. Example: p(x) = A + 074 * B*x + C*x^2 + D*x^3 + E*x^4 corresponds to coef[] = {A, B, C, D, E}. The number of 075 * elements in the coef array determines the order of the polynomial created. Number 076 * of elements be greater than 1. 077 * @return The chi-squared value of the polynomial returned. A value near zero 078 * indicates a perfect fit. 079 */ 080 public static double polynomial(double[] xarr, double[] yarr, double[] sig, double[] coef) throws RootException { 081 BasisFunction bf = new PolynomialFit(); 082 return fit(xarr, yarr, sig, bf, coef); 083 } 084 085 /** 086 * Method that returns the coefficients of an arbitrary basis function that best 087 * models or "fits" the supplied data values in a minimization of chi-squared sense. 088 * 089 * @param xarr Array of independent parameter data values to fit. 090 * @param yarr Array of dependent data values (associated with each X value) to be 091 * modeled. 092 * @param sig Array of individual standard deviations for each data point. 093 * @param func The basis function used to generate the data fit. 094 * @param coef An existing array that will be filled in with the coefficients of the 095 * basis function that best fits the data. 096 * @return The chi-squared value of the function is returned. A value near zero 097 * indicates a perfect fit. 098 */ 099 public static double fit(double[] xarr, double[] yarr, double[] sig, BasisFunction func, double[] coef) throws RootException { 100 if (xarr == null) 101 throw new NullPointerException("xarr == null"); 102 if (yarr == null) 103 throw new NullPointerException("yarr == null"); 104 if (sig == null) 105 throw new NullPointerException("sig == null"); 106 if (func == null) 107 throw new NullPointerException("func == null"); 108 if (coef == null) 109 throw new NullPointerException("coef == null"); 110 111 int ma = coef.length; 112 if (ma < func.getMinNumCoef()) 113 throw new IllegalArgumentException("Number of coefficients must be greater than " + (func.getMinNumCoef() - 1)); 114 115 int ndata = xarr.length; 116 if (ndata != yarr.length) 117 throw new IllegalArgumentException("xarr and yarr must be the same length."); 118 if (ndata != sig.length) 119 throw new IllegalArgumentException("sig must be the same length as xarr and yarr."); 120 121 // Allocate the arrays that we are going to need. 122 double[][] u = new double[ndata][ma]; 123 double[][] v = new double[ma][ma]; 124 double[] w = new double[ma]; 125 126 // Use Singular Value Decomposition to find the parameters of the basis function. 127 double chisq = svdfit(xarr, yarr, sig, coef, u, v, w, func); 128 129 return chisq; 130 } 131 132 /** 133 * Returns the slope of the line formed by a linear regression through the specified 134 * data arrays. 135 * 136 * @param x The independent array of data points. 137 * @param y The dependent array of data points (must be the same size as the x array. 138 * @return The slope of the line formed by a linear regression through the x and y 139 * arrays of data points. 140 */ 141 public static double linearSlope(double[] x, double[] y) { 142 143 // Find the mean of the input arrays. 144 int length = x.length; 145 double sumX = 0, sumY = 0; 146 for (int i = 0; i < length; ++i) { 147 sumX += x[i]; 148 sumY += y[i]; 149 } 150 double xbar = sumX / length; 151 double ybar = sumY / length; 152 153 // b = sum( (x - xbar)*(y - ybar) ) / sum ( (x - xbar)^2 ) 154 double denom = 0; 155 double num = 0; 156 for (int i = 0; i < 4; ++i) { 157 double a = x[i] - xbar; 158 num += a * (y[i] - ybar); 159 denom += a * a; 160 } 161 double b = num / denom; 162 163 return b; 164 } 165 166 /** 167 * Used to test the methods in this class. 168 */ 169 public static void main(String[] args) { 170 171 System.out.println(); 172 System.out.println("Testing ModelData..."); 173 174 try { 175 System.out.println(" Poly Curve Fit #1: linear"); 176 double[] x = {0, 1, 2, 3, 4, 5}; 177 double[] y = {1, 3, 5, 7, 9, 11}; 178 double[] coef = new double[2]; 179 double chisq = polynomial(x, y, coef); 180 System.out.println(" coef = " + (float)coef[0] + ", " + (float)coef[1] + ", chisq = " + (float)chisq); 181 System.out.println(" should be: p(x) = 1 + 2*x\n"); 182 183 System.out.println(" Poly Curve Fit #2: quadratic"); 184 double[] x1 = {3, 4, 5}; 185 double[] y1 = {-13.87, -0.09, -13.95}; 186 double[] coef1 = new double[3]; 187 chisq = polynomial(x1, y1, coef1); 188 System.out.println(" coef = " + (float)coef1[0] + ", " + (float)coef1[1] + ", " + (float)coef1[2] + ", chisq = " + (float)chisq); 189 System.out.println(" should be: p(x) = -221.05 + 110.52*x - 13.82*x^2"); 190 191 } catch (Exception e) { 192 e.printStackTrace(); 193 } 194 195 System.out.println("Done."); 196 197 } 198 199 //----------------------------------------------------------------------------------- 200 /** 201 * Given a set of data points x,y with individual standard deviations sig, use chi^2 202 * minimization to determine the coefficients "a" of the supplied basis fitting 203 * function y = sumi(ai*funci(x)). 204 */ 205 private static double svdfit(double[] x, double[] y, double[] sig, double[] a, double[][] u, double[][] v, double[] w, 206 BasisFunction f) throws RootException { 207 int ndata = x.length; 208 int ma = a.length; 209 210 double[] beta = new double[ndata]; 211 double[] afunc = new double[ma]; 212 213 // Accumulate coefficients of the fitting matrix. 214 for (int i = 0; i < ndata; ++i) { 215 f.func(x[i], afunc); 216 double sigi = 1.0 / sig[i]; 217 for (int j = 0; j < ma; ++j) { 218 u[i][j] = afunc[j] * sigi; 219 } 220 beta[i] = y[i] * sigi; 221 } 222 223 // Singular value decomposition. 224 svdcmp(u, w, v); 225 226 // Edit the singular values. 227 double wmax = 0; 228 for (int j = 0; j < ma; ++j) { 229 if (w[j] > wmax) 230 wmax = w[j]; 231 } 232 double thresh = TOL * wmax; 233 for (int j = 0; j < ma; ++j) { 234 if (w[j] < thresh) 235 w[j] = 0; 236 } 237 238 svdksb(u, w, v, beta, a); 239 240 // Evaluate chi-square. 241 double chisq = 0; 242 for (int i = 0; i < ndata; ++i) { 243 f.func(x[i], afunc); 244 double sum = 0; 245 for (int j = 0; j < ma; ++j) { 246 sum += a[j] * afunc[j]; 247 } 248 double tmp = (y[i] - sum) / sig[i]; 249 chisq += tmp * tmp; 250 } 251 252 return chisq; 253 } 254 255 /** 256 * Solves A*X = B for a vector X, where A is specified by the arrays u[][], w[], and 257 * v[][] as returned by svdcmp(). b is the right-hand side. x is the output solution 258 * vector. No input quantities are destroyed, so the routine can be called 259 * sequentially with different b's. 260 */ 261 private static void svdksb(double[][] u, double[] w, double[][] v, double b[], double x[]) { 262 int n = u[0].length; 263 int m = u.length; 264 265 double[] tmp = new double[n]; 266 for (int j = 0; j < n; ++j) { 267 // Calculate U^T*B. 268 double s = 0; 269 if (w[j] != 0) { 270 // Nonzero result only if w[j] is nonzero. 271 for (int i = 0; i < m; ++i) { 272 s += u[i][j] * b[i]; 273 } 274 s /= w[j]; 275 } 276 tmp[j] = s; 277 } 278 279 // Matrix multiply by V to get answer. 280 for (int j = 0; j < n; ++j) { 281 double s = 0; 282 for (int jj = 0; jj < n; ++jj) { 283 s += v[j][jj] * tmp[jj]; 284 } 285 x[j] = s; 286 } 287 } 288 289 /** 290 * Given a matrix a[][], this routine computes its singular value decomposition. A = 291 * U*W*V^T. The matrix U replaces "a" on output. The diagonal matrix of singular 292 * values W is output as a vector w[]. The matrix V (not the transpose V^T) is output 293 * as v[][]. 294 */ 295 private static void svdcmp(double[][] a, double[] w, double[][] v) throws RootException { 296 int n = a[0].length; 297 int m = a.length; 298 int k, l = 0; 299 300 double[] rv1 = new double[n]; 301 double g = 0; 302 double scale = 0; 303 double anorm = 0; 304 double s = 0; 305 306 // Householder reduction to bidiagonal form. 307 for (int i = 0; i < n; ++i) { 308 l = i + 1; 309 rv1[i] = scale * g; 310 s = g = scale = 0; 311 312 if (i < m) { 313 for (k = i; k < m; ++k) { 314 scale += Math.abs(a[k][i]); 315 } 316 if (scale != 0) { 317 for (k = i; k < m; ++k) { 318 a[k][i] /= scale; 319 s += a[k][i] * a[k][i]; 320 } 321 double f = a[i][i]; 322 g = -MathTools.sign(Math.sqrt(s), f); 323 double h = f * g - s; 324 a[i][i] = f - g; 325 for (int j = l; j < n; ++j) { 326 for (s = 0, k = i; k < m; ++k) { 327 s += a[k][i] * a[k][j]; 328 } 329 f = s / h; 330 for (k = i; k < m; ++k) { 331 a[k][j] += f * a[k][i]; 332 } 333 } 334 for (k = i; k < m; ++k) { 335 a[k][i] *= scale; 336 } 337 } 338 339 } 340 341 w[i] = scale * g; 342 g = s = scale = 0; 343 344 if (i < m && i + 1 != n) { 345 for (k = l; k < n; ++k) { 346 scale += Math.abs(a[i][k]); 347 } 348 if (scale != 0) { 349 for (k = l; k < n; ++k) { 350 a[i][k] /= scale; 351 s += a[i][k] * a[i][k]; 352 } 353 double f = a[i][l]; 354 g = -MathTools.sign(Math.sqrt(s), f); 355 double h = f * g - s; 356 a[i][l] = f - g; 357 for (k = l; k < n; ++k) { 358 rv1[k] = a[i][k] / h; 359 } 360 for (int j = l; j < m; ++j) { 361 for (s = 0, k = l; k < n; ++k) { 362 s += a[j][k] * a[i][k]; 363 } 364 for (k = l; k < n; ++k) { 365 a[j][k] += s * rv1[k]; 366 } 367 } 368 for (k = l; k < n; ++k) { 369 a[i][k] *= scale; 370 } 371 } 372 373 } 374 375 anorm = Math.max(anorm, (Math.abs(w[i]) + Math.abs(rv1[i]))); 376 } // Next i 377 378 // Accumulation of right-hand transformations. 379 for (int i = n - 1; i >= 0; --i) { 380 if (i + 1 < n) { 381 if (g != 0) { 382 // Double division to avoid possible underflow. 383 for (int j = l; j < n; ++j) { 384 v[j][i] = (a[i][j] / a[i][l]) / g; 385 } 386 for (int j = l; j < n; ++j) { 387 for (s = 0, k = l; k < n; ++k) { 388 s += a[i][k] * v[k][j]; 389 } 390 for (k = l; k < n; ++k) { 391 v[k][j] += s * v[k][i]; 392 } 393 } 394 } 395 for (int j = l; j < n; ++j) { 396 v[i][j] = v[j][i] = 0; 397 } 398 } 399 400 v[i][i] = 1; 401 g = rv1[i]; 402 l = i; 403 } // Next i 404 405 // Accumulation of left-hand transformations. 406 for (int i = Math.min(m, n) - 1; i >= 0; --i) { 407 l = i + 1; 408 g = w[i]; 409 for (int j = l; j < n; ++j) { 410 a[i][j] = 0; 411 } 412 if (g != 0) { 413 g = 1.0 / g; 414 for (int j = l; j < n; ++j) { 415 for (s = 0, k = l; k < m; ++k) { 416 s += a[k][i] * a[k][j]; 417 } 418 double f = (s / a[i][i]) * g; 419 for (k = i; k < m; ++k) { 420 a[k][j] += f * a[k][i]; 421 } 422 } 423 for (int j = i; j < m; ++j) { 424 a[j][i] *= g; 425 } 426 427 } else 428 for (int j = i; j < m; ++j) { 429 a[j][i] = 0; 430 } 431 432 ++a[i][i]; 433 } // Next i 434 435 // Diagonalization of the bidiagonal form: Loop over singular values, 436 // and over allowed iterations. 437 for (k = n - 1; k >= 0; --k) { 438 for (int its = 0; its < 30; ++its) { 439 boolean flag = true; 440 double x, y, z, c, f, h; 441 int nm = 0; 442 443 // Test for splitting. 444 for (l = k; l >= 0; --l) { 445 nm = l - 1; // Note that rv1[1] is always zero. 446 if ((Math.abs(rv1[l]) + anorm) == anorm) { 447 flag = false; 448 break; 449 } 450 if ((Math.abs(w[nm]) + anorm) == anorm) 451 break; 452 } 453 454 // Cancelation of rv1[l], if l > 1. 455 if (flag) { 456 c = 0; 457 s = 1; 458 for (int i = l; i <= k; ++i) { 459 f = s * rv1[i]; 460 rv1[i] = c * rv1[i]; 461 if ((Math.abs(f) + anorm) == anorm) 462 break; 463 g = w[i]; 464 h = pythag(f, g); 465 w[i] = h; 466 h = 1 / h; 467 c = g * h; 468 s = -f * h; 469 for (int j = 0; j < m; ++j) { 470 y = a[j][nm]; 471 z = a[j][i]; 472 a[j][nm] = y * c + z * s; 473 a[j][i] = z * c - y * s; 474 } 475 } 476 } 477 478 z = w[k]; 479 if (l == k) { 480 // Convergence 481 // Singular value is made nonnegative. 482 if (z < 0) { 483 w[k] = -z; 484 for (int j = 0; j < n; ++j) { 485 v[j][k] = -v[j][k]; 486 } 487 } 488 break; 489 } 490 491 if (its == 30) 492 throw new RootException("No convergence in 30 svdcmp() iterations."); 493 494 // Shift from bottom 2-by-2 minor. 495 x = w[l]; 496 nm = k - 1; 497 y = w[nm]; 498 g = rv1[nm]; 499 h = rv1[k]; 500 f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2 * h * y); 501 g = pythag(f, 1); 502 f = ((x - z) * (x + z) + h * ((y / (f + MathTools.sign(g, f))) - h)) / x; 503 c = s = 1; 504 505 // Next QR transformation. 506 for (int j = l; j <= nm; ++j) { 507 int i = j + 1; 508 g = rv1[i]; 509 y = w[i]; 510 h = s * g; 511 g = c * g; 512 z = pythag(f, h); 513 rv1[j] = z; 514 c = f / z; 515 s = h / z; 516 f = x * c + g * s; 517 g = g * c - x * s; 518 h = y * s; 519 y *= c; 520 for (int jj = 0; jj < n; ++jj) { 521 x = v[jj][j]; 522 z = v[jj][i]; 523 v[jj][j] = x * c + z * s; 524 v[jj][i] = z * c - x * s; 525 } 526 z = pythag(f, h); 527 w[j] = z; 528 529 // Rotation can be arbitrary if z == 0. 530 if (z != 0) { 531 z = 1 / z; 532 c = f * z; 533 s = h * z; 534 } 535 f = c * g + s * y; 536 x = c * y - s * g; 537 for (int jj = 0; jj < m; ++jj) { 538 y = a[jj][j]; 539 z = a[jj][i]; 540 a[jj][j] = y * c + z * s; 541 a[jj][i] = z * c - y * s; 542 } 543 544 } // Next j 545 546 rv1[l] = 0; 547 rv1[k] = f; 548 w[k] = x; 549 550 } // Next its 551 } // Next k 552 553 } 554 555 /** 556 * Method that calculates sqrt(a^2 + b^2) without destructive underflow or overflow. 557 */ 558 private static double pythag(double a, double b) { 559 double absa = Math.abs(a); 560 double absb = Math.abs(b); 561 562 if (absa > absb) { 563 double ratio = absb / absa; 564 return absa * Math.sqrt(1 + ratio * ratio); 565 566 } else { 567 if (absb == 0) 568 return 0.0; 569 else { 570 double ratio = absa / absb; 571 return absb * Math.sqrt(1 + ratio * ratio); 572 } 573 } 574 } 575}