How to fit integral function

5 views (last 30 days)
Luca Guida
Luca Guida on 6 Jul 2021
Answered: Luca Guida on 9 Jul 2021
Hi everyone, I have a function in two variables containing 5 parameters. I need to fit it to experimental data in order to find the values of parameters. The function is the following:
I have column vectors containing corresponding values of x, y and fun from experimental data.
The problem arises when i try to write the integral function:
%% fictitious experimental data
x=linspace(0,30,100);
y=linspace(0,50,100);
z=linspace(0,40,100);
% a=t(1)
% b=t(2)
% c=t(3)
% d=(4)
% e=t(5)
C=@(t) ((x - t(3))./t(1)).^2-2.*t(5).*((x - t(3))./t(1)).*((y - t(2))./t(4))+((y - t(2))./t(4)).^2;
k=@(t) 1./(t(1).*t(2)*sqrt(t(5))).*exp(-C(t)./t(5));
fun=@(t) integral2(k, 0, x, 0, y)-z;
t0=[1 1 1 1 1];
coeff = lsqnonlin(fun,t0);
What i get is
Error using integral2 (line 80)
XMAX must be a floating point scalar.
Does anyone know i can i express the integral in a better way and how to perform to the fitting phase?
Thank you in advance

Accepted Answer

Peter O
Peter O on 6 Jul 2021
Edited: Peter O on 6 Jul 2021
Hi,
A couple things I notice:
  1. integral2 expects max and min bounds, but you're passing it the entire grid. Instead of x and y, it needs 30 and 50 in the dummy data above.
  2. Because of the way MATLAB's compiler handles anonymous function handles, you'll probably see a significant speedup by creating a couple function files.
  3. For a minimization problem, consider targeting the sum of the squared error. I'm not sure of your actual data sensitivity, but also bear in mind that you want to try to keep the parameter weights close. With lsqnonlin it could be less of a concern depending on the solution topology, but it's not a very pretty surface. If there's broad variation in the magnitude of Z the algorithm might get the biggest parts right at the expense of the smaller magnitude regions, and you'd see that reflected in the parameter sensitivity. That is F(x,y) looks much better in some regions then other. With five parameters and an exponent implicit inside the integral, there's a strong likelihood that your initial guess has to be good.
Try setting it up this way. This is pseudocode, and unfortunately I'm without the optimization toolbox at the moment, so there might be a little error debugging to get it to work.
Edit: lsqnonlin works on a vector error function. My bad! Updated the target function output.
%% MyOptimizerScript.m
% fictitious experimental data: 100 sample points
X = linspace(0,30,100);
Y = linspace(0,50,100);
Z = linspace(0,40,100);
t0=[1 1 1 1 1];
% Define a match function with all the inputs that lsqnonlin can call. It
% will try driving this function to zero (minimize error).
matchFunction = @(t) myIntegralFit(t, x, y, z)
coeff = lsqnonlin(@matchFunction,t0);
%% myIntegralFit.m
function target = myIntegralFit(t,X,Y,Z)
% X and Y are measurement point inputs.
% Z is the measurement point output.
% We desire a set of constants, t, which minimize the difference between
% INT(F(x,y) and Z(x,y), for all input collections
% Given this set of parameters t, compute the guess for the output for each
% set of X and Y. Unfortunately, this does seem to require a lot of calls
% to integral2 (on first inspection), since the internal constants vary in
% a fairly ugly way that makes it tricky to vectorize.
IntGuess = zeros(length(x),1)
for ix=1:length(x)
% Loop through x/y combos
IntGuess(ix) = myIntegral(t, x(ix), y(ix))
end
% Compute the error
target = IntGuess - Z;
end
%% myIntegral.m
function z = myIntegral(t,x,y)
% Compute the integral for each "datapoint" set of an x and y bound
myFunc = @(x,y) kfunc(x,y,t) % Locks t to kfunc for this iteration
z = integral2(myFunc, 0, x, 0, y) % Returns a single Z for this x-y combo
end
function k = kfunc(x,y,t)
%KFUNC is going to get a grid of XY points from integral2. Use these to
%batch compute Cxy in a single call.
Cxy = ((x - t(3))./t(1)).^2-2.*t(5).*((x - t(3))./t(1)).*((y - t(2))./t(4))+((y - t(2))./t(4)).^2;
k = 1./(t(1).*t(2)*sqrt(t(5))).*exp(-Cxy./t(5));
end

More Answers (2)

Walter Roberson
Walter Roberson on 6 Jul 2021
fun=@(t) integral2(k, 0, x, 0, y)-z;
Where are you using the input value, t ?
No, the t parameter will not be passed into k just because k happens to be a function whose dummy parameter name is the same as t.
Remember, this is integral2(), so the first parameter (k here) must evaluate to the handle of a function that accepts x and y values. For example,
fun = @(t) integral2(@(x,y)k(x,y,t), 0, x, 0, y) - z;
Also
x=linspace(0,30,100);
Your x is a vector, but you are trying to use that vector as the upper bound for integral2(). integral2() can only accept one scalar bound at a time.
If you want to fit at several x and y then you can break the overall x range into pieces, and the y range into pieces, and integrate over the rectangles, and cumsum() in two directions in order to get the grid.

Luca Guida
Luca Guida on 9 Jul 2021
Thank you very much to both of you!

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!