Fitting data to custom integral function

Hi everyone,

I struggle with fitting my measured data to an integral function. I already tried this with the curve fitting toolbox and "lsqnonlin" from the Optimization Toolbox.

Let's say I have the following data set:

x = [228.1194 179.7485 149.5914 121.6736 91.7255 60.6427 40.8913 23.9615 14.4217 11.9658 9.7682 7.4930 5.4940 3.8771 2.6096];

y = [0.7440 0.7349 0.7276 0.7049 0.6939 0.6607 0.6417 0.6069 0.5868 0.5818 0.5781 0.5748 0.5704 0.5606 0.5611];

I want to fit the the following equation to the data set:

y =

where a,b,c,d are the coeffients to fit. i think my problem is the integration part of the fitting function

I tried different approaches.

1) With Curve Fitting Toolbox

%% data set
x = [228.1194  179.7485  149.5914  121.6736   91.7255   60.6427   40.8913   23.9615   14.4217   11.9658    9.7682    7.4930    5.4940    3.8771 2.6096];
y = [0.7440    0.7349    0.7276    0.7049    0.6939    0.6607    0.6417    0.6069    0.5868    0.5818    0.5781    0.5748    0.5704    0.5606 0.5611];
plot(x,y,'x'); 
fun= @(a,b,c,d,x) 1/2/x*(integral(b * (log(1 + a/b * (exp(x/c) - 1))) *exp(-x/d) * (1/(x/c)), 0, 2*x));
g = fittype('(1/2/x*integral(b * (log(1 + a/b * (exp(x/c) - 1))) *exp(-x/d) * (1/(x/c)), 0, 2*x))',...
            'dependent', {'y'},'independent', {'x'}, 'coefficients', {'a','b','c','d'});     
[fitobject,gof] = fit(x, y, g);

Thats the error message:

Error using fittype/testCustomModelEvaluation (line 12)

Expression (integral(b * (log(1 + a/b * (exp(x/c) - 1))) *exp(-x/d) * (1/(x/c)), 0, 2*x)) is not a valid MATLAB expression, has non-scalar

coefficients, or cannot be evaluated:

Error in fittype expression ==> (integral(b .* (log(1 + a./b .* (exp(x./c) - 1))) .*exp(-x./d) .* (1./(x./c)), 0, 2.*x))

??? First input argument must be a function handle.

[...]

2) "lsqnonlin"-approach

xdata = [228.1194  179.7485  149.5914  121.6736   91.7255   60.6427   40.8913   23.9615   14.4217   11.9658    9.7682    7.4930    5.4940    3.8771 2.6096];
ydata = [0.7440    0.7349    0.7276    0.7049    0.6939    0.6607    0.6417    0.6069    0.5868    0.5818    0.5781    0.5748    0.5704    0.5606 0.5611];
plot(xdata,ydata,'x'); 
fun= @(x,xdata) 1/2/xdata*(integral(x(2) * (log(1 + x(1)/x(2) * (exp(xdata/x(3)) - 1))) *exp(-xdata/x(4)) * (1/(xdata/x(3))), 0, 2*xdata)) - ydata;
x0 = [0.65, 0.8, 8, 100000];
x = lsqnonlin(fun, x0);

Error message:

Not enough input arguments.

Error in try2>@(x,xdata)1/2/xdata*(integral(x(2)*(log(1+x(1)/x(2)*(exp(xdata/x(3))-1)))*exp(-xdata/x(4))*(1/(xdata/x(3))),0,2*xdata))-ydata

Error in lsqnonlin (line 196)

initVals.F = feval(funfcn{3},xCurrent,varargin{:});

Error in try2 (line 18)

x = lsqnonlin(fun, x0);

Caused by:

Failure in initial objective function evaluation. LSQNONLIN cannot continue.

Another way I tried first was this, where I solve the problem with a loop and simply guess the parameters:

