How to fit imposing some conditions?

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

Answer 1

[Matt J]

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

[Matt J]

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

[Matt J]

The output of

 solution=lsqlin(C,d,Aineq,bineq);

is the estimate of the two unknown parameters.

A=solution(1);
lambda=solution(2);

You must use these to generate the fitted curve

fittedcurve=C*solution;

and then use that to generate your plots.

[Fc Fc]

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?

[Matt J]

(1) instead, they must be A=65% lambda=12%. Do you know the reason?

How could I? How could you, for that matter? If you already know the solution, why are you going through all this trouble? In any case, to convince yourself, you can check if A=65% lambda=12% gives a lower value of the objective function.

(2) I get negative values of R2 calculated as R2 = 1 - sse/SSt

Where does that expression come from? Are you sure it applies when the residuals are constrained to be positive?

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

 optimoptions(@lsqlin, 'Display','off');

[Fc Fc]

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]

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?

Probably what Walter said in his response.

[Fc Fc]

Hi in the laboratory, even if the code runs, matlab gives this warning:

Warning: The trust-region-reflective algorithm can handle bound constraints only;
    using active-set algorithm instead. 
> In lsqlin (line 303)

Do you know how to fix it? Thank you

[Walter Roberson]

Pass an options structure that tells it to use active-set . Or remove the constraints other than the upper and lower bound.

[Matt J]

Hmmm. What version of MATLAB are you using? In recent versions, lsqlin does not have an active-set algorithm. In any case, the warning is innocuous. It's just letting you know it had to deviate from the default algorithm selection.

[Fc Fc]

Thank you to the both of you. In the laboratory there is matlab 2015 ). I know it's innocuous but I'm sure that my supervisor (he can't use matlab) will makes problems looking the warning!

[Matt J]

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]

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]

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)

Answer 3

[Jeff Miller]

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.

Answer 4

[Fc Fc]

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...

[Jeff Miller]

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