How to fit imposing some conditions?

Question

Fc Fc 2018-4-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/392537-how-to-fit-imposing-some-conditions

编辑： Matt J 2018-4-13

Hi every one, I've a matlab program. In the program there is a call to an external program for the fit function. This is the external program

function [estimates, model] = fitcurvedemo(xdata, ydata)
% Call fminsearch with a random starting point.
global ckgraconv_n ckdian
global A lambda sse R2
start_point = rand(1,2);%%scegli un valore casuale per tutti i parametri
model = @fun;
estimates = fminsearch(model, start_point);
      function [sse, FittedCurve] = fun(params)
          A = params(1);
          lambda = params(2);
         FittedCurve =( A.*ckgraconv_n +  lambda.*ckdian); %%picco gaussiano allargato
          ymean = mean(ydata);
          SS = ydata-ymean;
          SSt =sum(SS.^2);
          ErrorVector = FittedCurve - ydata;
          sse = sum(ErrorVector .^ 2);
          R2 = 1 - sse/SSt;
      end
  end

The y data are called ydata and the fit function is:

A.*ckgraconv_n + lambda.*ckdian

I must edit the fit so that the residual values are all positive, that is

(ydata-A.*ckgraconv_n - lambda.*ckdian)>0.

Then I must fit the data (finding A and lambda values) and imposing the condition

(ydata-A.*ckgraconv_n - lambda.*ckdian)>0.

I tried writing

function [sse, FittedCurve] = fun(params)
        params=fminsearch(model,ydata-params(1).*ckgraconv_n -  params(2).*ckdian>0);
        A = params(1);
        lambda = params(2);

but I get the error

Error using fminsearch (line 96)
FMINSEARCH accepts inputs only of data type double.
Error in fitcurvedemo/fun (line 11)
        params=fminsearch(model,ydata-params(1)*ckgraconv_n -  params(2)*ckdian>0);
Error in fminsearch (line 200)
fv(:,1) = funfcn(x,varargin{:});
Error in fitcurvedemo (line 8)
estimates = fminsearch(model, start_point);
Error in sp2sp3diff_Titantah (line 659)
[estimates, model] = fitcurvedemo(xdata,ydata);

Please, can someone help me? Best regards

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2018-4-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/392537-how-to-fit-imposing-some-conditions#answer_313492

在 MATLAB Online 中打开

It is a linear fit, so lsqlin would be appropriate

C=[ckgraconv_n(:) ,  ckdian(:)];
d=ydata(:);
Aineq=C;
bineq=d;
solution=lsqlin(C,d,Aineq,bineq);

12 个评论
显示 10更早的评论隐藏 10更早的评论

Matt J 2018-4-7

在 MATLAB Online 中打开

Fc Fc commented:

Hi thank you for your reply, but maybe I didn't understand it... In the original program, this is the call to the fit external program (the external fit program source is in my first topic):

[estimates, model] = fitcurvedemo(xdata,ydata);
[sse, FittedCurve] = model(estimates);

After the fit, it plots the data, the fit function and the residuals

plot(xdata, ydata, 'b*'); hold on  %%%DATA
%  plot(xdata, FittedCurve, 'b');
  axis([ xminas   xmaxas  yminas  ymaxas])
  plot(xdata,ydata - A.*ckgraconv_n-lambda.*ckdian,'r'); %RESIDUALS
  plot(xdata, A.*ckgraconv_n+lambda.*ckdian,'g');  %%FIT FUNCTION

and add the legend in the plot

legend ('Data','Residuals','Fit');
  dim=[0.7 0.1 0.3 0.3];
  str={'FIT',['X_{sp^{2}}= ' num2str(A) ''] ['X_{sp^{3}}= ' num2str(lambda) '']};
  annotation('textbox',dim,'String',str,'FitBoxToText','on');

Morover, it generates another plot about all residuals in a loop (analysing more data files)

figure(3); hold on
  plot3=figure (3);
