Source Code for Module maths_fns.chi2

  1  ############################################################################### 
  2  #                                                                             # 
  3  # Copyright (C) 2003, 2004 Edward d'Auvergne                                  # 
  4  #                                                                             # 
  5  # This file is part of the program relax.                                     # 
  6  #                                                                             # 
  7  # relax is free software; you can redistribute it and/or modify               # 
  8  # it under the terms of the GNU General Public License as published by        # 
  9  # the Free Software Foundation; either version 2 of the License, or           # 
 10  # (at your option) any later version.                                         # 
 11  #                                                                             # 
 12  # relax is distributed in the hope that it will be useful,                    # 
 13  # but WITHOUT ANY WARRANTY; without even the implied warranty of              # 
 14  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               # 
 15  # GNU General Public License for more details.                                # 
 16  #                                                                             # 
 17  # You should have received a copy of the GNU General Public License           # 
 18  # along with relax; if not, write to the Free Software                        # 
 19  # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA   # 
 20  #                                                                             # 
 21  ############################################################################### 
 22   
 23   
 24  from Numeric import Float64, sum, transpose, zeros 
 25   
 26   
 27  # Chi-squared value. 
 28  #################### 
 29   
 30 -def chi2(data, back_calc_vals, errors): 
 31      """Function to calculate the chi-squared value. 
 32   
 33      The chi-sqared equation 
 34      ~~~~~~~~~~~~~~~~~~~~~~~ 
 35              _n_ 
 36              \    (yi - yi()) ** 2 
 37      Chi2  =  >   ---------------- 
 38              /__    sigma_i ** 2 
 39              i=1 
 40   
 41      where: 
 42          yi are the values of the measured data set. 
 43          yi() are the values of the back calculated data set. 
 44          sigma_i are the values of the error set. 
 45   
 46      The chi-squared value is returned. 
 47      """ 
 48   
 49      # Calculate the chi-squared statistic. 
 50      return sum((1.0 / errors * (data - back_calc_vals))**2) 
 51   
 52   
 53  # Chi-squared gradient. 
 54  ####################### 
 55   
 56 -def dchi2(data, back_calc_vals, back_calc_grad, errors): 
 57      """Function to create the chi-squared gradient. 
 58   
 59      The chi-sqared gradient 
 60      ~~~~~~~~~~~~~~~~~~~~~~~ 
 61                     _n_ 
 62       dChi2         \   /  yi - yi()      dyi()  \  
 63      -------  =  -2  >  | ----------  .  ------- | 
 64      dthetaj        /__ \ sigma_i**2     dthetaj / 
 65                     i=1 
 66   
 67      where: 
 68          yi are the values of the measured data set. 
 69          yi() are the values of the back calculated data set. 
 70          sigma_i are the values of the error set. 
 71   
 72      The chi-squared gradient vector is returned. 
 73      """ 
 74   
 75      # Calculate the chi-squared gradient. 
 76      return -2.0 * sum(1.0 / (errors**2) * (data - back_calc_vals) * back_calc_grad) 
 77   
 78   
 79  # Chi-squared Hessian. 
 80  ###################### 
 81   
 82 -def d2chi2(data, back_calc_vals, back_calc_grad_j, back_calc_grad_k, back_calc_hess, errors): 
 83      """Function to create the chi-squared Hessian. 
 84   
 85      The chi-squared Hessian 
 86      ~~~~~~~~~~~~~~~~~~~~~~~ 
 87                            _n_ 
 88           d2chi2           \       1      /  dyi()     dyi()                         d2yi()     \  
 89      ---------------  =  2  >  ---------- | ------- . -------  -  (yi - yi()) . --------------- | 
 90      dthetaj.dthetak       /__ sigma_i**2 \ dthetaj   dthetak                   dthetaj.dthetak / 
 91                            i=1 
 92   
 93      where: 
 94          yi are the values of the measured relaxation data set. 
 95          yi() are the values of the back calculated relaxation data set. 
 96          sigma_i are the values of the error set. 
 97      """ 
 98   
 99      # Calculate the chi-squared Hessian. 
100      #return 2.0 * sum(1.0 / (errors**2) * (back_calc_grad_j * back_calc_grad_k - (data - back_calc_vals) * back_calc_hess)) 
101      #return 2.0 * sum((back_calc_grad_j * back_calc_grad_k - (data - back_calc_vals) * back_calc_hess) / errors**2) 
102   
103      # Calculate the chi-squared Hessian. 
104      # This is faster than the above sums, and having the errors term first appears to minimise roundoff errors. 
105      d2chi2 = 0.0 
106      for i in xrange(len(data)): 
107          d2chi2 = d2chi2 + 2.0 / (errors[i]**2) * (back_calc_grad_j[i] * back_calc_grad_k[i] - (data[i] - back_calc_vals[i]) * back_calc_hess[i]) 
108   
109      return d2chi2 
110