maths_fns.chi2

1 ############################################################################### 2 # # 3 # Copyright (C) 2003, 2004, 2008 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 """Module containing functions for calculating the chi-squared value, gradient, and Hessian.""" 24 25 # Python module imports. 26 from numpy import dot, sum, transpose 27 28 29 # Chi-squared value. 30 #################### 31 32

33 -def chi2(data, back_calc_vals, errors):

34 """Function to calculate the chi-squared value. 35 36 The chi-squared equation 37 ======================== 38 39 The equation is:: 40 41 _n_ 42 \ (yi - yi(theta)) ** 2 43 chi^2(theta) = > --------------------- 44 /__ sigma_i ** 2 45 i=1 46 47 where 48 - i is the index over data sets. 49 - theta is the parameter vector. 50 - yi are the values of the measured data set. 51 - yi(theta) are the values of the back calculated data set. 52 - sigma_i are the values of the error set. 53 54 55 @param data: The vector of yi values. 56 @type data: numpy rank-1 size N array 57 @param back_calc_vals: The vector of yi(theta) values. 58 @type back_calc_vals: numpy rank-1 size N array 59 @param errors: The vector of sigma_i values. 60 @type errors: numpy rank-1 size N array 61 @return: The chi-squared value. 62 @rtype: float 63 """ 64 65 # Calculate the chi-squared statistic. 66 return sum((1.0 / errors * (data - back_calc_vals))**2, axis=0)

67 68 69 # Chi-squared gradient. 70 ####################### 71 72

73 -def dchi2(dchi2, M, data, back_calc_vals, back_calc_grad, errors):

74 """Calculate the full chi-squared gradient. 75 76 The chi-squared gradient 77 ======================== 78 79 The equation is:: 80 81 _n_ 82 dchi^2(theta) \ / yi - yi(theta) dyi(theta) \ 83 ------------- = -2 > | -------------- . ---------- | 84 dthetaj /__ \ sigma_i**2 dthetaj / 85 i=1 86 87 where 88 - i is the index over data sets. 89 - j is the parameter index of the gradient. 90 - theta is the parameter vector. 91 - yi are the values of the measured data set. 92 - yi(theta) are the values of the back calculated data set. 93 - dyi(theta)/dthetaj are the values of the back calculated gradient for parameter j. 94 - sigma_i are the values of the error set. 95 96 97 @param dchi2: The chi-squared gradient data structure to place the gradient elements 98 into. 99 @type dchi2: numpy rank-1 size M array 100 @param M: The dimensions of the gradient. 101 @type M: int 102 @param data: The vector of yi values. 103 @type data: numpy rank-1 size N array 104 @param back_calc_vals: The vector of yi(theta) values. 105 @type back_calc_vals: numpy rank-1 size N array 106 @param back_calc_grad: The matrix of dyi(theta)/dtheta values. 107 @type back_calc_grad: numpy rank-2 size MxN array 108 @param errors: The vector of sigma_i values. 109 @type errors: numpy rank-1 size N array 110 """ 111 112 # Calculate the chi-squared gradient. 113 grad = -2.0 * dot(1.0 / (errors**2) * (data - back_calc_vals), transpose(back_calc_grad)) 114 115 # Pack the elements. 116 for i in xrange(M): 117 dchi2[i] = grad[i]

118 119

120 -def dchi2_element(data, back_calc_vals, back_calc_grad_j, errors):

121 """Calculate the chi-squared gradient element j. 122 123 The chi-squared gradient 124 ======================== 125 126 The equation is:: 127 128 _n_ 129 dchi^2(theta) \ / yi - yi(theta) dyi(theta) \ 130 ------------- = -2 > | -------------- . ---------- | 131 dthetaj /__ \ sigma_i**2 dthetaj / 132 i=1 133 134 where 135 - i is the index over data sets. 136 - j is the parameter index of the gradient. 137 - theta is the parameter vector. 138 - yi are the values of the measured data set. 139 - yi(theta) are the values of the back calculated data set. 140 - dyi(theta)/dthetaj are the values of the back calculated gradient for parameter j. 141 - sigma_i are the values of the error set. 142 143 144 @param data: The vector of yi values. 145 @type data: numpy rank-1 size N array 146 @param back_calc_vals: The vector of yi(theta) values. 147 @type back_calc_vals: numpy rank-1 size N array 148 @param back_calc_grad_j: The vector of dyi(theta)/dthetaj values for parameter j. 149 @type back_calc_grad_j: numpy rank-1 size N array 150 @param errors: The vector of sigma_i values. 151 @type errors: numpy rank-1 size N array 152 @return: The chi-squared gradient element j. 153 @rtype: float 154 """ 155 156 # Calculate the chi-squared gradient. 157 return -2.0 * sum(1.0 / (errors**2) * (data - back_calc_vals) * back_calc_grad_j, axis=0)

158 159 160 # Chi-squared Hessian. 161 ###################### 162 163

164 -def d2chi2(d2chi2, M, data, back_calc_vals, back_calc_grad, back_calc_hess, errors):

