target_functions.chi2

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2014 Edward d'Auvergne # 4 # # 5 # This file is part of the program relax (http://www.nmr-relax.com). # 6 # # 7 # This program 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 3 of the License, or # 10 # (at your option) any later version. # 11 # # 12 # This program 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 this program. If not, see <http://www.gnu.org/licenses/>. # 19 # # 20 ############################################################################### 21 22 """Module containing functions for calculating the chi-squared value, gradient, and Hessian.""" 23 24 # Python module imports. 25 from numpy import dot, sum, transpose 26 27 28 # Chi-squared value. 29 #################### 30 31

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

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

66 67 68 # Chi-squared value for multi-dimensional axis. 69 #################### 70 71

72 -def chi2_rankN(data, back_calc_vals, errors):

73 """Function to calculate the chi-squared value for multiple numpy array axis. 74 75 The chi-squared equation 76 ======================== 77 78 The equation is:: 79 80 _n_ 81 \ (yi - yi(theta)) ** 2 82 chi^2(theta) = > --------------------- 83 /__ sigma_i ** 2 84 i=1 85 86 where 87 - i is the index over data sets. 88 - theta is the parameter vector. 89 - yi are the values of the measured data set. 90 - yi(theta) are the values of the back calculated data set. 91 - sigma_i are the values of the error set. 92 93 94 @param data: The multi dimensional vectors of yi values. 95 @type data: numpy multi dimensional array 96 @param back_calc_vals: The multi dimensional vectors of yi(theta) values. 97 @type back_calc_vals: numpy multi dimensional array 98 @param errors: The multi dimensional vectors of sigma_i values. 99 @type errors: numpy multi dimensional array 100 @return: The chi-squared value. 101 @rtype: float 102 """ 103 104 # Calculate the chi-squared statistic. 105 return sum((1.0 / errors * (data - back_calc_vals))**2)

106 107 108 # Chi-squared gradient. 109 ####################### 110 111

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

113 """Calculate the full chi-squared gradient. 114 115 The chi-squared gradient 116 ======================== 117 118 The equation is:: 119 120 _n_ 121 dchi^2(theta) \ / yi - yi(theta) dyi(theta) \ 122 ------------- = -2 > | -------------- . ---------- | 123 dthetaj /__ \ sigma_i**2 dthetaj / 124 i=1 125 126 where 127 - i is the index over data sets. 128 - j is the parameter index of the gradient. 129 - theta is the parameter vector. 130 - yi are the values of the measured data set. 131 - yi(theta) are the values of the back calculated data set. 132 - dyi(theta)/dthetaj are the values of the back calculated gradient for parameter j. 133 - sigma_i are the values of the error set. 134 135 136 @param dchi2: The chi-squared gradient data structure to place the gradient elements 137 into. 138 @type dchi2: numpy rank-1 size M array 139 @param M: The dimensions of the gradient. 140 @type M: int 141 @param data: The vector of yi values. 142 @type data: numpy rank-1 size N array 143 @param back_calc_vals: The vector of yi(theta) values. 144 @type back_calc_vals: numpy rank-1 size N array 145 @param back_calc_grad: The matrix of dyi(theta)/dtheta values. 146 @type back_calc_grad: numpy rank-2 size MxN array 147 @param errors: The vector of sigma_i values. 148 @type errors: numpy rank-1 size N array 149 """ 150 151 # Calculate the chi-squared gradient. 152 grad = -2.0 * dot(1.0 / (errors**2) * (data - back_calc_vals), transpose(back_calc_grad)) 153 154 # Pack the elements. 155 for i in range(M): 156 dchi2[i] = grad[i]

157 158

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

160 """Calculate the chi-squared gradient element j. 161 162 The chi-squared gradient 163 ======================== 164 165 The equation is:: 166 167 _n_ 168 dchi^2(theta) \ / yi - yi(theta) dyi(theta) \ 169 ------------- = -2 > | -------------- . ---------- | 170 dthetaj /__ \ sigma_i**2 dthetaj / 171 i=1 172 173 where 174 - i is the index over data sets. 175 - j is the parameter index of the gradient. 176 - theta is the parameter vector. 177 - yi are the values of the measured data set. 178 - yi(theta) are the values of the back calculated data set. 179 - dyi(theta)/dthetaj are the values of the back calculated gradient for parameter j. 180 - sigma_i are the values of the error set. 181 182 183 @param data: The vector of yi values. 184 @type data: numpy rank-1 size N array 185 @param back_calc_vals: The vector of yi(theta) values. 186 @type back_calc_vals: numpy rank-1 size N array 187 @param back_calc_grad_j: The vector of dyi(theta)/dthetaj values for parameter j. 188 @type back_calc_grad_j: numpy rank-1 size N array 189 @param errors: The vector of sigma_i values. 190 @type errors: numpy rank-1 size N array 191 @return: The chi-squared gradient element j. 192 @rtype: float 193 """ 194 195 # Calculate the chi-squared gradient. 196 return -2.0 * sum(1.0 / (errors**2) * (data - back_calc_vals) * back_calc_grad_j, axis=0)

