Nonlinear Curve fitting with integrals

I encountered a nonlinear fitting problem, and the fitting formula is shown in Equation (1), which includes two infinite integrals (in practice, the integration range can be set from 0.01E-6 to 200E-6).

In these formulas, except for x and y being vectors, all other variables are scalars, and Rmedian and sigma are the parameters to be fitted.

I found a related post and tried to write the code based on it, but it keeps reporting errors. The error message seems to indicate that the vector dimensions are inconsistent, preventing the operation from proceeding. However, these functions are all calculations for individual scalars.

Error using  / 
Matrix dimensions must agree.

Error in dsdmain>@(r)1/(2*r*sigma*sqrt(2*pi))*exp(-(log(2*r)-log(2*Rmean)^2)/(2*sigma^2)) (line 13)
    gauss = @(r) 1/(2*r*sigma*sqrt(2*pi))* exp( -(log(2*r)-log(2*Rmean)^2)/(2*sigma^2) );

My question is: Can I refer to the content of this post to solve my problem? If yes, what does this error message mean? If not, how should I resolve my problem? (Note: The range of Rmedian is 1E-6 to 5E-6)

After modifying the code according to Walter Roberson and Torsten's suggestion, the program no longer throws errors. But no matter how I set the initial values on my end, it always prompts:

    Initial point is a local minimum.
Optimization completed because the size of the gradient at the initial point 
is less than the default value of the optimality tolerance.
<stopping criteria details>
theta = initial point
   1.0e-05 *
    0.2000
    0.1000
Optimization completed: The final point is the initial point.
The first-order optimality measure, 0.000000e+00, is less than
options.OptimalityTolerance = 1.000000e-06.
Optimization Metric                                       Options
relative first-order optimality =   0.00e+00  OptimalityTolerance =   1e-06 (default)

I have checked all the formulas and the units of the variables, and I didn't find any problems.

-------------------------------------- Beblow is my code for issue reproduction ----------------------------------------------------------

function testmain()
clc
    function Kvec = model(param,xdata)
        % Vector of Kd for every xdata:
        Kvec = zeros(size(xdata));
        Rmean = param(1);
        Rstd = param(2);
        
        for i = 1:length(xdata)
            model = @(r) unified(xdata(i),r,Rmean,Rstd,delta,Delta,D,lambda);
            Kvec(i) = integral(model,0.01E-6,200E-6);
        end
    end
    function s = unified(g,R,Rmean,Rstd,delta,Delta,D,lambda)
        %unified Unified fitting model for DSD
        %   exponentional combined
        
        factor = 1./(2*R*Rstd*sqrt(2*pi)) *2 ; % int(P(r)) = 0.5,1/0.5=2
        p1 = -(log(2*R)-log(2*Rmean)).^2/(2*Rstd^2);
        c = -2*2.675E8^2*g.^2/D;
        
        tmp = 0;
        for il = 1:length(lambda)
            a2 = (lambda(il)./R).^2;
            an4 = (R/lambda(il)).^4;
            Psi = 2+exp(-a2*D*(Delta-delta))-2*exp(-a2*D*delta)-2*exp(-a2*D*Delta)+exp(-a2*D*(Delta+delta));
            tmp = tmp+an4./(a2.*R.*R-2).*(2*delta-Psi./(a2*D));
        end
        p2 = c*tmp;
        
        s = factor.*exp(p1+p2);
    end
Delta = 0.075;
delta = 0.002;
D = 0.098E-9;
lambda = [2.0816    5.9404    9.2058   12.4044   15.5792   18.7426   21.8997   25.0528   28.2034   31.3521];

g = [ 0.300616, 0.53884, 0.771392, 1.009616, 1.24784, 1.480392, 1.718616, 1.95684, 2.189392, 2.427616, 2.66584, 2.898392 ];
xdata = g;
ydata = [100, 91.16805426, 80.52955192, 67.97705378, 55.1009735,41.87307917, 30.39638776, 21.13515607, 13.7125649, 8.33083767, 5.146756077, 2.79768609];
ydata = ydata/ydata(1); % normalize