165 """Calculate the full chi-squared Hessian. 166 167 The chi-squared Hessian 168 ======================= 169 170 The equation is:: 171 172 _n_ 173 d2chi^2(theta) \ 1 / dyi(theta) dyi(theta) d2yi(theta) \ 174 --------------- = 2 > ---------- | ---------- . ---------- - (yi-yi(theta)) . --------------- | 175 dthetaj.dthetak /__ sigma_i**2 \ dthetaj dthetak dthetaj.dthetak / 176 i=1 177 178 where 179 - i is the index over data sets. 180 - j is the parameter index for the first dimension of the Hessian. 181 - k is the parameter index for the second dimension of the Hessian. 182 - theta is the parameter vector. 183 - yi are the values of the measured data set. 184 - yi(theta) are the values of the back calculated data set. 185 - dyi(theta)/dthetaj are the values of the back calculated gradient for parameter j. 186 - d2yi(theta)/dthetaj.dthetak are the values of the back calculated Hessian for the 187 parameters j and k. 188 - sigma_i are the values of the error set. 189 190 191 @param d2chi2: The chi-squared Hessian data structure to place the Hessian elements 192 into. 193 @type d2chi2: numpy rank-2 size MxM array 194 @param M: The size of the first dimension of the Hessian. 195 @type M: int 196 @param data: The vector of yi values. 197 @type data: numpy rank-1 size N array 198 @param back_calc_vals: The vector of yi(theta) values. 199 @type back_calc_vals: numpy rank-1 size N array 200 @param back_calc_grad: The matrix of dyi(theta)/dtheta values. 201 @type back_calc_grad: numpy rank-2 size MxN array 202 @param back_calc_hess: The matrix of d2yi(theta)/dtheta.dtheta values. 203 @type back_calc_hess: numpy rank-3 size MxMxN array 204 @param errors: The vector of sigma_i values. 205 @type errors: numpy rank-1 size N array 206 """ 207 208 # Calculate the chi-squared Hessian. 209 for j in xrange(M): 210 for k in xrange(M): 211 d2chi2[j, k] = 0.0 212 for i in xrange(len(data)): 213 d2chi2[j, k] = d2chi2[j, k] + 2.0 / (errors[i]**2) * (back_calc_grad[j, i] * back_calc_grad[k, i] - (data[i] - back_calc_vals[i]) * back_calc_hess[j, k, i])

214 215

216 -def d2chi2_element(data, back_calc_vals, back_calc_grad_j, back_calc_grad_k, back_calc_hess_jk, errors):

217 """Calculate the chi-squared Hessian element {j, k}. 218 219 The chi-squared Hessian 220 ======================= 221 222 The equation is:: 223 224 _n_ 225 d2chi^2(theta) \ 1 / dyi(theta) dyi(theta) d2yi(theta) \ 226 --------------- = 2 > ---------- | ---------- . ---------- - (yi-yi(theta)) . --------------- | 227 dthetaj.dthetak /__ sigma_i**2 \ dthetaj dthetak dthetaj.dthetak / 228 i=1 229 230 where 231 - i is the index over data sets. 232 - j is the parameter index for the first dimension of the Hessian. 233 - k is the parameter index for the second dimension of the Hessian. 234 - theta is the parameter vector. 235 - yi are the values of the measured data set. 236 - yi(theta) are the values of the back calculated data set. 237 - dyi(theta)/dthetaj are the values of the back calculated gradient for parameter j. 238 - d2yi(theta)/dthetaj.dthetak are the values of the back calculated Hessian for the 239 parameters j and k. 240 - sigma_i are the values of the error set. 241 242 243 @param data: The vector of yi values. 244 @type data: numpy rank-1 size N array 245 @param back_calc_vals: The vector of yi(theta) values. 246 @type back_calc_vals: numpy rank-1 size N array 247 @param back_calc_grad_j: The vector of dyi(theta)/dthetaj values for parameter j. 248 @type back_calc_grad_j: numpy rank-1 size N array 249 @param back_calc_grad_k: The vector of dyi(theta)/dthetak values for parameter k. 250 @type back_calc_grad_k: numpy rank-1 size N array 251 @param back_calc_hess_jk: The vector of d2yi(theta)/dthetaj.dthetak values at {j, k}. 252 @type back_calc_hess_jk: numpy rank-1 size N array 253 @param errors: The vector of sigma_i values. 254 @type errors: numpy rank-1 size N array 255 @return: The chi-squared Hessian element {j,k}. 256 @rtype: float 257 """ 258 259 # Calculate the chi-squared Hessian. 260 #return 2.0 * sum(1.0 / (errors**2) * (back_calc_grad_j * back_calc_grad_k - (data - back_calc_vals) * back_calc_hess), axis=0) 261 #return 2.0 * sum((back_calc_grad_j * back_calc_grad_k - (data - back_calc_vals) * back_calc_hess) / errors**2, axis=0) 262 263 # Calculate the chi-squared Hessian. 264 # This is faster than the above sums, and having the errors term first appears to minimise roundoff errors. 265 d2chi2 = 0.0 266 for i in xrange(len(data)): 267 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_jk[i]) 268 269 # Return the {j, k} element. 270 return d2chi2

271

Source Code for Module maths_fns.chi2