xdata = [228.1194  179.7485  149.5914  121.6736   91.7255   60.6427   40.8913   23.9615   14.4217   11.9658    9.7682    7.4930    5.4940    3.8771 2.6096];
ydata = [0.7440    0.7349    0.7276    0.7049    0.6939    0.6607    0.6417    0.6069    0.5868    0.5818    0.5781    0.5748    0.5704    0.5606 0.5611];
plot(xdata,ydata,'x'); 
hold on
x = [0.558, 0.78, 25, 100000];
fun = @(xdata) (x(2) .* (log(1 + x(1)./x(2) .* (exp(xdata./x(3)) - 1))) .* exp(-xdata./x(4)) .* (1./(xdata./x(3))));
for i=1:length(xdata)
    y_fit(i) = integral(fun, 0, 2*xdata(i));   
    y_fit2(i) = 1/(2*xdata(i)) *y_fit(i);         
end
plot(xdata,y_fit2,'o');

I found some similar questions in the forum but I wasn't able to adapt these to my case. I'm very thankful for any tip / hint.

Thanks in advance

Answer 1

[Ameer Hamza]

Try the following code. Compare the equations with your code to see the differences. Also, the reason for using arrayfun() is a bit intricate, and I am feeling a bit lazy to write a long explanation :P, so I leave it to you to read the documentation and closely study this code to find out why everything is being done like this. If you have any confusion, then you can ask in the comment.

xdata = [228.1194  179.7485  149.5914  121.6736   91.7255   60.6427   40.8913   23.9615   14.4217   11.9658    9.7682    7.4930    5.4940    3.8771 2.6096];
ydata = [0.7440    0.7349    0.7276    0.7049    0.6939    0.6607    0.6417    0.6069    0.5868    0.5818    0.5781    0.5748    0.5704    0.5606 0.5611];

integral_term = @(a,b,c,d,X) integral(@(x) b.*(log(1 + a./b.*(exp(x./c)-1))).*exp(-x/d).*(1./(x/c)), 0, 2*X);
fun = @(x,xdata) 1/(2*xdata)*integral_term(x(1),x(2),x(3),x(4), xdata);
fun2 = @(x,Xdata) arrayfun(@(xdata) fun(x,xdata), Xdata);

x0 = [0.65, 0.8, 8, 100000];
x = lsqcurvefit(fun2, x0, xdata.', ydata.');

plot(xdata, ydata, '+', xdata, fun2(x, xdata), '-');

[Alex Sha]

I got the result below, seems a little better:

Root of Mean Square Error (RMSE): 0.00278880437155963
Sum of Squared Residual: 0.000116661447342452
Correlation Coef. (R): 0.999135085255474
R-Square: 0.998270918588464

Parameter          	Best Estimate
--------------------	-------------
a	0.554127311813213
b	0.760373675750371
c	29.9709830682305
d	-10085.5078219474

[Ameer Hamza]

It depends on the starting guess and the range of the parameters. In your initial guess, the value of d is positive. If negative values are acceptable, then this solution is also good.

Is this result from MATLAB or some other tool like 1stOpt.

[Alex Sha]

Yes, by using 1stOpt, initial guess is not needed, all were done by softweare itself.

[Van Trung Tin HUYNH]

Thank you for your suggestion. This code works for my problem with Gaussian cumulative distribution function.

Answer 2

[Carlo Vinchi]

Hi Ameer,

thank you very much for your answer! It work's without problems. The arrayfunction was the part I missed.

Like Alex already identified, the fit quality really depends on the starting parameters.

Nevertheless really helpful code. In just 5 lines of code.

Awesome!!

Thanks to both of you ,

Carlo

[Ameer Hamza]

I am glad to be of help.

[Chandra Mouli]

Hi Ameer bro, i am really in chaos in fitting a customized integral based equation.

The main aim is to find the best values of the unknown coefficients in the equation by using the given data and the equation.

In the above equation we have T vs M(T) data set, we should find the best value for Tc and other remaining coefficients by fitting the above equation.

Can u please help me solve this ameer bro ?

thanks,

chandra.