nlinfit with modelfun as an integral2 and variable integral limits

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]

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?

Answer 1

[Peter Seibold]

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]

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);