%% This demonstrates Hypothesis Testing using the chi-squared statistic

%This also demostrates the use of 'polyfit', which is a linear
%least-squares fit (see HW 1)


%% Create a data set

A=3.5;
B=7;
C=.1;

N=50;
sigma=10; %this is the sigma in normally distributed noise model

x=(1:N)';
y=A+B*x+C*x.^2; 
y_noise=y+sigma*randn(size(y));

ebars=sigma*ones(size(y));
figure
hold on
plot(x,y,'k','linewidth',2)
errorbar(x,y_noise,ebars,'r.','linewidth',1)
legend('underlying model','noisy data')
xlabel('X (m)','fontsize',16)
ylabel('Y (m)','fontsize',16)
set(gca,'fontsize',14)

%% Fit the data to a line 

n=1; % the polynomial order to use in the fit. n=1 is a line. 
P = polyfit(x,y_noise,n)

B_found=P(1);
A_found=P(2);

y_found_1=A_found+B_found*x; %this can also be done using polyval, but this is more explicit

%look at the fit
figure
hold on
plot(x,y_found_1,'k','linewidth',2)
errorbar(x,y_noise,ebars,'r.','linewidth',1)
legend('first order polynomial fit','noisy data')
xlabel('X (m)','fontsize',16)
ylabel('Y (m)','fontsize',16)
set(gca,'fontsize',14)


%% Calculate the p-value

%Null Hypothesis: A first order polynomical with the coefficeints found 
%above is the correct underlying model for the data.

%calculate the chi-squared test statistic
X2=sum(1/sigma^2*(y_found_1-y_noise).^2)

%The degrees of freedom are the number of data points minus the number of
%fit parameters

k=N-(n+1)

%the statistics toolbox has these functions, but here we'll just make our
%own

%this is the integral of the chi-squared PDF from 0 to x.
X2_CDF=inline('gammainc(x/2,k/2)','k','x')

%The p-value is the probability that X2 would be found to be more extreme
%than that of our data, assuming that our model (Null Hypothesis) is
%correct

pvalue=1-X2_CDF(k,X2)

siglevel=.01;
if pvalue<siglevel
    fprintf('The null hypothesis is rejected at the %g level\n',siglevel)
    fprintf('Found p-value: %g \n',pvalue)
else
    fprintf('Null hypothesis not rejected at the %g level\n',siglevel)
    fprintf('Found p-value: %f \n',pvalue)
end




%% Try again with a n=2 polynomial


n=2; % the polynomial order to use in the fit. n=1 is a line. 
P = polyfit(x,y_noise,n)

C_found=P(1);
B_found=P(2);
A_found=P(3);

y_found_1=A_found+B_found*x+C_found*x.^2; %this can also be done using polyval, but this is more explicit
%look at the fit

figure
hold on
plot(x,y_found_1,'k','linewidth',2)
errorbar(x,y_noise,ebars,'r.','linewidth',1)
legend('second order polynomial fit','noisy data')
xlabel('X (m)','fontsize',16)
ylabel('Y (m)','fontsize',16)
set(gca,'fontsize',14)

% Calculate the p-value

%Null Hypothesis: A first order polynomical with the coefficeints found 
%above is the correct underlying model for the data.

%calculate the chi-squared test statistic
X2=sum(1/sigma^2*(y_found_1-y_noise).^2)

%The degrees of freedom are the number of data points minus the number of
%fit parameters

k=N-(n+1)

%The p-value is the probability that X2 would be found to be more extreme
%than that of our data, assuming that our model (Null Hypothesis) is
%correct

pvalue=1-X2_CDF(k,X2)

siglevel=.01;
if pvalue<siglevel
    fprintf('The null hypothesis is rejected at the %g level\n',siglevel)
    fprintf('Found p-value: %g \n',pvalue)
else
    fprintf('Null hypothesis not rejected at the %g level\n',siglevel)
    fprintf('Found p-value: %g \n',pvalue)
end














