8 views (last 30 days)
Show older comments
Luca Guida on 6 Jul 2021
-
-
Link
Direct link to this question
https://support.mathworks.com/matlabcentral/answers/873343-how-to-fit-integral-function
Answered: Luca Guida on 9 Jul 2021
Accepted Answer: Peter O
Open in MATLAB Online
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
0 Comments Show -2 older commentsHide -2 older comments
Show -2 older commentsHide -2 older comments
Sign in to comment.
Sign in to answer this question.
Accepted Answer
Peter O on 6 Jul 2021
Edited: Peter O on 6 Jul 2021
Open in MATLAB Online
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
0 Comments Show -2 older commentsHide -2 older comments
Show -2 older commentsHide -2 older comments
Sign in to comment.
More Answers (2)
Walter Roberson on 6 Jul 2021
Open in MATLAB Online
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 c*msum() in two directions in order to get the grid.
0 Comments Show -2 older commentsHide -2 older comments
Show -2 older commentsHide -2 older comments
Sign in to comment.
Luca Guida on 9 Jul 2021
Thank you very much to both of you!
0 Comments Show -2 older commentsHide -2 older comments
Show -2 older commentsHide -2 older comments
Sign in to comment.
Sign in to answer this question.
See Also
Categories
Mathematics and OptimizationOptimization ToolboxOptimization ResultsSolver Outputs and Iterative Display
Find more on Solver Outputs and Iterative Display in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- Deutsch
- English
- Français
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)
Contact your local office
