Curve fitting with an integral involved

55 views (last 30 days)
I have a 2 parameter model that includes an integral function and I´m quite lost how to solve this.
My data consists of a set of given molecule sizes, r_m, and a calculated response, K.
I want to fit this K with a theoretical model, where each K_i comes from integrating a particle size distribiution function, f(r), for instance a Gaussian distribution. The equations then look like this:
as you can see the integral function has the parameteres that I want to fit, but also the lower limit of the integral in the numerator involves a variable, r_m. I´m not sure if I can implement this situation to one of Matlab´s built-in fitting procedures, or do I have to code my own fitting algorithm? Here´s something I tried but I know I´m just playing Frankenstein here:
% Data = ...
[0.5 1
1.1 0.83
1.6 0.74
2.2 0.55
2.5 0.28
3.5 0];
r = Data(:,1);
K_exp = Data(:,2);
F = @(x,xdata) quad( (exp(-1/2.*((xdata - x(1))/x(2)).^2).*(1 - (xdata(1)/xdata).^2)),xdata(1),120)./quad(exp(-1/2.*((xdata - x(1))/x(2)).^2,0,120) ; ;
x0 = [6 0.5] ;
[x,resnorm,~,exitflag,output] = lsqcurvefit(F,x0,r,K_exp)
Any help on how to approach this problem is greatly appreciated!

Accepted Answer

Are Mjaavatten
Are Mjaavatten on 17 Jul 2018
Edited: Are Mjaavatten on 17 Jul 2018
function theta = Gil
data = [
0.5 1
1.1 0.83
1.6 0.74
2.2 0.55
2.5 0.28
3.5 0];
xdata = data(:,1);
ydata = data(:,2);
% Inital guess for parameters:
rp0 = 1;
rs0 = 1;
theta0 = [rp0;rs0];
% lsqcurvefit is in the optimization toolbox.
% fit, from the curve fitting toolbox may be an alternative
theta = lsqcurvefit(@Kdvec,theta0,xdata,ydata);
% Check fit:
figure;
plot(xdata,ydata,'or')
hold on
xplot = linspace(min(xdata),max(xdata));
plot(xplot,Kdvec(theta,xplot))
hold off
end
function Kvec = Kdvec(theta,xdata)
% Vector of Kd for every xdata:
Kvec = zeros(size(xdata));
for i = 1:length(xdata)
Kvec(i) = Kd(theta,xdata(i));
end
end
function K = Kd(theta,rm)
% Calculates integral. Assumes integrand(1000) is negligible
rp = theta(1);
sp = theta(2);
gauss = @(r) exp(-0.5*((r-rp)/sp).^2);
integrand = @(r) gauss(r).*(1-(rm./r).^2);
K = integral(integrand,rm,1000)/integral(gauss,0,1000);
end
  18 Comments
Are Mjaavatten
Are Mjaavatten on 28 Jun 2021
Edited: Are Mjaavatten on 28 Jun 2021
Gianrico: Take a look at the anwer by Matt J to Calculate Uncertainty for fitted parameter from least squares fit. Also browse through the discussion in the comments to that answer.
Are Mjaavatten
Are Mjaavatten on 30 Jun 2021
You want to find the uncertainty in the parameters found. Thinking a bit more closely about the expression by Matt J, I now realise that this may not be so relevant, since it is independent of the ydata values and says little about the goodness of fit. I fear that the short answer to your question may be that a closed form expression for the uncertainty does not exist.
Your question is an interesting one, though. I recommend that you formulate it as a separate question to Matlab Answers, in the hope that some clever expert may have some insight to offer.
One possible approach may be a Monte Carlo method where you introduce random errors to ydata and find some statistics from that.

Sign in to comment.

More Answers (5)

Costas Papadopoulos
Costas Papadopoulos on 2 Sep 2019
Hi,
I am trying to modify the code provided by Are Mjaavatten to apply weighted Langevin fitting on M(H)/Ms vs H data. The model I try to fit is:
weighted_langevin.JPG
where Fv(μ) a lognormal distribution:
Fv(μ) = 1/(μ*σ*(2*pi)^(1/2))*exp(-ln(μ/μ0)^2/(2*σ^2))
I actually used the same code but changed the variable names and the equations but I get several erros. Thanks in advance for any help
function MH = MomentDist
load ('magnetization.txt', '-ascii');
xdata = magnetization(:,1);
ydata = magnetization(:,2);
% Inital guess for parameters:
mi0 = 8e-14;
sigma0 = 2;
MH0 = [mi0;sigma0];
% lsqcurvefit is in the optimization toolbox.
% fit, from the curve fitting toolbox may be an alternative
MH = lsqcurvefit(@MMSvec,MH0,xdata,ydata);
% Check fit:
figure;
plot(xdata,ydata,'or')
hold on
xplot = linspace(min(xdata),max(xdata));
plot(xplot,MMSvec(MH,xplot))
hold off
end
function MMvec = MMSvec(MH,xdata)
% Vector of Kd for every xdata:
MMvec = zeros(size(xdata));
for i = 1:length(xdata)
MMvec(i) = MMS(MH,xdata(i));
end
end
function WLM = MMS(MH,mu)
% Calculates integral. Assumes integrand(1000) is negligible
mi = MH(1);
sigma = MH(2);
lognormal = @(mu) 1/(mu.*sigma.*(2*pi).^(1/2)).*exp((-(log(mu)-log(mi)).^2/(2.*sigma.^2)));
integrand = @(mu) lognormal(mu).*(coth(mu*xdata)-(1/(mu*xdata)));
WLM = integral(integrand,0,1000);
end
  6 Comments
Costas Papadopoulos
Costas Papadopoulos on 2 Sep 2019
Yes, I think that the problem is within the Langevin function. Pleas may you look at the model that I attatched in my question and suggest a possible right modification? Thanks

Sign in to comment.


Costas Papadopoulos
Costas Papadopoulos on 5 Sep 2019
Edited: Costas Papadopoulos on 5 Sep 2019
Hi again the code now runs but the fitting is far from good. It is more as if a single Langevin fitting has been applied. Any Ideas? Thanks in advance.
function MH = MomentDist
load ('magnetization.txt', '-ascii');
xdata = magnetization(:,1);
ydata = magnetization(:,2);
% Inital guess for parameters:
mi0 = 0.003;
sigma0 = 2;
MH0 = [mi0;sigma0];
% lsqcurvefit is in the optimization toolbox.
% fit, from the curve fitting toolbox may be an alternative
MH = lsqcurvefit(@MMSvector,MH0,xdata,ydata);
% Check fit:
figure;
plot(xdata,ydata,'or')
hold on
xplot = linspace(min(xdata),max(xdata));
plot(xplot,MMSvector(MH,xplot))
hold off
end
function MMSvec = MMSvector(MH,xdata)
% Vector of MMS for every xdata:
MMSvec = zeros(size(xdata));
for i = 1:length(xdata)
MMSvec(i) = MMS(MH,xdata(i));
end
end
function WLM = MMS(MH, x)
% Calculates integral. Assumes integrand(1000) is negligible
mi = MH(1);
sigma = MH(2);
lognormal = @(m) (m*sigma*(2*pi).^(1/2)).*exp(-(log(m/mi).^2/(2*sigma.^2)));
integrand = @(m) lognormal(m).*(coth(x)-(1/(x)));
WLM = integral(integrand,0,1000);
end
  2 Comments
Torsten
Torsten on 5 Sep 2019
Edited: Torsten on 5 Sep 2019
As written, your integral evaluates to
WLM = coth(x) - 1/x
since the integral over the lognormal density function from m=0 to m=Inf is equal to 1.
Thus WLM does not depend on the fitting parameters mi and sigma at all.
I doubt that this is what you want.
Costas Papadopoulos
Costas Papadopoulos on 5 Sep 2019
How should I write it to depend on distribution parameters? x=μΗ/kT but H/T is the xdata axis so I think it shoulf be x=μ/k however this doesn't work.

Sign in to comment.


Christina MacAskill
Christina MacAskill on 15 Feb 2021
Edited: Christina MacAskill on 15 Feb 2021
Hi,
I am trying to modify the code that @Are Mjaavatten posted to fit the following equation:
Here, I need to fit the equation for Ktrans,TOI and ve,TOI. The following parameters have the following set values:
t = 0:1:40;
T = 40;
Ktrans,RR = 0.1;
Ve,RR = 0.1;
R10,RR = 0.7;
R10,TOI = 0.55;
R = Ktrans,TOI/Ktrans,RR;
R1,RR (t) = [0.700000000000000 0.943337115297993 1.10117152595811 1.19865999687980 1.25353183962319 1.27843925419341 1.28249973305886 1.27232729639958 1.25273560941394 1.22722692370035 1.19833857557102 1.16789282643654 1.13717974605949 1.10709276431873 1.07823012430127 1.05097135016971 1.02553514095478 1.00202329296918 0.980454017107267 0.960787153232550 0.942943166951052 0.926817364737644 0.912290430224259 0.899236133648157 0.887526875270124 0.877037576393890 0.867648317476160 0.859246033826186 0.851725509770204 0.844989857576874 0.838950624624483 0.833527638716565 0.828648675160487 0.824249008679563 0.820270897225652 0.816663032341132 0.813379981129823 0.810381637534921 0.807632695011100 0.805102148438438 0.802762829957032];
R1,TOI(t) = [0.550000000000000 1.52334846117846 2.15468610381014 2.54463998749148 2.76412735846199 2.86375701674150 2.87999893220310 2.83930918556651 2.76094243762505 2.65890769477209 2.54335430225640 2.42157130572015 2.29871898421367 2.17837105725229 2.06292049718408 1.95388540065936 1.85214056380102 1.75809317185995 1.67181606841349 1.59314861291571 1.52177266779071 1.45726945893797 1.39916172088524 1.34694453458156 1.30010750107008 1.25815030556572 1.22059326989533 1.18698413529590 1.15690203907239 1.12995943029944 1.10580249849021 1.08411055485884 1.06459470063480 1.04699603471135 1.03108358889593 1.01665212935805 1.00351992451299 0.991526550133554 0.980530780038419 0.970408593747913 0.961051319822417];
ydata = [0.0275000000000000 0.0761674230589228 0.107734305190507 0.127231999374574 0.138206367923100 0.143187850837075 0.143999946610155 0.141965459278326 0.138047121881252 0.132945384738605 0.127167715112820 0.121078565286007 0.114935949210683 0.108918552862615 0.103146024859204 0.0976942700329678 0.0926070281900511 0.0879046585929977 0.0835908034206745 0.0796574306457857 0.0760886333895356 0.0728634729468987 0.0699580860442622 0.0673472267290780 0.0650053750535039 0.0629075152782862 0.0610296634947663 0.0593492067647949 0.0578451019536194 0.0564979715149721 0.0552901249245107 0.0542055277429422 0.0532297350317401 0.0523498017355675 0.0515541794447963 0.0508326064679024 0.0501759962256497 0.0495763275066777 0.0490265390019210 0.0485204296873957 0.0480525659911209];
I am unable to get any good curve fitting results. Any help would be greatly appreciated. I started to modify the code as shown below.
function theta = Gil
% Inital guess for parameters:
Kt0 = 0.25;
Ve0 = 0.4;
theta0 = [Kt0;Ve0];
% lsqcurvefit is in the optimization toolbox.
% fit, from the curve fitting toolbox may be an alternative
theta = lsqcurvefit(@Kdvec,theta0,time,NoisySig_0)
K1 = theta(1)
K2 = theta(2)
% Check fit:
figure
plot(t,ydata,'or')
hold on
xplot = linspace(min(t),max(t));
plot(xplot,Kdvec(theta,xplot))
hold off
end
function Kvec = Kdvec(theta,time)
% Vector of Kd for every xdata:
Kvec = zeros(size(time));
for i = 1:length(time)
Kvec(i) = Kd(theta,time(i));
end
end
function K = Kd(theta,T)
Kt_TOI_fit = theta(1);
Ve_TOI_fit = theta(2);
integrand = @(T) (R1_RR(i)-R10_RR)*exp((-Kt_TOI_fit/Ve_TOI_fit)*(max(time)-T));
Fun1 = @(Kt_TOI_fit,Ve_TOI_fit) ((Kt_TOI_fit/Kt_RR)*(R1_RR(i)-R10_RR));
Fun2= ((Kt_TOI_fit/Kt_RR)*((Kt_RR/Ve_RR)-(Kt_ROI_fit/Ve_TOI_fit)));
K = Fun1 + (Fun2.*integral(integrand,0,max(time))) + R10_TOI;
end
  1 Comment
Are Mjaavatten
Are Mjaavatten on 22 Feb 2021
Edited: Are Mjaavatten on 22 Feb 2021
Programming your model requires a bit of Matlab experience, so it is understandable that you struggled. I took the liberty of modifying your variable names to make them valid Matlab names. Note also that in newer Matlab versions it is possible to combine a script and functions in one file. (Not allowed when I wrote the original answer in Matlab R2014b).
Note that there is nothing special with having an integral in the function to optimize. You use the integral function to calculate the integral and include it with whatever else is needed to evaluate your expression.
t = 0:1:40;
T = 40;
K_RR = 0.1;
v_RR = 0.1;
R10_RR = 0.7;
R10_TOI = 0.55;
R1_RR = [0.700000000000000 0.943337115297993 1.10117152595811 1.19865999687980 1.25353183962319 1.27843925419341 1.28249973305886 1.27232729639958 1.25273560941394 1.22722692370035 1.19833857557102 1.16789282643654 1.13717974605949 1.10709276431873 1.07823012430127 1.05097135016971 1.02553514095478 1.00202329296918 0.980454017107267 0.960787153232550 0.942943166951052 0.926817364737644 0.912290430224259 0.899236133648157 0.887526875270124 0.877037576393890 0.867648317476160 0.859246033826186 0.851725509770204 0.844989857576874 0.838950624624483 0.833527638716565 0.828648675160487 0.824249008679563 0.820270897225652 0.816663032341132 0.813379981129823 0.810381637534921 0.807632695011100 0.805102148438438 0.802762829957032];
R1_TOI = [0.550000000000000 1.52334846117846 2.15468610381014 2.54463998749148 2.76412735846199 2.86375701674150 2.87999893220310 2.83930918556651 2.76094243762505 2.65890769477209 2.54335430225640 2.42157130572015 2.29871898421367 2.17837105725229 2.06292049718408 1.95388540065936 1.85214056380102 1.75809317185995 1.67181606841349 1.59314861291571 1.52177266779071 1.45726945893797 1.39916172088524 1.34694453458156 1.30010750107008 1.25815030556572 1.22059326989533 1.18698413529590 1.15690203907239 1.12995943029944 1.10580249849021 1.08411055485884 1.06459470063480 1.04699603471135 1.03108358889593 1.01665212935805 1.00351992451299 0.991526550133554 0.980530780038419 0.970408593747913 0.961051319822417];
% ydata = R1_TOI*20, and are superfluous
% Anonymous function for interpolating in R1_RR: Could also be calculated inside model.
R1 = @(T) interp1(t,R1_RR,T);
% Inital guess for parameters:
Kt0 = 0.25;
Ve0 = 0.4;
theta0 = [Kt0;Ve0];
% Create anonymous function for the right-hand side of the expression,
% with extra paramters built in:
fvec = @(theta,time) modelvec(theta,time,R1,K_RR,R10_RR,R10_TOI,v_RR);
theta = lsqcurvefit(fvec,theta0,t,R1_TOI);
K_toi = theta(1);
v_toi = theta(2);
% Check fit:
figure
plot(t,R1_TOI,'o')
hold on
tplot = linspace(t(1),t(end));
plot(tplot,fvec(theta,tplot));
hold off
function yvec = modelvec(theta,time,R1,K_RR,R10_RR,R10_TOI,v_RR)
% Vector of y model for a vector of time values:
yvec = zeros(size(time));
for i = 1:length(time)
yvec(i) = model(theta,time(i),R1,K_RR,R10_RR,R10_TOI,v_RR);
end
end
function y = model(theta,T,R1,K_RR,R10_RR,R10_TOI,v_RR)
% Model value, given theta and T
% Fixed parameters, plus function R1, are transferred in the call
K_toi = theta(1);
v_toi = theta(2);
phi = K_toi/v_toi;
R = K_RR/K_toi;
integrand = @(tau) (R1(tau)-R10_RR).*exp(-phi*(T-tau));
I = integral(integrand,0,T);
y = R*((R1(T)-R10_RR) + (K_RR/v_RR-phi)*I)+R10_TOI;
end

Sign in to comment.


Zhao
Zhao on 15 Jun 2021
Hello,
I have a similar problem that share the same background as Chandra and Costas, about the M(H) Curve fitting. I have tried to modify the code that @Are Mjaavatten posted to my problem, but there are erros and can't get any results. I need help to work it out, thanks so much for the suggestions in advance!
The formulation that I try to fit is:
where
and
There are 6 parameters needed, A, Sigma, Ms, Mu, d and C.
The data for fitting are enclosured, the xdata is H and ydata is m(H) in the formulation.
The Codes are below:
xdata = X20210416_Perimag_0_5mg_pro_ml_export(:,1);
ydata = X20210416_Perimag_0_5mg_pro_ml_export(:,2)-X20210521_Water_1_export(:,2);
A_0 = 3e5;
sigma_0 = 0.25;
mu_0 = 4.86;
Ms_0 = 3e5;
C_0 = 2e-9;
para_0 = [A_0:sigma_0;mu_0;Ms_0:C_0];
[para,y] = MH(xdata,ydata,para_0);
plot(xdata,ydata,'ro',xdata,y,'b-');
function [para,y] = MH(xdata,ydata,para_0)
para = lsqcurvefit(@func1,para_0,xdata,ydata);
y = func1(para,xdata);
end
function y = func1(para,xdata)
y = zeros(size(xdata));
for i = 1:length(xdata)
y(i) = func2(para,xdata(i));
end
end
function y = func2(para,H)
intterm = @(d) func3(d,para,H);
C = para(5);
y = integral(intterm,0,200e-9)+C.*H;
end
function ydata = func3(d,H,para)
A = para(1);
sigma = para(2);
mu = para(3);
Ms = para(4);
kb=1.38065e-23;
T=295;
ydata = (A*(pi/6).*d.^3)*(1./(sqrt(2*pi).*d.*sigma)*exp(-log(d./mu).^2/(2*sigma^2)))*(coth((pi/6).*d.^3*Ms*H./kb*T)-(kb*T)./((pi/6).*d.^3*Ms*H));
end
  7 Comments
Are Mjaavatten
Are Mjaavatten on 19 Jun 2021
The reason only mu0 and C0 are modified is obviously the this is all it takes to fit the model to the data. This means that you probably cannot determine more than two parameters from your data. To demonstrate, try another initial value for A0, say 2e17.
There are two possible explanations for this:
Either there are coupling is your model that mean that changing one parameter can always be compensated by changing one of the others. Example: if a and b always appear in the combination a+b, then you cannot determine both from the data. I am not sure that mu and sigma are sufficiently independent, especially since you only use the integral, which says nothing about the shape of the curve.
Or: your experiment / data do not explore all relevant combinations of inputs.
My guess would be a combination.
Zhao
Zhao on 19 Jun 2021
Hello Mr. Mjaavatten,
thank you for your explanation! I think you are right, some parameters may not be independent, I willl have a more detailed research about the formulation and also have a discussion with my superviser about the problem, if there's anything change in the formulation and when I face problems, I will ask you for help again.
Thank you so much for your effort, your time on solving my problem and thank you for your help!

Sign in to comment.


Joor Driez
Joor Driez on 7 Apr 2022
Hello all,
Would someone be able to explain me how to change the code from Mr. Are Mjaavatten if the function 'f(r)' is changing for each molecule size. So basically to alter the to be integrated function ( gauss = @(r) exp(-0.5*((r-rp)/sp).^2); ) in each iteration?
Best regards,
Joor
  3 Comments
Joor Driez
Joor Driez on 8 Apr 2022
The function form is staying constant, but will have more parameters, such that f(r) becomes f(r, a1, a2, a3 ,a4). An examplary function would:
I have added 4 more columns to the data set, containing the parameters a1-a4 for each molecule size. 'gauss' therefore needs to change every i-increment with respect to these parameters, but there I get conflicts.
Are Mjaavatten
Are Mjaavatten on 8 Apr 2022
Edited: Are Mjaavatten on 8 Apr 2022
If a1,...,a4 are functions of the radius rm, I guess you can treat them as extra parameters. Add columns for the a values in the input and pass them to your modified¨Gauss function. lsqcurvefit assumes only one independent variable x (rm in this case), so you must use an anonymous function to 'hide' the parameters.
Example: (will not run, because my modified Gauss does not make sense with my arbitrary a parameters. Hopefully, it wokrs better with realistic function and paramters.
%Driezen.m
data = [
% rm a1 a2 y
0.5 3.7 4 1
1.1 2.4 4 0.83
1.6 0.5 3 0.74
2.2 0.2 4 0.55
2.5 0.4 3.5 0.28
3.5 0.5 0 0
];
xdata = data(:,1);
ydata = data(:,2);
adata = data(:,3:4);
% Inital guess for parameters:
rp0 = 1;
rs0 = 1;
theta0 = [rp0;rs0];
% 'hide' the a parameters for use an anonymous function thuse an anonymous function in the anonymous function fun
fun = @(theta,xdata) Kdvec(theta,xdata,adata);
theta = lsqcurvefit(fun,theta0,xdata,ydata);
function Kvec = Kdvec(theta,xdata,adata)
% Vector of Kd for every xdata:
Kvec = zeros(size(xdata));
for i = 1:length(xdata)
Kvec(i) = Kd(theta,xdata(i),adata(i,:));
end
end
function K = Kd(theta,rm,a)
% Calculates integral. Assumes integrand(1000) is negligible
% Use a lower upper limit if appropriate
rp = theta(1);
sp = theta(2);
gauss = @(r) exp(a(1)*((a(2)*r-rp)/sp).^2);
integrand = @(r) gauss(r).*(1-(rm./r).^2);
K = integral(integrand,rm,1000)/integral(gauss,0,1000);
end

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!