Curve fitting with an integral involved

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!

Answer 1

[Are Mjaavatten]

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

[Eugenio Gil]

This is great! thank you very much!

[lorrenzo]

Hi,

I have similar problem but the expression has three derived parameter and xdata is upper limit of the integral. I am attaching the expression for reference. I tried to use similar code uploaded by Are Mjaavatten but the code stops saying it did not coverge. Can you please help GIL or Are Mjaavatten.

Code:

function theta = Gil

load 150K.txt;

data=X150K;

xdata = data(:,1);

ydata = data(:,2);

% Inital guess for parameters:

EC = 4.7;

ED = 0.5;

A = 1E5;

theta0 = [EC;ED;A];

% 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,E)

% Calculates integral. Assumes integrand(1000) is negligible

EC = theta(1);

ED = theta(2);

A= theta(3);

gauss = @(E) exp(((E-EC)/ED).^2); %r=E and rp= E0, sp=E1;

integrand = @(E) A.*gauss(E).*((xdata-E).^1.5/(xdata.^3)); %gauss(E).*

K = integral(integrand,0,xdata);

end

[Are Mjaavatten]

Lorrenzo:

In function Kd, the second argument should be h*nu (which is an element in xdata) and not E, which only exists inside the integral. You also forgot a minus sign in gauss. I attach a corrected version that converges, but the fit is quite poor. It seems that the model cannot decribe the data very well.

[lorrenzo]

Thanks Are Mjaavatten. Actually this model is still incomplete. Really appricieate your help.

[Are Mjaavatten]

A comment on the code in my answer:

When I ran this code in Matlab2014b the fitted curve looked OK. I have now obtained a newer version (2020a) and there the fitted curve shows a spurious dip to almost zero for x in the range from x = 1.1 to 1.17. The explanation is:

1) In 2020a lsqcurvefit finds marginally better values for theta.

2) The integrand function sent to integral (in function Kd) is very close to zero everywhere except for a narrow range around 3. It seems that with the new value of theta, the integral function happens to not sample the integrand in this range and thus gives an integral value of around 1e-10.

To fix the problem, change the upper integration limit from 1000 to 10. (Indeed, integrand(x) yields exactly 0 for x >= 7.)

[Chandra Mouli]

I have a 5 parameter model that includes integral function and i started working to fit it by using the above problem, but i am lost somewhere and i am in chaos.

we have the dataset of M(T) vs T , and our main aim is to fit this theoratical fit by using the above formulae and finding the 5 unknown parameters. Especially, i have to find the ''Tc'' for this fit.

My dataset contains of 205 data points so, neither i can paste here nor i can attach the excel sheet. someone who wants to help me and give a solution to this, please revert back to me ASAP. The solution for this is very urgent for me. I hope you guys understand my situation and help me.

Thanks in advance.

Regards,

Chandra Mouli.

[Are Mjaavatten]

In your expression, Tc is the variable of integration and not a parameter to be determined. This does not make sense to me. I think you need to reformulate your problem:

Which are the unknown parameters? (M_0, theta, sigma, mu, Tc?)

What is your variable of integration?

It would help if you also could attach a table of corresponding values of T and M.

[Chandra Mouli]

Yes sir, i think the printing of formulae was wrong the variable of integration is ' T ' not Tc.

The unknown parameters are : M_0 , sigma , mu , Tc and beta.

Theta function is a heaviside function which is if (Tc-T)>0 its ' 1 ' if not (Tc-T)<0 its ' 0 '

T(K) Mag for cooling (M)

749.9995 0.00494

748.0492 0.00646

746.0591 0.00662

744.0491 0.00581

742.0703 0.00829

740.0906 0.00656

738.109 0.00649

736.1371 0.00483

734.1043 0.00527

I hope this data is enough to implement the fit. Thanks in advance sir.

Regards,

Chandra Mouli.

[Are Mjaavatten]

Still a bit unclear: If T is the variable of integration, the right-hand side is not a function of T. Should it be M(Tc) instead of M(T) on the left hand side? And what is beta? Do you mean sigma?

[Chandra Mouli]

sir, the above is the modified formulae and beta is a unknown parameter (marked with an arrow for clarity).

[Are Mjaavatten]

Sorry, didn't see the beta. But the problem remains: the integral is not a function of T. Just as the integral from a to b of x*dx = (b^2-a^2)/2 is not a function of x.

The optimization routine needs intial values for the search for parameters. So it would be useful if you have some thoughts about reasonable initial values for at least some of the parameters,

[Chandra Mouli]

Hi sir, i spoke with my supervisor regarding this actually and he said that we should find the distribution of the Tc values not the single Tc parameter it seems. And i have the initial parameters for all of the unknowns. So, can we find the distribution of Tc from the below formulae sir ?

