nlinfit with modelfun as an integral2 and variable integral limits

1 view (last 30 days)
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);
  1 Comment
Peter Seibold
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?

Sign in to comment.

Accepted Answer

Peter Seibold
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]
  1 Comment
Peter Seibold
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);

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!