Help finding R^2 value for curves

40 views (last 30 days)
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
  2 Comments
Image Analyst
Image Analyst on 10 Jan 2021
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
BobbyJoe on 10 Jan 2021
The whole code is API2, it should run if you copy and paste into matlab. There is no other code

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 10 Jan 2021
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.)
  4 Comments
BobbyJoe
BobbyJoe on 11 Jan 2021
Thank you so so much star strider! extremely helpful as always!!!
Star Strider
Star Strider on 11 Jan 2021
As always, my pleasure!
Thank you!

Sign in to comment.

More Answers (0)

Tags

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!