nlinfit with modelfun as an integral2 and variable integral limits

1 次查看(过去 30 天)
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 个评论
Peter Seibold
Peter Seibold 2020-7-6
编辑:Peter Seibold 2020-7-7
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
Peter Seibold 2020-7-7
编辑:Peter Seibold 2020-7-7
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 个评论
Peter Seibold
Peter Seibold 2020-7-7
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);

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Chemistry 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by