1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22 """Module containing functions for calculating the chi-squared value, gradient, and Hessian."""
23
24
25 from numpy import dot, sum, transpose
26
27
28
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
65 return sum((1.0 / errors * (data - back_calc_vals))**2, axis=0)
66
67
68
69
70
71
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
105 return sum((1.0 / errors * (data - back_calc_vals))**2)
106
107
108
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
152 grad = -2.0 * dot(1.0 / (errors**2) * (data - back_calc_vals), transpose(back_calc_grad))
153
154
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
196 return -2.0 * sum(1.0 / (errors**2) * (data - back_calc_vals) * back_calc_grad_j, axis=0)
197
198
199
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
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
299
300
301
302
303
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
309 return d2chi2
310