Hi,
A couple things I notice:
- 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.
- Because of the way MATLAB's compiler handles anonymous function handles, you'll probably see a significant speedup by creating a couple function files.
- 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