[Chandra Mouli]

The initial values of unknown paramters are :

sigma = 70 Kelvin

Beta = 1.12

mu = 330 Kelvin

M0 = 0.99

[Are Mjaavatten]

I attach two files. chandra.m fits parameters to the data. candra_demo.m is a script setting up the data and running chandra for two sets of initial parameter values. Note that the parameter set found by lsqcurvefit differ in the two cases and that neither fit is very good. This is not surprising with such a small data set with a lot of noise. With a larger data set this may possibly change.

Note that the effect of the heaviside function is to change the lower integration limit to T. I have used 1200 as the upper integration limit. Depending on the parameters (most notably mu) this may not always be sufficient.

Hope this helps, and good luck!

[Gianrico LAMURA]

Dear Are,

I succesfully used your example to make a fit with an integral function (the superfluid density of a superconductor). I thank you so much.

My problem is that I do not understand how to display the errorbars on the fitted parameters. How could I do?

thank you

[Are Mjaavatten]

Dear Gianrico

To get a single measure of the goodness of your fit, you can use

[theta,resnorm] = lsqcurvefit(___)

resnorm is defined as sum((fun(x,xdata)-ydata).^2). The root mean square error is also often used:

rms = sqrt(resnorm/numel(x))

One way of illustrating the error is of course to make a plot of your model as a function of x, together with xdata and ydata:

theta = lsqcurvefit(fun,theta0,xdata,ydata);
x = linspace(min(xdata),max(xdata));
plot(x,fun(theta,x),xdata,ydata,'or')

An example is given in my comment to Christina MacAskill, futher down on this page..

[Are Mjaavatten]

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]

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.

Answer 2

[Costas Papadopoulos]

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:

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

[Torsten]

What is "xdata" in "MMS" ? You don't define it.

[Costas Papadopoulos]

xdata represents the values of H/(kT) in the experimental data. Even when I define it, it still doesn't work. In addition I am not sure abbout the use of the special characters .* and .^. Obviously the problems are in the function MMS where I define the model, since the remaining code is as suggested by Are Mjaavatten.

[Torsten]

Please show your modified code.

[Costas Papadopoulos]

The only function that I actually edited compared to the provided is the last one. (in other parts I just changed names):

function WLM = MMS(MH,mu,xadata)
  mi = MH(1);
  sigma = MH(2);
  lognormal = @(mu) 1/(mu.*sigma.*(2*pi).^(1/2)).*exp((-ln(mu/mi).^2/(2.*sigma.^2)));  
  integrand = @(mu) lognormal(mu).*(coth(mu*xdata)-(1/(mu*xdata)));  
  WLM = integral(integrand,0,1000);
end

[Torsten]

xadata is not xdata.

[Costas Papadopoulos]

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

Answer 3

[Costas Papadopoulos]

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

[Torsten]

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]

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.

Answer 4

[Christina MacAskill]

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

[Are Mjaavatten]

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

Answer 5

[Zhao]

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

[Are Mjaavatten]

I am sorry, but I am not able to debug your code within reasonable effort,

The trick when creating complex code like this is to divide and conquer: Start with the low-level functions and test that they run without error and return reasonable results. So I recommend to start bottom-up instead of top-down.

A few items:

You can only fit one data column at a time. Use xdata(:,1), or xdata(:). The three columns differ vastly in value, so the latter does not seem so smart...

Your f(d) function looks strange. With the initial parameters it gives Inf for d = 0 and 0 for all other values.

Use the debugger. It is easy to use, and a great help when you try to figure out what went wrong.

[Zhao]

Hello Mr. Mjaavatten,

thank you so much for your detailed advices and your effort on solving my probelm! I also found that the codes that I wrote was a mess, it's my first time to try to fit curve without cft toolbox because of the integral term, so I just imitated your codes but didn't know the reason to write so. Now I have rewritten a new codes, it runs but there's wrong with the result. The new 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))/1000;
     A0 = 1e17;
     mu0 = 130e-9;
     Ms0 = 3e5;
     sigma0 = 0.25;
     C0 = 2e-9;
     para0 = [A0;mu0;Ms0;sigma0;C0];
     para = lsqcurvefit(@func,para0,xdata,ydata);
     y_fit = func(para,xdata);
     figure;
     plot(xdata,ydata,'or');
     hold on;
     plot(xdata,y_fit,'-b');
     
     
     function m = func(para,xdata)
          m = zeros(size(xdata));
          for i = 1:length(xdata)
              m = MH(para,xdata(i));
          end
     end
     
