target_functions.chi2

1 ############################################################################### 2 # # 3 # Copyright (C) 2001-2004,2008 Edward d'Auvergne # 4 # Copyright (C) 2014 Troels E. Linnet # 5 # # 6 # This file is part of the program relax (http://www.nmr-relax.com). # 7 # # 8 # This program is free software: you can redistribute it and/or modify # 9 # it under the terms of the GNU General Public License as published by # 10 # the Free Software Foundation, either version 3 of the License, or # 11 # (at your option) any later version. # 12 # # 13 # This program is distributed in the hope that it will be useful, # 14 # but WITHOUT ANY WARRANTY; without even the implied warranty of # 15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # 16 # GNU General Public License for more details. # 17 # # 18 # You should have received a copy of the GNU General Public License # 19 # along with this program. If not, see <http://www.gnu.org/licenses/>. # 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 value for multi-dimensional axis. 70 #################### 71 72

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

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

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

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

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

158 159

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

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

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

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

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

254 255

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

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

311

Source Code for Module target_functions.chi2