% Class 18 - Regression day 2
% Coded by NFR on 10.24
%
% Warm up problem
% In class 18
%Coded by  Erin, Lance, and Sarah
%10.24.2019
close all
x = [1 2 3 4 5 6 7 8 9 10]';
y = [5 9.6 11.7 15.5 19.4 21 25.9 27.2 30.8 34.1]';
%you HAVE to put those stupid little thingys after the brackets
p = polyfit(x,y,1);
f1 = @(x) p(1)*x +p(2);
figure
plot(x,y,'ro')
hold on
fplot(f1,[1 10])
[fit2, gof]  = fit(x,y,'poly1');
confint(fit2,0.95)
plot(fit2,x,y,'predobs')
xlabel('X')
ylabel('Y')
title('Y vs X and line fit and 95% confidence bounds')
rsqr = gof.rsquare
%
% Example of using the fit command
% for a non-linear equation
x = [1 2 3 4 5];
y = [0.5 2 2.9 3.5 4];
% Reminder: plot if you don't know the form
% check to see if it is linear
%
ft = fittype('exp((y-b)/a)','independent','y','dependent','x');
% We have the form of the equation
% now we fit data to this fit equation
model2 = fit(y',x',ft,'StartPoint',[1 1],'Lower',[0 0],'Upper',[100 100])
model1 = fit(y',x',ft)
% When using fit, put your independent variable first
% followed by dependent
% NOTE: model 1 and 2 are the same.
hold off
% NOTE: parameters are given in alphabetical order
plot(model1,y,x,'predobs')
% If you wanted to use the model directly
% in this case a y is given and you want 
% to solve for x
xsol = feval(model1,3.2)
% There is also a grapical user interface (GUI)
% tool available called cftool.  You must
% first have your data available as saved
% variables.
%
% Cftool can be useful when Matlab
% is having a hard time finding parameters
% you can set bounds and visible see what 
% the fit looks like.
%
% Limitation: can only be used for x,y or x,y,z problems
%
% Example with a loop
% Lots of data, repetitive plotting + fitting
% sounds like a great time to use matlab
D = readmatrix('Slopes');
% First determine number of trials and time points
[ntp,nTrial] = size(D);
timevec = (1:ntp)*10; % minutes
% Set up a loop to go through each reactor
% Initialize a vector to save fill rates
FR = zeros(nTrial,1);
RSvec = zeros(nTrial,1); % to record Rsqr values
for i = 1:nTrial
    ydata = D(:,i);
    % Determine the slope, fit linear model
    p = polyfit(timevec',ydata,1);
    slope = p(1);
    FR(i,1) = slope;
    % Let's call the fit command to get r squared
    [fit1,gof] = fit(timevec',ydata,'poly1');
    RS = gof.rsquare;
    RSvec(i,1) = RS;
    % Plot the raw data with fit and display
    % rsquared value in the title of the plot
    subplot(4,5,i)
    plot(fit1,timevec',ydata)
    xlabel('Time(min)')
    ylabel('Volume (gal)')
    %title(['Rsqr = ',num2str(RS)])
    title(['Flow Rate = ',num2str(p(1)),'gal/min'])
end
% What is the average flow rate
AvgFlowRate = mean(FR)
StandardDev_FlowRate = std(FR)
% What about cases where you need to read
% from a bunch of different files
% 
% Example = Data1.xlsx, Data2.xlsx, etc.
% Have a loop that counts 1, 2, 3...
% within the loop, call the file name like this
% D = readmatrix(['Data',int2str(i),'.xlsx'])