target_functions.chi2

1 ############################################################################### 2 # # 3 # Copyright (C) 2003-2012 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 gradient. 69 ####################### 70 71

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

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

117 118

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

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

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

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

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

213 214

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

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

270

Source Code for Module target_functions.chi2