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
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
112 grad = -2.0 * dot(1.0 / (errors**2) * (data - back_calc_vals), transpose(back_calc_grad))
113
114
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
156 return -2.0 * sum(1.0 / (errors**2) * (data - back_calc_vals) * back_calc_grad_j, axis=0)
157
158
159
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
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
259
260
261
262
263
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
269 return d2chi2
270