set(plot3, 'Visible', 'off');
title('Energy- Residuals');
  axis([ xminas   xmaxas  yminas  ymaxas])
plot(xdata,ydata - A.*ckgraconv_n-lambda.*ckdian); %Residuals
xlabel('Energy (eV)');
  ylabel('Relative counts');
  legend_str{count}=nomespetanal;
  legend (legend_str);

Morover, it calculates the values percent of the fit results to print it on a file

A100=100*A;
  lambda100=100*lambda;
  sp2e=100*[1-(A+lambda)];
   Nosp2=100*[[1-(A+lambda)]+lambda];

It prints also the residuals values on a file

f2=fopen(path_fileresiduo,'w+');
  Residuo=[xdata; ydata - A.*ckgraconv_n-lambda.*ckdian];
  fprintf (f2, '%2.2f \t %2.4f \r\n', Residuo);

Lastly it prints the R2 values on a file (The R2 values are calculated in the external fit program)

   f3=fopen(path_filer2sp2,'a+');
   end
    Matsp2e=[R2; sp2e];
    fprintf (f3, '%2.2f \t %2.4f \r\n', Matsp2e);

It works but, as you can see by next plot, there are negative residuals (red line)

then I must impose that all residuals have to be positive:

ydata - A.*ckgraconv_n-lambda.*ckdian>0

By your reply I deleted the call to the external program and I wrote:

C=[ckgraconv_n(:) ,  ckdian(:)];
d=ydata(:);
Aineq=C;
bineq=d;
solution=lsqlin(C,d,Aineq,bineq);
plot(xdata, ydata, 'b*'); hold on  %%%DATA
  axis([ xminas   xmaxas  yminas  ymaxas])
  plot(solution,'g'); 
  plot(ydata - solution,'r'); %residuals

To start to generate the plot, but this is the generated plot

there is just the data plot without the fit and residuals ones. I must obtain this plot

Thank you

Fc Fc 2018-4-7

编辑：Fc Fc 2018-4-7

在 MATLAB Online 中打开

Hi matt and thanks. I deleted the call to the external program and I wrote (in the main program):

C=[ckgraconv_n(:) ,  ckdian(:)];
d=ydata(:);
Aineq=C;
bineq=d;
solution=lsqlin(C,d,Aineq,bineq);
A=solution(1);
lambda=solution(2);
FittedCurve=C*solution;
ymean = mean(ydata);
        SS = ydata-ymean;
        SSt =sum(SS.^2);
        FittedCurve=FittedCurve';
        ErrorVector = FittedCurve - ydata;
        sse = sum(ErrorVector .^ 2);
        R2 = 1 - sse/SSt;
plot(xdata, ydata, 'b*'); hold on
plot(xdata,ydata - A.*ckgraconv_n-lambda.*ckdian,'r'); 
plot(xdata, A.*ckgraconv_n+lambda.*ckdian,'g');

Now I get positive residuals,

but

1) unfortunately, if I check the A and lambda values are (called Xsp2 and Xsp3 in the plot), for example for the sample 1

A=36% lambda=37%

instead, they must be

A=65% lambda=12%

Do you know the reason?

2) I get negative values of R2 calculated as

R2 = 1 - sse/SSt;

3) By using lsqlin, at the end of the fit, matlab writes in the command window:

Minimum found that satisfies the constraints.
Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.
<stopping criteria details>

Is there a way do not show it?

Fc Fc 2018-4-8

编辑：Fc Fc 2018-4-8

在 MATLAB Online 中打开

1. Maybe I didn't explain good! I've a fellowship and for it, my supervisors gave me this program to modify (adding new functions) and to analyse some data. Before analysing new data, to be sure that the new program works fine, I must analys old data analysed in their scientific article. The original source that they gave me is based on the external program that I showed you. As I showed you, the original source gave negatives residuals, but my supervisors said that I had to get positive ones (they got positive residuals in their article, but unfortunately they didn't help me to get them!). By using your code, I finally got positive residuals, but unfortunately the lambda and A values are different than lambda and A values written in the article, so I just asked you if there are differences between the fit made by using the external program that I posted here and lsqlin matlab function!

2) R2 expression comes from

