How to fit integral function

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

Answer 1

[Peter O]

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

Answer 2

[Walter Roberson]

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.

Answer 3

[Luca Guida]

Thank you very much to both of you!