# User:Trettels:Session 20 - 25MAR

Jump to navigation Jump to search

# Pre-class Activities

## To Calculate

### 1.

Show that if $y(\mathit{b}|\mathbf{x})$ is linear then minimizing $\chi^2$ produces a set of linear equations for $b$.

First we will define a function $S$ such that

$S_i=\frac{y_i - \sum_k b_k X_k(x_i)}{\sigma_i}$ which gives us

$\chi^2 = \sum_i S_i^2$

Now, in order to minimize $\chi^2$ we need to find it's derivative:

$\frac{\delta\chi^2}{\delta b} = \sum_i S_i^2 \frac{\delta}{\delta S} \cdot \frac{\delta S}{\delta b}$ $= 2 \sum_i S_i\cdot \frac{\delta S}{\delta b} $ $= 2 \sum_i \frac{y_i - \sum_k b_k X_k(x_i)}{\sigma_i}\cdot\frac{\delta S}{\delta b} $

As $S_i$ is linear in $b$ we can see that $\frac{\delta S}{\delta b_k} = \frac{-X_k(x_i)}{\sigma_i}$, leaving us with

$\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} $

Finally, to find the minima we set the derivative equal to zero and solve for $b$:

$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}$

and taking $\mathbf{r_i} = \frac{1}{\sigma_i^2}$

$2\mathbf{r}\mathbf{X^T}\mathbf{y} = 2\mathbf{r}\mathbf{X^T}\mathbf{X}\mathbf{b}$

and by multiplying both sides by $2\mathbf{r^{-1}}$:

$\mathbf{X^Ty} = \mathbf{X^tXb}$

### 2.

Given that $\mathbf{b}=\left[ \begin{array}{c}b_0 \\ b_1\end{array}\right]$ and $\mathbf{x} = \left[\begin{array}{c} 1 \\ x_i \end{array}\right]$ and isolating $\mathbf{b}$ then we can see that

$\mathbf{b} =\mathbf{X^{-1}(X^T)^{-1}X^Ty} = \mathbf{(X)^{-1}y}$

so:

$\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$

# 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;