필터 지우기
필터 지우기

Curve fitting with an integral involved

조회 수: 24 (최근 30일)
Eugenio Gil
Eugenio Gil 2018년 7월 17일
편집: Are Mjaavatten 2022년 4월 8일
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!

채택된 답변

Are Mjaavatten
Are Mjaavatten 2018년 7월 17일
편집: Are Mjaavatten 2018년 7월 17일
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
Are Mjaavatten
Are Mjaavatten 2021년 6월 28일
편집: Are Mjaavatten 2021년 6월 28일
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 2021년 6월 30일
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.

댓글을 달려면 로그인하십시오.

추가 답변 (5개)

Costas Papadopoulos
Costas Papadopoulos 2019년 9월 2일
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
Torsten
Torsten 2019년 9월 2일
xadata is not xdata.
Costas Papadopoulos
Costas Papadopoulos 2019년 9월 2일
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

댓글을 달려면 로그인하십시오.


Costas Papadopoulos
Costas Papadopoulos 2019년 9월 5일
편집: Costas Papadopoulos 2019년 9월 5일
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
Torsten
Torsten 2019년 9월 5일
편집: Torsten 2019년 9월 5일
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 2019년 9월 5일
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.

댓글을 달려면 로그인하십시오.


Christina MacAskill
Christina MacAskill 2021년 2월 15일
편집: Christina MacAskill 2021년 2월 15일
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
Are Mjaavatten
Are Mjaavatten 2021년 2월 22일
편집: Are Mjaavatten 2021년 2월 22일
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

댓글을 달려면 로그인하십시오.


Zhao
Zhao 2021년 6월 15일
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
Are Mjaavatten
Are Mjaavatten 2021년 6월 19일
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 2021년 6월 19일
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!

댓글을 달려면 로그인하십시오.


Joor Driez
Joor Driez 2022년 4월 7일
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
Joor Driez
Joor Driez 2022년 4월 8일
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 2022년 4월 8일
편집: Are Mjaavatten 2022년 4월 8일
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

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 MATLAB에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by