197 198 199 # Chi-squared Hessian. 200 ###################### 201 202

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

204 """Calculate the full chi-squared Hessian. 205 206 The chi-squared Hessian 207 ======================= 208 209 The equation is:: 210 211 _n_ 212 d2chi^2(theta) \ 1 / dyi(theta) dyi(theta) d2yi(theta) \ 213 --------------- = 2 > ---------- | ---------- . ---------- - (yi-yi(theta)) . --------------- | 214 dthetaj.dthetak /__ sigma_i**2 \ dthetaj dthetak dthetaj.dthetak / 215 i=1 216 217 where 218 - i is the index over data sets. 219 - j is the parameter index for the first dimension of the Hessian. 220 - k is the parameter index for the second dimension of the Hessian. 221 - theta is the parameter vector. 222 - yi are the values of the measured data set. 223 - yi(theta) are the values of the back calculated data set. 224 - dyi(theta)/dthetaj are the values of the back calculated gradient for parameter j. 225 - d2yi(theta)/dthetaj.dthetak are the values of the back calculated Hessian for the 226 parameters j and k. 227 - sigma_i are the values of the error set. 228 229 230 @param d2chi2: The chi-squared Hessian data structure to place the Hessian elements 231 into. 232 @type d2chi2: numpy rank-2 size MxM array 233 @param M: The size of the first dimension of the Hessian. 234 @type M: int 235 @param data: The vector of yi values. 236 @type data: numpy rank-1 size N array 237 @param back_calc_vals: The vector of yi(theta) values. 238 @type back_calc_vals: numpy rank-1 size N array 239 @param back_calc_grad: The matrix of dyi(theta)/dtheta values. 240 @type back_calc_grad: numpy rank-2 size MxN array 241 @param back_calc_hess: The matrix of d2yi(theta)/dtheta.dtheta values. 242 @type back_calc_hess: numpy rank-3 size MxMxN array 243 @param errors: The vector of sigma_i values. 244 @type errors: numpy rank-1 size N array 245 """ 246 247 # Calculate the chi-squared Hessian. 248 for j in range(M): 249 for k in range(M): 250 d2chi2[j, k] = 0.0 251 for i in range(len(data)): 252 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])

253 254

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

256 """Calculate the chi-squared Hessian element {j, k}. 257 258 The chi-squared Hessian 259 ======================= 260 261 The equation is:: 262 263 _n_ 264 d2chi^2(theta) \ 1 / dyi(theta) dyi(theta) d2yi(theta) \ 265 --------------- = 2 > ---------- | ---------- . ---------- - (yi-yi(theta)) . --------------- | 266 dthetaj.dthetak /__ sigma_i**2 \ dthetaj dthetak dthetaj.dthetak / 267 i=1 268 269 where 270 - i is the index over data sets. 271 - j is the parameter index for the first dimension of the Hessian. 272 - k is the parameter index for the second dimension of the Hessian. 273 - theta is the parameter vector. 274 - yi are the values of the measured data set. 275 - yi(theta) are the values of the back calculated data set. 276 - dyi(theta)/dthetaj are the values of the back calculated gradient for parameter j. 277 - d2yi(theta)/dthetaj.dthetak are the values of the back calculated Hessian for the 278 parameters j and k. 279 - sigma_i are the values of the error set. 280 281 282 @param data: The vector of yi values. 283 @type data: numpy rank-1 size N array 284 @param back_calc_vals: The vector of yi(theta) values. 285 @type back_calc_vals: numpy rank-1 size N array 286 @param back_calc_grad_j: The vector of dyi(theta)/dthetaj values for parameter j. 287 @type back_calc_grad_j: numpy rank-1 size N array 288 @param back_calc_grad_k: The vector of dyi(theta)/dthetak values for parameter k. 289 @type back_calc_grad_k: numpy rank-1 size N array 290 @param back_calc_hess_jk: The vector of d2yi(theta)/dthetaj.dthetak values at {j, k}. 291 @type back_calc_hess_jk: numpy rank-1 size N array 292 @param errors: The vector of sigma_i values. 293 @type errors: numpy rank-1 size N array 294 @return: The chi-squared Hessian element {j,k}. 295 @rtype: float 296 """ 297 298 # Calculate the chi-squared Hessian. 299 #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) 300 #return 2.0 * sum((back_calc_grad_j * back_calc_grad_k - (data - back_calc_vals) * back_calc_hess) / errors**2, axis=0) 301 302 # Calculate the chi-squared Hessian. 303 # This is faster than the above sums, and having the errors term first appears to minimise roundoff errors. 304 d2chi2 = 0.0 305 for i in range(len(data)): 306 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]) 307 308 # Return the {j, k} element. 309 return d2chi2

310

Source Code for Module target_functions.chi2