% Inital guess for parameters:
Rmean0 = 2E-6;
Rstd0 = 1E-6;
p0 = [Rmean0;Rstd0];
% lsqcurvefit is in the optimization toolbox.
% fit, from the curve fitting toolbox may be an alternative
theta = lsqcurvefit(@model,p0,xdata,ydata,[0.01E-6;0.1E-6],[10E-6,2E-6])
end

Answer 1

[Torsten]

This code works for your supplied data:

Delta = 0.075;
delta = 0.002;
D = 0.098E-9;
gamma = 2.675E8;
Rmin = 0.01e-6;
Rmax = 200e-6;
npoints = 100000;
lambda = [2.0816    5.9404    9.2058   12.4044   15.5792   18.7426   21.8997   25.0528   28.2034   31.3521];
xdata = [ 0.300616, 0.53884, 0.771392, 1.009616, 1.24784, 1.480392, 1.718616, 1.95684, 2.189392, 2.427616, 2.66584, 2.898392 ];
ydata = [100,91.16805426,80.52955192,67.97705378,55.1009735,41.87307917,30.39638776,21.13515607,13.7125649,8.33083767,5.146756077,2.79768609]/100;
% Optimize
mu0 = (Rmin + Rmax)/2;
sigma0 = 0.2;
p0 = [mu0,sigma0];
format long
p = lsqcurvefit(@(p,xdata)fun_trapz(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda,npoints),p0,xdata,ydata,[Rmin,0],[Rmax,Inf])

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

p = 1×2

0.000002705618161 0.001475482837009

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% Plot results
y = fun_trapz(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda,npoints);
hold on
plot(xdata,ydata)
plot(xdata,y)
hold off
grid on

%
function y = fun_trapz(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda,npoints)
  mu = p(1);
  sigma = p(2);
  y = zeros(size(xdata));
  n1 = ceil(npoints*(p(1)-Rmin)/(Rmax-Rmin));
  n2 = npoints - n1;
  r1 = linspace(Rmin,mu,n1);
  r2 = linspace(mu,Rmax,n2);
  r = [r1,r2(2:end)];
  for i = 1:numel(xdata)
     F = fun(r,xdata(i),Delta,delta,D,gamma,Rmin,Rmax,lambda,mu,sigma);
     y(i) = trapz(r,F);
  end