ymean = mean(ydata);
        SS = ydata-ymean;
        SSt =sum(SS.^2);
        FittedCurve=FittedCurve';
        ErrorVector = FittedCurve - ydata;
        sse = sum(ErrorVector .^ 2);
        R2 = 1 - sse/SSt;

In my supervisor's article, they just got positive values of R2 (and now I'm analysing the article data, so I must get same results!)

3) Thank you!

---------------------EDIT------------------------ Don't worry! I understood! A,lambda,R2 values also depend by other parameters (energy range, etc.)...this is the reason of differences! Thank you for your help!

P.s. Do you know the reason because using the lsqline I don't get negative residuals instead, using the external program I got them? What is the problem in the fitcurvedemo.m program?

Matt J 2018-4-13

As Walter said, just use the options input to select a different algorithm

https://www.mathworks.com/help/optim/ug/lsqlin.html#buchrig-1-options

Fc Fc 2018-4-13

编辑：Matt J 2018-4-13

Thank you, I will try it! Please, can you help me also in this https://it.mathworks.com/matlabcentral/answers/394080-calculating-and-plotting-fwhm please?

请先登录，再进行评论。

Answer 2

Walter Roberson 2018-4-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/392537-how-to-fit-imposing-some-conditions#answer_313872

在 MATLAB Online 中打开

You are passing in

ydata-params(1).*ckgraconv_n - params(2).*ckdian>0

In MATLAB, > has lower priority than the arithmetic operations, so this expression is

(ydata-params(1).*ckgraconv_n - params(2).*ckdian)>0

which is of data type logical.

If you want to check params(2).*ckdian and subtract 1 in the locations it is greater than 0, then you would use

ydata-params(1).*ckgraconv_n - (params(2).*ckdian>0)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 3

Jeff Miller 2018-4-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/392537-how-to-fit-imposing-some-conditions#answer_313875

在 MATLAB Online 中打开

If you want to keep the residuals positive, trying building in an extra penalty for negative residuals. In your original model fun, for example, you might use:

sse = sum(ErrorVector .^ 2) + 1000*sum(ErrorVector(ErrorVector<0) .^ 2);

Seems like that would quickly drive fminsearch away from solutions that give any negative residuals.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 4

Fc Fc 2018-4-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/392537-how-to-fit-imposing-some-conditions#answer_313892

编辑：Fc Fc 2018-4-11

在 MATLAB Online 中打开

Thank you both for your reply. @Walter

(ydata-params(1).*ckgraconv_n - params(2).*ckdian)>0

it is exactly what I must obtain. To understnand it, look my plot

Residuals are plotted by red line. As you can see, in the ranges [283;285] and [290.5;294] residuals are negative. Instead, I must obtain this plot

Here you can see that in all energy range, residuals (red line) are positive.

@Jeff I edited the fitcurvedemo.m file writing your code

sse = sum(ErrorVector .^ 2) + 1000*sum(ErrorVector(ErrorVector<0) .^ 2);

instead of sse = sum(ErrorVector .^ 2) but, as you can see by next plot, in this case I obtain negative residuals in most of the range.

Thanks

PS. About your comment

Seems like that would quickly drive fminsearch away from solutions that give any negative residuals.

To obtain positive residuals, It isn't my demand, but my supervisor ones...

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Jeff Miller 2018-4-7

Looks like I got the sign wrong and it should have been 1000*sum(ErrorVector(ErrorVector>0) .^ 2);

请先登录，再进行评论。

How to fit imposing some conditions?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

12 个评论
显示 10更早的评论隐藏 10更早的评论

更多回答（3 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

How to fit imposing some conditions?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

12 个评论 显示 10更早的评论隐藏 10更早的评论

更多回答（3 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

12 个评论
显示 10更早的评论隐藏 10更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论