1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23 """Module containing functions for calculating the chi-squared value, gradient, and Hessian."""
24
25
26 from numpy import dot, sum, transpose
27
28
29
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
66 return sum((1.0 / errors * (data - back_calc_vals))**2, axis=0)
67
68
69
70
71
72
73 -def dchi2(dchi2, M, data, back_calc_vals, back_calc_grad, errors):
74 """Calculate the full chi-squared gradient.
75
76 The chi-squared gradient
77 ========================
78
79 The equation is::
80
81 _n_
82 dchi^2(theta) \ / yi - yi(theta) dyi(theta) \
83 ------------- = -2 > | -------------- . ---------- |
84 dthetaj /__ \ sigma_i**2 dthetaj /
85 i=1
86
87 where
88 - i is the index over data sets.
89 - j is the parameter index of the gradient.
90 - theta is the parameter vector.
91 - yi are the values of the measured data set.
92 - yi(theta) are the values of the back calculated data set.
93 - dyi(theta)/dthetaj are the values of the back calculated gradient for parameter j.
94 - sigma_i are the values of the error set.
95
96
97 @param dchi2: The chi-squared gradient data structure to place the gradient elements
98 into.
99 @type dchi2: numpy rank-1 size M array
100 @param M: The dimensions of the gradient.
101 @type M: int
102 @param data: The vector of yi values.
103 @type data: numpy rank-1 size N array
104 @param back_calc_vals: The vector of yi(theta) values.
105 @type back_calc_vals: numpy rank-1 size N array
106 @param back_calc_grad: The matrix of dyi(theta)/dtheta values.
107 @type back_calc_grad: numpy rank-2 size MxN array
108 @param errors: The vector of sigma_i values.
109 @type errors: numpy rank-1 size N array
110 """
111
112
113 grad = -2.0 * dot(1.0 / (errors**2) * (data - back_calc_vals), transpose(back_calc_grad))
114
115
116 for i in xrange(M):
117 dchi2[i] = grad[i]
118
119
120 -def dchi2_element(data, back_calc_vals, back_calc_grad_j, errors):
121 """Calculate the chi-squared gradient element j.
122
123 The chi-squared gradient
124 ========================
125
126 The equation is::
127
128 _n_
129 dchi^2(theta) \ / yi - yi(theta) dyi(theta) \
130 ------------- = -2 > | -------------- . ---------- |
131 dthetaj /__ \ sigma_i**2 dthetaj /
132 i=1
133
134 where
135 - i is the index over data sets.
136 - j is the parameter index of the gradient.
137 - theta is the parameter vector.
138 - yi are the values of the measured data set.
139 - yi(theta) are the values of the back calculated data set.
140 - dyi(theta)/dthetaj are the values of the back calculated gradient for parameter j.
141 - sigma_i are the values of the error set.
142
143
144 @param data: The vector of yi values.
145 @type data: numpy rank-1 size N array
146 @param back_calc_vals: The vector of yi(theta) values.
147 @type back_calc_vals: numpy rank-1 size N array
148 @param back_calc_grad_j: The vector of dyi(theta)/dthetaj values for parameter j.
149 @type back_calc_grad_j: numpy rank-1 size N array
150 @param errors: The vector of sigma_i values.
151 @type errors: numpy rank-1 size N array
152 @return: The chi-squared gradient element j.
153 @rtype: float
154 """
155
156
157 return -2.0 * sum(1.0 / (errors**2) * (data - back_calc_vals) * back_calc_grad_j, axis=0)
158
159
160
161
162
163
164 -def d2chi2(d2chi2, M, data, back_calc_vals, back_calc_grad, back_calc_hess, errors):
165 """Calculate the full chi-squared Hessian.
166
167 The chi-squared Hessian
168 =======================
169
170 The equation is::
171
172 _n_
173 d2chi^2(theta) \ 1 / dyi(theta) dyi(theta) d2yi(theta) \
174 --------------- = 2 > ---------- | ---------- . ---------- - (yi-yi(theta)) . --------------- |
175 dthetaj.dthetak /__ sigma_i**2 \ dthetaj dthetak dthetaj.dthetak /
176 i=1
177
178 where
179 - i is the index over data sets.
180 - j is the parameter index for the first dimension of the Hessian.
181 - k is the parameter index for the second dimension of the Hessian.
182 - theta is the parameter vector.
183 - yi are the values of the measured data set.
184 - yi(theta) are the values of the back calculated data set.
185 - dyi(theta)/dthetaj are the values of the back calculated gradient for parameter j.
186 - d2yi(theta)/dthetaj.dthetak are the values of the back calculated Hessian for the
187 parameters j and k.
188 - sigma_i are the values of the error set.
189
190
191 @param d2chi2: The chi-squared Hessian data structure to place the Hessian elements
192 into.
193 @type d2chi2: numpy rank-2 size MxM array
194 @param M: The size of the first dimension of the Hessian.
195 @type M: int
196 @param data: The vector of yi values.
197 @type data: numpy rank-1 size N array
198 @param back_calc_vals: The vector of yi(theta) values.
199 @type back_calc_vals: numpy rank-1 size N array
200 @param back_calc_grad: The matrix of dyi(theta)/dtheta values.
201 @type back_calc_grad: numpy rank-2 size MxN array
202 @param back_calc_hess: The matrix of d2yi(theta)/dtheta.dtheta values.
203 @type back_calc_hess: numpy rank-3 size MxMxN array
204 @param errors: The vector of sigma_i values.
205 @type errors: numpy rank-1 size N array
206 """
207
208
209 for j in xrange(M):
210 for k in xrange(M):
211 d2chi2[j, k] = 0.0
212 for i in xrange(len(data)):
213 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])
214
215
216 -def d2chi2_element(data, back_calc_vals, back_calc_grad_j, back_calc_grad_k, back_calc_hess_jk, errors):
217 """Calculate the chi-squared Hessian element {j, k}.
218
219 The chi-squared Hessian
220 =======================
221
222 The equation is::
223
224 _n_
225 d2chi^2(theta) \ 1 / dyi(theta) dyi(theta) d2yi(theta) \
226 --------------- = 2 > ---------- | ---------- . ---------- - (yi-yi(theta)) . --------------- |
227 dthetaj.dthetak /__ sigma_i**2 \ dthetaj dthetak dthetaj.dthetak /
228 i=1
229
230 where
231 - i is the index over data sets.
232 - j is the parameter index for the first dimension of the Hessian.
233 - k is the parameter index for the second dimension of the Hessian.
234 - theta is the parameter vector.
235 - yi are the values of the measured data set.
236 - yi(theta) are the values of the back calculated data set.
237 - dyi(theta)/dthetaj are the values of the back calculated gradient for parameter j.
238 - d2yi(theta)/dthetaj.dthetak are the values of the back calculated Hessian for the
239 parameters j and k.
240 - sigma_i are the values of the error set.
241
242
243 @param data: The vector of yi values.
244 @type data: numpy rank-1 size N array
245 @param back_calc_vals: The vector of yi(theta) values.
246 @type back_calc_vals: numpy rank-1 size N array
247 @param back_calc_grad_j: The vector of dyi(theta)/dthetaj values for parameter j.
248 @type back_calc_grad_j: numpy rank-1 size N array
249 @param back_calc_grad_k: The vector of dyi(theta)/dthetak values for parameter k.
250 @type back_calc_grad_k: numpy rank-1 size N array
251 @param back_calc_hess_jk: The vector of d2yi(theta)/dthetaj.dthetak values at {j, k}.
252 @type back_calc_hess_jk: numpy rank-1 size N array
253 @param errors: The vector of sigma_i values.
254 @type errors: numpy rank-1 size N array
255 @return: The chi-squared Hessian element {j,k}.
256 @rtype: float
257 """
258
259
260
261
262
263
264
265 d2chi2 = 0.0
266 for i in xrange(len(data)):
267 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])
268
269
270 return d2chi2
271