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

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!

Answer 1

[William Rose]

@alfredo scigliani,

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!

[Alfredo Scigliani]

Thank you so much! This definitely helps. I appreciate the help.

Just a follow up, if I wanted to include a term in front of the integral, such as:

to compare which fit of the two equations is better. I tried to just include the term in front like this: (in bold and underlined), but it does not like the modification.

fun2 = @(z,x,D) exp(-x.*D*z.^2); %function inside the integral

%fun2(z,x,D): z=variable of integration, x=data, D=parameter to be estimated

D0 = 1; %initial guess for D

fun = @(D,xdata) exp(-x.*D) .* integral(@(z) fun2(z,xdata,D),0,1,'ArrayValued',true);

%Above: fun(D,xdata)=array-valued integral from z=0 to 1

D = lsqcurvefit(fun, D0, xdata, ydata);

%Above: find D that minimizes squared error between fun(D,xdata) and ydata

fprintf('Best fit: D=%.5f\n',D);

plot(xdata, ydata, '+', xdata, fun(D, xdata), '-');

Any suggestion?

[William Rose]

@alfredo scigliani,

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:

Answer 2

[John D'Errico]

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.