Help finding R^2 value for curves

Hi all,

The following code produces a graph which compares simulated experimental results against actual experimental data points:

function API2

function C=kinetics(theta,t)
c0=[0.752;1.278;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(3,1);
    dcdt(1)=-theta(1).*c(1).^2+theta(2).*c(2)-theta(3).*c(1).*c(2);
    dcdt(2)=theta(1).*c(1).^1-theta(2).*c(2).^1-theta(3).*c(1).^1.*c(2);
    dcdt(3)=theta(3).*c(1).*c(2).^3;
    dC=dcdt;
    end
C=Cv;
end

T = [0 1 2 5 10 15];
t = T';

%y values for a
a_ydata = [0.752 0.0596 0.0596 0.0596 0.0502 0.0424];
A_Ydata = a_ydata';

%y values for b
b_ydata = [1.278 0.378 0.101 0.101 0.085 0.072];
B_Ydata = b_ydata';

c_ydata = [0 0.692 0.692 0.692 0.702 0.71];
C_Ydata = c_ydata';

c = [A_Ydata B_Ydata C_Ydata];
theta0=[0.5;0.5;0.5];

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure(1)
h = plot(t, c,'.');
set(h,{'Marker'},{'s';'d';'^'},{'MarkerFaceColor'},{'r';'b';'k'},{'MarkerEdgeColor'},{'r';'b';'k'});
hold on
hlp = plot(tv, Cfit,'LineWidth',1.5);
set(hlp,{'Color'},{'r';'b';'k'});
hold off
grid
xlabel('Time (min)')
ylabel('Concentration (M)')
legend(hlp, 'Rifamycin Oxazine', 'Piperazine', 'Rifampicin', 'Location','N')

end

I want to know what code should be added so that I can see the R^2 value for each curve to see how accurate the simulated curve is. Thanks

[Image Analyst]

You forgot to give us the code for API2 that sets the values for theta and t and actually calls

C=kinetics(theta,t)

Please do so, so we can run your code.

[BobbyJoe]

The whole code is API2, it should run if you copy and paste into matlab. There is no other code

Answer 1

[Star Strider]

I appreciate your quoting my code!

The value is not difficult to calculate.

One example using a much simpler single-variable curve fit (not the posted problem):

yfit = @(b,x) exp(b(1)) .* x.^b(2);                         % Power Function
SSECF = @(b) sum((yfit(b,x) - y).^2);                       % Sum-Squared-Error Cost Function
B = fminsearch(SSECF, [1; 1]);

ypred = yfit(B,x);
SSE=sum((y-ypred).^2);
SST=sum((y-mean(y)).^2);
Rsq   = 1 - (SSE/SST);

Here, the idea is to compute the value for the fit of this particular power function given by ‘yfit’ (where ‘x’ and ‘y’ are the data vectors) and the data, using the fminsearch function to do the nonlinear regression. Here ‘B’ is the coefficients vector comkputed by fminsearch.

Following that, the code calculates ‘ypred’ (the predicted value, corresponding to your model vector) and the actual data. The ‘SSE’ value is explained sum-of-squares, and ‘SST’ is the total sum-of-squares, as described in the Wikipedia article. .

So in the oiriginal problem, ‘ypred’ is the model, ‘y’ the data vector, and ‘Rsq’ the desired result.

In the context of the posted parameter estimation problem, replace ‘yfit’ wiith ‘Cfit’, and ‘y’ with ‘c’, however do this with each curve, not all of them simultaneously, since may not be defined for the sort of problem you are posing here. (I have never seen it defined for a multiple-fit problem such as this, however there may be ways to extend it to such problems.)

[BobbyJoe]

It's an excellent code that has been extremely helpful thus far!!

Thank you so much Star Strider for the very helpful answer! I tried to do as you followed and trying to repeat your code for each curve but I'm getting negative infinity or a value that is definitely not correct.. I'm not sure what I'm doing wrong :/ Do you mind showing me the code for one curve to see what I'm missing? Thank you so much!

[Star Strider]

Appending this code:

Cfit_mtx = kinetics(theta, t);                                      % Calculate R² For Each Compartment
for k = 1:size(Cfit,2)
    ypred  = Cfit_mtx(:,k);
    SSE    = sum((c(:,k)-ypred).^2);
    SST    = sum((c(:,k)-mean(c(:,k))).^2);
    Rsq(k) = 1 - (SSE/SST);
    fprintf('\t\tR² c(%d) = %7.4f\n',k, Rsq(k))
end

to the end (after the plots) of the ‘BobbyJo_GA.m’ code I posted earlier I get:

	Rate Constants:
		Theta(1) =  7.47489
		Theta(2) =  0.89174
		Theta(3) =  0.97884
		Theta(4) =  0.01652

        R² c(1) =  0.9161
        R² c(2) =  0.8991
        R² c(3) =  0.7007

all of which appear to be finite. The ‘Cfit_mtx’ values are evaluated at the independent variable ‘t’ points so that a direct calculation is possible.

If my Answer helped you solve your problem, please Accept it!

.

[BobbyJoe]

Thank you so so much star strider! extremely helpful as always!!!

[Star Strider]

As always, my pleasure!

Thank you!