User:Trettels:Session 20 - 25MAR
Pre-class Activities
To Calculate
1.
Show that if <math>y(\mathit{b}|\mathbf{x})</math> is linear then minimizing <math>\chi^2</math> produces a set of linear equations for <math>b</math>.
First we will define a function <math>S</math> such that
<math> S_i=\frac{y_i - \sum_k b_k X_k(x_i)}{\sigma_i} </math> which gives us
<math> \chi^2 = \sum_i S_i^2 </math>
Now, in order to minimize <math>\chi^2</math> we need to find it's derivative:
<math> \frac{\delta\chi^2}{\delta b} = \sum_i S_i^2 \frac{\delta}{\delta S} \cdot \frac{\delta S}{\delta b} </math> <math>
= 2 \sum_i S_i\cdot \frac{\delta S}{\delta b}
</math> <math>
= 2 \sum_i \frac{y_i - \sum_k b_k X_k(x_i)}{\sigma_i}\cdot\frac{\delta S}{\delta b}
</math>
As <math>S_i</math> is linear in <math>b</math> we can see that <math> \frac{\delta S}{\delta b_k} = \frac{-X_k(x_i)}{\sigma_i}</math>, leaving us with
<math> \frac{\delta\chi^2}{\delta b} = 2 \sum_i \frac{y_i - \sum_k b_k X_k(x_i)}{\sigma_i}\cdot\frac{-X_k(x_i)}{\sigma_i}
= 2\sum_{i,k} \frac{-y_i X_k(x_i) + b_k X_k^2(x_i)}{\sigma_i^2}
</math>
Finally, to find the minima we set the derivative equal to zero and solve for <math>b</math>:
<math> 2\sum_{i,k} \frac{y_i X_k(x_i)}{\sigma_i^2} = 2\sum_{i,k} \frac{b_k X_k^2(x_i)}{\sigma_i^2} </math>
and taking <math> \mathbf{r_i} = \frac{1}{\sigma_i^2}</math>
<math> 2\mathbf{r}\mathbf{X^T}\mathbf{y} = 2\mathbf{r}\mathbf{X^T}\mathbf{X}\mathbf{b} </math>
and by multiplying both sides by <math>2\mathbf{r^{-1}}</math>:
<math> \mathbf{X^Ty} = \mathbf{X^tXb} </math>
2.
Given that <math>\mathbf{b}=\left[ \begin{array}{c}b_0 \\ b_1\end{array}\right]</math> and <math>\mathbf{x} = \left[\begin{array}{c} 1 \\ x_i \end{array}\right]</math> and isolating <math>\mathbf{b}</math> then we can see that
<math> \mathbf{b} =\mathbf{X^{-1}(X^T)^{-1}X^Ty} = \mathbf{(X)^{-1}y} </math>
so:
<math> \left[\begin{array}{c} b_0 \\ b_1\end{array}\right] = \left[\begin{array}{c} 1 \\ \frac{1}{x_i}\end{array}\right] y_i </math>
In Class
MATLAB code (done with Lori): Computes the linear, parabolic, exponential and full-quadratic models for the given set of data.
clear all;
imp=importdata('Chisqfitdata.txt')
x=imp.data(:,1);
y=imp.data(:,2);
sigma=imp.data(:,3);
errorbar(x,y,sigma,'*');
hold on;
init=[0.01, 2.8];
%quadratic
ymodel=@(x,b) (b(1) + b(2)*x.^2);
chisqfun = @(b) sum(((ymodel(x,b)-y)./sigma).^2);
b_hat = fminsearch(chisqfun, init);
b_hat
fprintf('\nQuad X^2: %f \n',chisqfun(b_hat));
plot(x,ymodel(x,b_hat),'r')
%linear
ymodel=@(x,b) (b(1) + b(2)*x);
chisqfun = @(b) sum(((ymodel(x,b)-y)./sigma).^2);
b_hat_1 = fminsearch(chisqfun,init);
b_hat_1
fprintf('\nLin X^2: %f\n', chisqfun(b_hat_1));
plot(x,ymodel(x,b_hat_1),'g')
%exponential
ymodel=@(x,b) (b(1) *exp( b(2)*x));
chisqfun = @(b) sum(((ymodel(x,b)-y)./sigma).^2);
b_hat_exp = fminsearch(chisqfun,init);
b_hat_exp
fprintf('\nExp X^2: %f\n', chisqfun(b_hat_exp));
plot(x,ymodel(x,b_hat_exp),'k');
%quad2
init=[0.01, 1, 2.8];
ymodel=@(x,b) (b(1) + b(2)*x + b(3)*x.^2);
chisqfun = @(b) sum(((ymodel(x,b)-y)./sigma).^2);
b_hat_q2 = fminsearch(chisqfun, init);
b_hat_q2
fprintf('\nQuad2 X^2: %f \n',chisqfun(b_hat_q2));
plot(x,ymodel(x,b_hat_q2),'c')
[b_hat_q2_2,~,~,~,~,hessian] = fminunc(chisqfun,init);
b_hat_q2_2
fprintf('\nQuad2 X^2: %f \n',chisqfun(b_hat_q2_2));
stderrormat_1=sqrt(diag(inv(0.5*hessian)))
plot(x,ymodel(x,b_hat_q2_2),'m')
%cubic_test
init=[0.01, 1, 2.8, 1];
ymodel=@(x,b) (b(1) + b(2)*x + b(3)*x.^2 + b(4)*x.^3);
chisqfun = @(b) sum(((ymodel(x,b)-y)./sigma).^2);
b_hat_cube = fminsearch(chisqfun, init);
b_hat_cube
fprintf('\nQuad2 X^2: %f \n',chisqfun(b_hat_cube));
plot(x,ymodel(x,b_hat_cube),'-cx')
h = 0.1;
unit = @(i) (1:4) == i;
hess = zeros(4,4);
for i=1:4, for j=1:4,
bpp = b_hat_cube + h*(unit(i)+unit(j));
bmm = b_hat_cube + h*(-unit(i)-unit(j));
bpm = b_hat_cube + h*(unit(i)-unit(j));
bmp = b_hat_cube + h*(-unit(i)+unit(j));
hess(i,j) = (chisqfun(bpp)+chisqfun(bmm)-chisqfun(bpm)-chisqfun(bmp))./(2*h)^2;
end
end
stderrmat_2 = sqrt(diag(inv(0.5*hess)))
legend('data','quadratic','linear', 'exponential','quadratic 2','quad2-fminunc', 'cubic', 'Location','NorthWest');
hold off;
