# nlinfit with modelfun as an integral2 and variable integral limits

2 views (last 30 days)
Peter Seibold on 5 Jul 2020
Commented: Peter Seibold on 7 Jul 2020
I have several volumes over specific areas and want to know which 3D bell shape fits best for those partial volumes. I try to estimate the peak value and sigma of the 3D bell shape by a double integral (‘integral2’) within ‘nlinfit’. I have seen a solution for a single integral ("nlinfit with modelfun as an integral"). and tried to modify it.
But I do not know how to pass over the integral limits.
%Given e.g. 3 volumes
%Unknown: peak and sigma for bell shape that fits best those volume parts
Vin=[0.1503 0.0949 0.0238];%Volume at xin,yin
xin=[0 1 2];%x-position of volume part under the 3D bell shape
yin=[0 0 0];%y-position of volume part under the 3D bell shape
parIn=[0.1632 1];%initial guess for peak and sigma of bell shape
%In this case it is also the result.
%Formula for bell volume, double integral. par(1)=peak, par(2)=sigma:
Integrand = @(par,x,y) par(1)*exp(-(x.^2+y.^2)/(2*par(2)^2));
Integralfun=@(par,x,y,xin,yin) integral2(@(x,y) Integrand(par,x,y),xin-0.5,xin+0.5,yin-0.5,yin+0.5);
parOut= nlinfit(xin,Vin,Integralfun,parIn);
Peter Seibold on 6 Jul 2020
Edited: Peter Seibold on 7 Jul 2020
I hope following explains better my problem. I am looking for P and σ.
A solution could be to minimize the sum of absolute differences of the given volumes Vin and the calculated volume-parts under the 3D bell shape. The example below is only for three Vin, but more Vin are possible.
Please observe that only the integral limits change with every Vin.
To simplify the problem, the y-limits remain constant.
The integral x-limits are:
xin+-0.5 with xin=0, 1, 2
How can I find P and σ by using nlinfit or lsqcurvefit?

Peter Seibold on 7 Jul 2020
Edited: Peter Seibold on 7 Jul 2020
I could not figure out how to solve the problem with 'nIinfit' , therefore, I solved the problem by using fminsearch’.
%Given e.g. 3 volumes
%Unknown: peak and sigma for bell shape that fits best those volume parts
Vin=[0.1503 0.0949 0.0238];%Volume at xin,yin
xin=[0 1 2];%x-position of volume part under the 3D bell shape
yin=[0 0 0];%y-position of volume part under the 3D bell shape
% parIn=[0.1632 1];%initial guess for peak and sigma of bell shape
parIn=[1 2];%initial guess for peak and sigma of bell shape
%Formula for bell volume, double integral. par(1)=peak, par(2)=sigma:
Integrand = @(par,x,y) exp(-(x.^2+y.^2)/(2*par(2)^2));
%Build the function for the sum of absolute differences:
funStr='@(par) ';
for i=1:numel(Vin)
is=num2str(i);%i as string
funStr=[funStr 'abs(par(1)*integral2(@(x,y) Integrand(par,x,y),'...
'xin(' is ')-0.5,xin(' is ')+0.5,yin(' is ')-0.5,yin(' is ')+0.5)-Vin(' is '))+'];
end
funStr=funStr(1:end-1);%remove last '+'
mystr2func=@(Str)evalin('base',Str);%dirty trick to include workspace variables
fun=mystr2func(funStr);%convert string to function
[parOut,fval,exitflag,output] = fminsearch(fun,parIn)
%parOut expected: [0.1632 1]
Peter Seibold on 7 Jul 2020
If your variables (Vin, xin,yin) are not in the base workspace, then replace 'mystr2func=@(Str)evalin('base',Str)' with
mystr2func=CreateMystr2func;
and add this function:
function mystr2func=CreateMystr2func
mystr2func=@(Str)evalin('caller',Str);

### Categories

Find more on Calculus 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!