% Class 17 - Regression
% Coded by NFR on 10.22.2019
%
% Warm Up
xsol = fsolve(@SF,[1,1,1]) % FROM SPAM group
% Example of linear regression with polyfit
D = readmatrix('660nmData_reg');
plot(D(:,1),D(:,2),'.','MarkerSize',20)
xlabel('Protein Concentration (ug/ml)')
ylabel('ABS 660nm')
p = polyfit(D(:,1),D(:,2),1)
% P1 = slope, p2 = intercept
F1 = @(x) p(1)*x+p(2);
hold on
fplot(F1,[0 2000])
hold off
% Reminder that you can use polyval to figure
% out an unknown value
% Another (more robust) built in function
% is fit command, you can also get goodness
% of fit (gof) data from this.
[Fit2,gof] = fit(D(:,1),D(:,2),'poly1')
% Cool, we can see it in the command prompt
% what can we do with these within the code?
% 1) Plotting a fit directly
hold off
figure
plot(Fit2,D(:,1),D(:,2),'predobs')
% 2) Pull off the R square value
rsqr = gof.rsquare;
% 3) We can also look at custom confidence intervals
confint(Fit2,0.999)
% This prints out confidence interval for fit parameters
% in this case slope and intercept.
% 
% Be careful transforming data (linearizing)
% it can misrepresent your residuals
% best to use a nonlinear model directly
% we will talk about this on Thursday
%
% BONUS - how to make animations with matlab
% 
% Bonus time - countour plot as function of time
k = 23;
H = 100;
L = 200;
%t = 10;
Xvec = linspace(0,200,1000);
Yvec = linspace(0,100,1000);
[xmat, ymat] = meshgrid(Xvec,Yvec);
% Create a loop for multiple times
tvec = 1:20;
% If you want to collect the frames into a video
v = VideoWriter('TestVid.avi');
open(v);

for i = tvec
    F1 = @(x,y)   1600/(pi^2)*sin(pi*x./L)...
    .*sin(pi.*y/(2.*H))*exp(-(1/(4.*H.^2)...
    +1./L.^2).*k.*pi^2.*i);
    zmat = F1(xmat,ymat);
    contourf(Xvec,Yvec,zmat,'ShowText','on')
    %pause(.01)
    frame = getframe(gcf);
    writeVideo(v,frame);
end
close(v)


function F = SF(x)
% x is an array of the input variables
% x(1) = x
% x(2) = y
% x(3) = z
F(1) = 3*x(1).*x(2)+x(2)-x(3)-12;
F(2) = x(1) + x(2).*x(1).^2+x(3)-12;
F(3) = x(1)-x(2)-x(3)+4;
end