function y = MH(para,x)
  A = para(1);
  mu = para(2);
  Ms = para(3);
  sigma = para(4);
  C = para(5);
  kb=1.38065e-23;
  T=295;
  lognormal =@(d) exp(-(log(d/mu).^2/(2*sigma^2)))./(d*sigma*sqrt((2*pi)));
  langevin =@(d,x) coth(((d.^3).*pi/6*Ms*x)/(kb*T))-(kb*T)./(pi/6*Ms*x.*d.^3);
  integrand = @(d) ((pi/6)*d.^3).*langevin(d,x).*lognormal(d); 
  y = A*integral(integrand,0,inf)+C*x;
end

I had a thought again, d should be a variable, and the limit of integration is (0,inf), H is the xdata and in this code it is replaced as x. The first column of the data is the value of H, second column is m(H), the third is not needed. The unknown parameter are A, Ms, Sigma, Mu and C.

The problem is, the part 'intergrand' isn't integrated at all and it isn't fitted, only the part C*x is fitted, the returned parameter Ms, Sigma and Mu doesn't change from the estimated parameter at beginning, only A and C changed, the result is as below:

I am that I had taken much your effort on my question, but if it's possible, may I ask you to have a look of my new problem of the codes? Thank you so much for your help!

[Zhao]

Hello Mr.Mjaavatten,

I am not sure about the estimated initial parameter, according to the literature the A should be a value of 10^17, Ms is 10^5, Sigma and C should be the same order as the estimated, Mu is around 10^-9, namely nanoscale.

I am sorry that, the value of the data need to be changed as:

xdata = X20210416_Perimag_0_5mg_pro_ml_export(:,1);

ydata = (X20210416_Perimag_0_5mg_pro_ml_export(:,2)-X20210521_Water_1_export(:,2))/1000;

The tested data are initially without units, so after a convertion I think this order should be right. And after running I get a result as below:

Now the parameter A changes but it becomes a complex number and Ms doesn't change, maybe there's still erro by the written codes of 'lognormal' and 'langevin'?

May I also have a ask, why the quality after fitting isn't so good? Are there other function to get a better result? Thank you so much for your help!

[Are Mjaavatten]

I found a couple of problems with your revised code:

It is wise to introduce bounds on the parameter to keep the search from wandering into unphysical values that give nonsense results. In line 9 I therefore use

para = lsqcurvefit(@func,para0,xdata,ydata,zeros(5,1),[1e6;20;1e6;20;20]);

func must return a vector of values, so I modify line 20 to

m(i) = MH(para,xdata(i)); 

Integrate has ploblems with the very wide limits, so I use (line 35):

y = A*integral(integrand,.1,50)+C*x;  

The integrand is negligible outside these limits. Now I get good convergence.

I prefer to use ydata directly instead of the difference. This converges slightly faster and maybe more robustly. The converged parameters are somewhat different, but the results are almost identical.

  xdata = csvread('20210416_Perimag_0.5mg pro ml export.txt');
  zdata = csvread('20210521 Water_1 export.txt');
  xdata = xdata(:,1);
  ydata = zdata(:,1); 

It turns out that ydata is very close to a linear function of xdata:

p = polyfit(xdata,ydata,2);
p(1)
ans =
   1.1185e-11

The deviation is small and seemingly random. Are you sure you are using the correct data set?

[Zhao]

Hello Mr. Mjaavatten,

thank you so much for you help and advices! Actually what I measured is the response of magnetic nanoparticle under magnetic field, so the particle is nanoscale. I had a caculation of the value by hand, A should have a magnitude of 1e17, Mu is median of the partcle and schould be around 130e-9, Sigma should be 1e-2 or 1e-3, Ms should be around 1e5 and C is 1e-9. Ya, maybe I need to have a discussion with my superviser about the current problem.

The curve fitting should be carried out between the first column and the second column, the first column is the magnetic field and the second is the response of the particle, so when only compare the 1st column of the two documents, it should be a line.

May I ask you for a advice, why it's the parameter A, Ms and Sigma don't change after the fitting process? I had get a fitting curve like below:

But even so, A, Ms and Sigma don't varify from the estimated initial parameter, what's the possible reason from the code for this problem?

I am sorry that, the litertaure that I reference has failure, it lost a mv=4*pi*10^(-7) in the formulation, the 'langevin' should be:

mv=4*pi*10^(-7);
langevin =@(d,x) coth(((d.^3).*pi/6*Ms*mv*x)/(kb*T))-(kb*T)./(pi/6*Ms*mv*x.*d.^3);

But still, after adding this part, the A, Ms and Sigma don't change.

Thank you so much for your help!

[Are Mjaavatten]

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]

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!

Answer 6

[Joor Driez]

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

[Are Mjaavatten]

Please explain your problem in a little more detail. How does your function vary with moldecule size? Does it contain a size-dependent parameter or does the function form also vary? During iterations (inside lsqcurvefit) the function cannot change, so presumably you mean "for each molecure size"?

[Joor Driez]

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]

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