% Modified version of the MATLAB template file FitNonlinSR.m used to solve
% Example 2 of Supplemental Unit 3 of "A First Course on Kinetics and
% Reaction Engineering."
%
function S4_Example_1(p_guess)
    % Set variables
    x = [1	0
    1.1	1
    1.2	2
    1.3	3
    1.4	4
    1.5	5
    1.6	6
    1.7	7
    1.8	8
    1.9	9
    2	10];
    % response variables
    y_hat = [7.8
    10.3
    12.1
    15.0
    16.4
    18.2
    20.0
    23.7
    24.6
    28.3
    29.6];
    
    % Function that calculates the responses using the model
    function y = nlmodel(p,x)
        n_data = size(x,1);
        y = zeros(n_data,1);
        for i = 1:n_data
            % p(i) is the value of parameter i
            % x(i,j) is the value of set variable xj for data point i
            y(i) = (p(1)*x(i,1)+1)*sqrt(p(2)*x(i,2)+p(3));
        end
    end % of internal function nlmodel

    % Perform the regression analysis
    [pf,resid,J,sigma,mse] = nlinfit(x,y_hat,@nlmodel,p_guess);
    confint = nlparci(pf,resid,'covar',sigma);
    pf_u = (confint(:,2) - confint(:,1))/2.0;
    
    % Calculate the model-predicted responses
    y = nlmodel(pf,x);
    
    % Display the correlation coefficient, parameters and uncertainties
    format shortE;
    r_squared = 1.0 - (sum(resid.^2))/(sum((y_hat-mean(y_hat)).^2))
    display('Best Values for the Parameters:')
    pf
    display('95% Confidence Intervals for the Parameters:')
    pf_u

    % Generate the graphical output
    n_set = size(x,2);
    if n_set == 1
        % Model plot
        % Use many data points for the model so the curve is smooth
        xx = linspace(min(x),max(x))';
        yy = nlmodel(pf,xx);
        figure;
        plot(x,y,'r',x,y_hat,'o')
        xlabel('X')
        ylabel('Y')
        title('Model Plot')
        legend('Model','Data')
    else
        % Parity plot
        figure;
        plot(y_hat,y_hat,'r',y_hat,y,'o')
        xlabel('Measured Response')
        ylabel('Predicted Response')
        title('Parity Plot')

        % Residuals plots for each set variable
        for n = 1:n_set
            figure;
            plot(x(:,n),zeros(length(y),1),'r',x(:,n),resid,'o')
            xlabel(['x',num2str(n)])
            ylabel('Residual')
            title(['Residuals Plot vs. x',num2str(n)])
        end
    end
end % of S4_Example_1.m