How to fit a function with an integral term to a set of data?

2 views (last 30 days)
I have this set of data of 32 points (xdata and ydata) and a function that I am trying to fit. This function includes an integral term with one parameter (D), which I think is causing the difficulty.
This is the function that I am trying to fit to the data.
clear; clc;
xdata=[10.30 29.88 59.64 99.58 149.66 209.96 280.44 361.03 451.87 552.89 664.10 785.38 916.94 1058.68 1210.48 1372.58 1544.86 1727.33 1919.81 2122.64 2335.65 2558.64 2792.01 3035.55 3289.29 3552.97 3827.05 4111.33 4405.52 4710.15 5024.96 5349.96];
ydata=[1.00 0.96 0.93 0.87 0.82 0.78 0.72 0.67 0.61 0.56 0.49 0.47 0.43 0.39 0.37 0.34 0.31 0.30 0.29 0.27 0.26 0.26 0.25 0.25 0.24 0.24 0.25 0.23 0.24 0.24 0.21 0.21];
fun = integral(@(x) exp(-xdata.*x(1).*x.^2),0,1);
x0 = 1; %initial guess
x = lsqcurvefit(fun, x0, xdata, ydata);
plot(xdata, ydata, '+', xdata, fun2(x, xdata), '-');
This code runs unsuccessfully and this is the message I get:
Matrix dimensions must agree.
Error in Untitled>@(x)exp(-xdata.*x(1).*x.^2) (line 5)
fun = integral(@(x) exp(-xdata.*x(1).*x.^2),0,1);
Error in integralCalc/iterateScalarValued (line 314)
fx = FUN(t);
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 Untitled (line 5)
fun = integral(@(x) exp(-xdata.*x(1).*x.^2),0,1);
final comments: I am not sure if the problem is trying to use lsqcurvefit with a integral function with one parameter, or if I should be using a different method. I am open to any suggestion. Thank you so much in advance!

Accepted Answer

William Rose
William Rose on 2 Nov 2021
This fitting problem is complicated, since the function to be fitted is an integral.
You can pass the parameter D and the xdata to the integral by defining the function that is inside the integral. You must specify that the integral is "ArrayValued" so that it returns a vector with the same size as ydata.
There are different meanings of x in this problem: x=parameter to fit, x=variable of integration, and xdata=x values of the data. In order to avoid confusion, I have defined the parameter to be fitted as D, and the variable of integration is z. See attached code.
The code produces some warning messages. Despite the warnings, the code sucessfully finds a best-fit solution for D.
>> modelFitWithIntegral
...
Best fit: D=0.00438
The graph shows that the fit is reasonable.
Good luck!
  2 Comments
William Rose
William Rose on 3 Nov 2021
Do
fun = @(D,xdata) exp(-xdata*D).*integral(@(z) fun2(z,xdata,D),0,1,'ArrayValued',true);
This works. D=0.00052. The fit is not as good.
I agree with @John D'Errico that it is more elegant to use the erf function instead of the integral. I compliment him for noticing that. You can prove that Matlab's symbolic math result is correct by substitution of variables:

Sign in to comment.

More Answers (1)

John D'Errico
John D'Errico on 2 Nov 2021
Edited: John D'Errico on 2 Nov 2021
Certainly for this problem, it is easiest to just perform the integration in advance. The integral you show will have a solution in terms of the erf function (often known as a special function in mathematics.)
Thus, we can do:
syms Xd D x
I = int(exp(-Xd*D*x.^2),[0,1])
I = 
As I said, we see an erf in there. Next, can we fit this form, given a vector xdata and ydata, can we find the optimal value for D? Of course. The curve fitting toolbox is best, because it makes the problem simple. And then we get nice plots, etc.
Ifun = @(D,X) sqrt(pi)*erf(sqrt(D*X))./(2*sqrt(D*X))
Ifun = function_handle with value:
@(D,X)sqrt(pi)*erf(sqrt(D*X))./(2*sqrt(D*X))
mdl = fittype(Ifun,'indep','X')
mdl =
General model: mdl(D,X) = sqrt(pi)*erf(sqrt(D*X))./(2*sqrt(D*X))
xdata=[10.30 29.88 59.64 99.58 149.66 209.96 280.44 361.03 451.87 552.89 664.10 785.38 916.94 1058.68 1210.48 1372.58 1544.86 1727.33 1919.81 2122.64 2335.65 2558.64 2792.01 3035.55 3289.29 3552.97 3827.05 4111.33 4405.52 4710.15 5024.96 5349.96];
ydata=[1.00 0.96 0.93 0.87 0.82 0.78 0.72 0.67 0.61 0.56 0.49 0.47 0.43 0.39 0.37 0.34 0.31 0.30 0.29 0.27 0.26 0.26 0.25 0.25 0.24 0.24 0.25 0.23 0.24 0.24 0.21 0.21];
estmdl = fit(xdata',ydata',mdl,'start',.001)
estmdl =
General model: estmdl(X) = sqrt(pi)*erf(sqrt(D*X))./(2*sqrt(D*X)) Coefficients (with 95% confidence bounds): D = 0.004376 (0.004178, 0.004574)
plot(estmdl)
hold on
plot(xdata,ydata,'bo')
xlabel X
ylabel Y
grid on
The fit seems reasonable.
Could I have done it by performing an integration in the call itself, thus, using integral? Well, yes. But that would be far less efficient. As well, we would now need to use tools like arrayfun. Far better to recognize the integral is integrable as a special function.

Products


Release

R2019b

Community Treasure Hunt

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

Start Hunting!