end
%
function I = fun(r,x,Delta,delta,D,gamma,Rmin,Rmax,lambda,mu,sigma)
  alpha = lambda./r.';
  Psi = 2 + exp(-alpha.^2*D*(Delta-delta))-...
        2*(exp(-alpha.^2*D*Delta)+exp(-alpha.^2*D*delta))+...
        exp(-alpha.^2*D*(Delta+delta));
  lnE = -2*(gamma*x)^2/D*sum(alpha.^(-4)./(alpha.^2.*(r.').^2-2).*(2*delta-Psi./(alpha.^2*D)),2);
  E = exp(lnE);
  P = 1./(2*r.'*sigma*sqrt(2*pi)).*exp(-log(r.'/mu).^2/(2*sigma^2));
  I = (2*E.*P).';
end

[Shannon]

In your code, you used trapz to approximate the integral, but the textbook says Gaussian quadrature is more accurate. I modified the code you provided by replacing trapz with Gauss-Legendre integration, but the results turned out worse. Why is this happening?

http://www.deltaquants.com/gaussian-quadrature-gauss-legendre-Integration

weights and abscissa were adopted from https://pomax.github.io/bezierinfo/legendre-gauss.html (n=24)

[Sam Chak]

Hi @Shannon, Have you heard about the Gauss-Kronrod quadrature, quadgk()?

[Torsten]

I modified the code you provided by replacing trapz with Gauss-Legendre integration, but the results turned out worse. Why is this happening?

Available codes try to use a small number of grid points together with a high-order method to get a good approximation of an integral. This assumes that the grid points are placed at the correct positions in the integration interval. But this so called "adaptive gridding" is difficult in case of your function: for small values of sigma, the interesting part of your function is concentrated in a very narrow region. Without further help for the integrator (e.g. by prescribing a very fine grid (as I did in trapz) or by defining Rmean as a Waypoint (as I tested with "integral" and which worked, too), it is very difficult for the integrator to find this region and will fail.

So summarizing: in your case, the integration method doesn't matter, but it's essential that enough grid points are placed at the relevant positions. This failed for MATLAB's "integral" without further support, and it will most probably fail with every other integrator as well without giving it a helping hand (surprisingly, "quadgk" worked even without the waypoints option).

Here is the "Waypoints" code for the use of "integral":

Delta = 0.075;
delta = 0.002;
D = 0.098E-9;
gamma = 2.675E8;
Rmin = 0.01e-6;
Rmax = 200e-6;
npoints = 100000;
lambda = [2.0816    5.9404    9.2058   12.4044   15.5792   18.7426   21.8997   25.0528   28.2034   31.3521];
xdata = [ 0.300616, 0.53884, 0.771392, 1.009616, 1.24784, 1.480392, 1.718616, 1.95684, 2.189392, 2.427616, 2.66584, 2.898392 ];
ydata = [100,91.16805426,80.52955192,67.97705378,55.1009735,41.87307917,30.39638776,21.13515607,13.7125649,8.33083767,5.146756077,2.79768609]/100;
% Optimize
mu0 = (Rmin + Rmax)/2;
sigma0 = 0.2;
p0 = [mu0,sigma0];
format long
p = lsqcurvefit(@(p,xdata)fun_integral(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda),p0,xdata,ydata,[Rmin,0],[Rmax,Inf])

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 4.1e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 8.7e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 4.7e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.8e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 5.4e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.2e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 2.1e-09. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 7.6e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 8.3e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 4.6e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.9e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 5.5e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.3e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 2.1e-09. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.2e-10. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 4.1e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 8.7e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 4.7e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.8e-07. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 5.4e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 1.2e-08. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 2.1e-09. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

p = 1×2

0.000002705623324 0.002304341838077

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% Plot results
y = fun_integral(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda);
hold on
plot(xdata,ydata)
plot(xdata,y)
hold off
grid on

%
function y = fun_integral(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda)
  mu = p(1);
  sigma = p(2);
  for i = 1:numel(xdata)
     F = @(r)fun(r,xdata(i),Delta,delta,D,gamma,Rmin,Rmax,lambda,mu,sigma);
     y(i) = integral(F,Rmin,Rmax,'Waypoints',mu);
  end
end
%
function I = fun(r,x,Delta,delta,D,gamma,Rmin,Rmax,lambda,mu,sigma)
  alpha = lambda./r.';
  Psi = 2 + exp(-alpha.^2*D*(Delta-delta))-...
        2*(exp(-alpha.^2*D*Delta)+exp(-alpha.^2*D*delta))+...
        exp(-alpha.^2*D*(Delta+delta));
  lnE = -2*(gamma*x)^2/D*sum(alpha.^(-4)./(alpha.^2.*(r.').^2-2).*(2*delta-Psi./(alpha.^2*D)),2);
  E = exp(lnE);
  P = 1./(2*r.'*sigma*sqrt(2*pi)).*exp(-log(r.'/mu).^2/(2*sigma^2));
  I = (2*E.*P).';
end

[Shannon]

The code I'm using adopted from P313 of https://porousmedia.rice.edu/resources/PhD%20Thesis.pdf , which implements 24-point Gauss-Legendre integration. But my data isn't producing correct results.

[Shannon]

@Sam Chak How amazing, the Gauss–Kronrod quadrature actually works!

Answer 2

[Walter Roberson]

integral() always passes in a vector of locations to calculate at, unless you use the 'ArrayValued', true option.

    sig = integral(f,0.01E-6,200E-6)/integral(gauss,0.01E-6,200E-6);

refers to

    f = @(r) gauss(r)*kernel(r,delta,Delta,D)*2;

with the input r being a vector, gauss( r) can be expected to be a vector and kernal(...) can be expected to be a vector, so you have set up for two vectors to be "*" together. As the two vectors are likely the same size, you can expect problems from the * operator which does algebraic matrix multiplication for which the "inner dimensions" must match. You need to use the .* operator.

Meanwhile,

    gauss = @(r) 1/(2*r*sigma*sqrt(2*pi))* exp( -(log(2*r)-log(2*Rmean)^2)/(2*sigma^2) );

with the input r being a vector, the 1/(2*r*sigma*sqrt(2*pi)) component is 1/(vector) which is going to fail because the / operator is the matrix-right-divide operator in which the number of columns on the two sides must match. You need to use the ./ operator.

Then the exp() is going to be a vector and you have a * operator again, and that's going to fail because the inner dimensions must be the same size... again you need to be using the .* operator instead of the * operator.

[Shannon]

If that's the case, Equation (3) cannot be implemented in MATLAB because both g and R are vectors with different dimensions.

[Torsten]

Equation (3) (and thus Equation (1)) must be evaluated for each x-value separately. After integration, Equation (1) will give a vector with the same length as x (or y).

[Shannon]

Can you get the fitting results? No matter how I set the initial values on my end, it always prompts :

Initial point is a local minimum.

Optimization completed because the size of the gradient at the initial point 
is less than the default value of the optimality tolerance.

<stopping criteria details>


theta =

   1.0e-05 *

    0.2000
    0.1000


Optimization completed: The final point is the initial point.
The first-order optimality measure, 0.000000e+00, is less than
options.OptimalityTolerance = 1.000000e-06.

Optimization Metric                                       Options
relative first-order optimality =   0.00e+00  OptimalityTolerance =   1e-06 (default)

[Torsten]

p1 is in the order of -1e13. This makes s = 0 for all values of r and thus Kvec identically 0. So lsqcurvefit does not notice any change in the y-values you return when it varies the parameters to be optimized and consequently stops the optimization. Thus my advice: Check your model and your model parameters for validity.

[Shannon]

I checked all the formulas and all the variables within the formulas, including the units of the variables, but I did not find any issues.

[Shannon]

The units of all variables are in the International System of Units (SI), so the unit of the parameter RR to be fitted is meters (m). Since the optimal solution is near 2 micrometers (2 µm), this problem inherently suffers from an imbalance in magnitude that cannot be resolved through scaling. How should such a situation be addressed?

[Torsten]

But you see that giving sigma an initial value of 1e-6 gives a value for P(R) = 0 because you divide by 2*1e-12 in the exp() expression ?

[Shannon]

It appears that the integral function discretizes R into 150 points between the integration limits. Due to the very narrow distribution, if these discretized points miss Rmean, the computed p1p1 becomes entirely zero.

[Walter Roberson]

That's what the 'waypoints' option is for.

[Shannon]

I don't know why setting waypoints to Rmean doesn't work. Could the presence of ln() cause the waypoints to fail?

[Torsten]

It doesn't help if your function is different from 0 in only one point - the numerical value of the integral will remain 0. If you really expect "sigma" in the order of magnitude of 1e-6, then it is impossible to use a numerical integrator for your task. If you try to plot P(r), you will see that you almost get the dirac distribution that is Inf at r = Rmean and 0 else.

Go one step back and try plotting P(r)*E(r,x,delta,s,D) over r for some of your x-values and see up to which value of "sigma" you get a reasonable function that you think MATLAB could successfully integrate. Then see if this value is still acceptable for your fitting task to start with.

Just a sidenote: Are you sure about the 2 in the denominator of the prefactor for P(R) ? The integral evaluates to 0.5 whereas I think it should evaluate to 1 as the probability density function for a lognormal distribution.

On Mathematics Stack Exchange you write

αm is the positive root of the equation J3/2(αmR)=αmRJ5/2(αmR)

So you mean that for each value of R, you get a different vector αm with the first 10 positive roots of the above equation ?

Maybe you could give a link to a publication where the above model has been used successfully to fit the parameters sigma and Rmedian.

[Shannon]

https://www.sciencedirect.com/science/article/pii/S0167732224008493

This processing procedure is widely adopted in the literature. The only distinction lies in the slight variation of the constant term preceding the summation symbol in Eq. (3), which depends on the sampling equipment.

[Shannon]

So you mean that for each value of R, you get a different vector αm with the first 10 positive roots of the above equation ?

Yes, that's why I substituted ‌am‌ with ‌R‌ via ‌am = lambda / R

[Torsten]

Yes, that's why I substituted am with R via am = lambda / R

Ah, you solve for the first 10 positive roots for lambda of the equation J_3/2(lambda) = lambda*J_5/2(lambda) .

Is y(i)/y(1) = P(D>=x(i)) where D is the droplet size in mu-m ?

[Shannon]

R is the droplet size (radius) and the distribution P® is supposed to be a lognornal distribution. Eq(2)

Answer 3

[Torsten]

It takes too long to run the code online - try if it produces reasonable results. Maybe you have to set lower and upper bounds for the parameters in "lsqcurvefit".

Note that "mu" and "sigma" are not the expected value and the standard deviation of the Lognormal Distribution. Further I think that "sigma" is dimensionless.

Delta = 0.075;
delta = 0.002;
D = 0.098E-9;
gamma = 2.675E8;
Rmin = 0.01e-6;
Rmax = 200e-6;
lambda = [2.0816    5.9404    9.2058   12.4044   15.5792   18.7426   21.8997   25.0528   28.2034   31.3521];
xdata = [ 0.300616, 0.53884, 0.771392, 1.009616, 1.24784, 1.480392, 1.718616, 1.95684, 2.189392, 2.427616, 2.66584, 2.898392 ];
ydata = [100,91.16805426,80.52955192,67.97705378,55.1009735,41.87307917,30.39638776,21.13515607,13.7125649,8.33083767,5.146756077,2.79768609]/100;

mu0 = (Rmin + Rmax)/2;
sigma0 = 0.2;

% Plot integrand for initial values for mu and sigma
N = 100;
r = linspace(Rmin,Rmax,N);
I = arrayfun(@(x)fun_int(r,x,Delta,delta,D,gamma,Rmin,Rmax,lambda,mu0,sigma0),xdata,'UniformOutput',0);
plot(r,cell2mat(I))
grid on

% Optimize mu and sigma
p0 = [mu0,sigma0];
p = lsqcurvefit(@(p,xdata)fun(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda),p0,xdata,ydata)

function y = fun(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda)
  mu = p(1);
  sigma = p(2);
  y = zeros(numel(xdata),1);
  for i = 1:numel(xdata)
    F = @(r) fun_int(r,xdata(i),Delta,delta,D,gamma,Rmin,Rmax,lambda,mu,sigma);
    y(i) = integral(F,Rmin,Rmax);
  end
end

function I = fun_int(r,x,Delta,delta,D,gamma,Rmin,Rmax,lambda,mu,sigma)
  I = zeros(numel(r),1);
  for i = 1:numel(r)
    alpha = lambda/r(i);
    Psi = 2 + exp(-alpha.^2*D*(Delta-delta))-...
          2*(exp(-alpha.^2*D*Delta)+exp(-alpha.^2*D*delta))+...
          exp(-alpha.^2*D*(Delta+delta));
    lnE = -2*(gamma*x)^2/D*sum(alpha.^(-4)./(alpha.^2*r(i)^2).*(2*delta-Psi./(alpha.^2*D)));
    E = exp(lnE);
    P = 1./(2*r(i)*sigma*sqrt(2*pi))*exp(-log(r(i)/mu).^2/(2*sigma^2));
    I(i) = 2*E*P;
  end
end

[Shannon]

Error using integralCalc/finalInputChecks (line 526)
Output of the function must be the same size as the input. If FUN is an array-valued integrand, set the 'ArrayValued' option to true.

Error in integralCalc/iterateScalarValued (line 315)
                finalInputChecks(x,fx);

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 88)
Q = integralCalc(fun,a,b,opstruct);

Error in Untitled2>fun (line 30)
    y(i) = integral(F,Rmin,Rmax);

Error in Untitled2>@(p,xdata)fun(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda)

Error in lsqcurvefit (line 213)
            initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});

Error in Untitled2 (line 22)
p = lsqcurvefit(@(p,xdata)fun(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda),p0,xdata,ydata)

Caused by:
    Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.

[Walter Roberson]

Delta = 0.075;
delta = 0.002;
D = 0.098E-9;
gamma = 2.675E8;
Rmin = 0.01e-6;
Rmax = 200e-6;
lambda = [2.0816    5.9404    9.2058   12.4044   15.5792   18.7426   21.8997   25.0528   28.2034   31.3521];
xdata = [ 0.300616, 0.53884, 0.771392, 1.009616, 1.24784, 1.480392, 1.718616, 1.95684, 2.189392, 2.427616, 2.66584, 2.898392 ];
ydata = [100,91.16805426,80.52955192,67.97705378,55.1009735,41.87307917,30.39638776,21.13515607,13.7125649,8.33083767,5.146756077,2.79768609]/100;

mu0 = (Rmin + Rmax)/2;
sigma0 = 0.2;

% Plot integrand for initial values for mu and sigma
N = 100;
r = linspace(Rmin,Rmax,N);
I = arrayfun(@(x)fun_int(r,x,Delta,delta,D,gamma,Rmin,Rmax,lambda,mu0,sigma0),xdata,'UniformOutput',0);
whos r I
size(I{1})
plot(r,cell2mat(I(:)))
grid on
% Optimize mu and sigma
p0 = [mu0,sigma0];
p = lsqcurvefit(@(p,xdata)fun(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda),p0,xdata,ydata)

function y = fun(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda)
  mu = p(1);
  sigma = p(2);
  y = zeros('like',xdata);
  for i = 1:numel(xdata)
    F = @(r) fun_int(r,xdata(i),Delta,delta,D,gamma,Rmin,Rmax,lambda,mu,sigma);
    y(i) = integral(F,Rmin,Rmax);
  end
end

function I = fun_int(r,x,Delta,delta,D,gamma,Rmin,Rmax,lambda,mu,sigma)
  I = zeros('like', r);
  for i = 1:numel(r)
    alpha = lambda/r(i);
    Psi = 2 + exp(-alpha.^2*D*(Delta-delta))-...
          2*(exp(-alpha.^2*D*Delta)+exp(-alpha.^2*D*delta))+...
          exp(-alpha.^2*D*(Delta+delta));
    lnE = -2*(gamma*x)^2/D*sum(alpha.^(-4)./(alpha.^2*r(i)^2).*(2*delta-Psi./(alpha.^2*D)));
    E = exp(lnE);
    P = 1./(2*r(i)*sigma*sqrt(2*pi))*exp(-log(r(i)/mu).^2/(2*sigma^2));
    I(i) = 2*E*P;
  end
end

However, running this generates a long list of messages similar to

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is 9.3e+45. The integral may not exist, or it may be difficult to approximate numerically to the requested accuracy.

> In integralCalc>iterateScalarValued (line 393)

In

integralCalc>vadapt (line 148)

In

integralCalc (line 77)

In

integral (line 87)

In

q2175948>fun (line 32)

In

q2175948>@(p,xdata)fun(p,xdata,Delta,delta,D,gamma,Rmin,Rmax,lambda) (line 24)

In

lsqcurvefit/objective (line 314)

In finDiffEvalAndChkErr

In finitedifferences

In computeFinDiffGradAndJac

In

sfdnls (line 51)

In

snls (line 345)

In

lsqncommon (line 178)

In

lsqcurvefit (line 299)

In

q2175948 (line 24)

Notice how large the error bound is. There is no hope of fixing error bounds that large with small tweaks.

[Torsten]

@Shannon Walter Roberson corrected the problem with your error message (y and I were column instead of row vectors), but the integration seems to be unsuccessful with the new suggested parameters for mu and sigma by "lsqcurvefit" after the first optimization step. If mu, e.g., becomes negative, do as I suggested and impose bounds in the parameters. I don't know if it helps. If you look at the plot of the integrand over r for xdata(end), e.g., it just looks like a peak - thus almost impossible to integrate with numerical methods.

I forgot to subtract 2 here, but I think this won't change much:

lnE = -2*(gamma*x)^2/D*sum(alpha.^(-4)./(alpha.^2*r(i)^2-2).*(2*delta-Psi./(alpha.^2*D)));

instead of

lnE = -2*(gamma*x)^2/D*sum(alpha.^(-4)./(alpha.^2*r(i)^2).*(2*delta-Psi./(alpha.^2*D)));

[Shannon]

Origin can perform the fitting, but MATLAB fails? ‌

From one of the authors:

We did fit the dataset you provided to test in the way presented in the paper. Assuming a monomodal log-normal distribution (volume based) and using 10 roots (and not including free component or other effects), we get RMedian = 2.706 μm +/- 0.031 μm and log normal width factor of 0.0024 +/- 2.8231

• The integration was included in the Origin Fitting function and defined using a LabTalk Script but see the help file in Origin help for integration notes. If you go to Origin help – All Books – Tutorials – Data Analysis (i.e., 4 Data Analysis) and look in the section on nonlinear fitting there are some sheets there that should help you create your function in Origin. Look at Fitting with Integral using LabTalk Function and Fitting with Two Integrals using LabTalk Function. There are some other files or ways in there, that may also help.
• The help indicates that Origin also allows for +inf/-inf for limits, however, as mentioned we limited our integral as discussed after consideration of the probed range. Depending on the limits you use and the integral functions you may also need to split the range as indicated in the help file.
• The origin help files we have indicate that using integrals needs version 8.6SR0 (or 8.0SR6 depending on the help file page) and up.
• So, we defined the two integrals and then called on them in the function and instead of using +inf we use an upper limit as mentioned. The integrals themselves are done using Origin capability (numerically) (see also the help file for the math functions in LabTalk for integral and other files there). The function summation (with roots) was written in expanded form that way the roots were entered (fixed parameters) up to the number used but there is likely other ways as well.

[Torsten]

So you sent your data to the authors of the article and got Rmean and sigma for your data ? And how good is the fit of the y-vector when you use these parameters ? The uncertainty in the width-factor is strange, isn't it ? Especially it can become negative ...

[Shannon]

I send my data in this post to the authors and they give me the fitting result back. The uncertainty in the width-factor is strange, but at least it can complete the fitting, and the fitted y matches the original y-vector very well.

[Shannon]

xdata ydata [%] Fitting [%]

0.300616 100 100

0.53884 91.16805426 91.83196659

0.771392 80.52955192 80.6936292

1.009616 67.97705378 67.4689156

1.24784 55.1009735 53.86768785

1.480392 41.87307917 41.40418259

1.718616 30.39638776 30.28783311

1.95684 21.13515607 21.24404364

2.189392 13.7125649 14.45360593

2.427616 8.33083767 9.377799625

2.66584 5.146756077 5.865511151

2.898392 2.79768609 3.587601962