Problem with FMINCON, please help!

Hi,

I am trying to use FMINCON to solve a fitting problem. I cannot use the fitting toolbox since I want the fit to pass through a specific point, and after trying a lot I realized reading several posts that this is not possible.

However, I am not able to solve the problem with FMINCON either.
Here is what I need to do.

I have two column vectors

xdata= [-4.22; -3.64; -3.35; -2.79; -0.0076; 1.71]
ydata= [0.34; 0.28; 0.26; 0.25; 0.1424, 0.137]

and I need the fitting to pass through the point xdata(5,1) ydata(5,1) that is (-0.0076 , 0.1424).

I must use a custom objective function, which is given by

ydata= z*(1+beta(1)*xdata+beta(2)*xdata.^2)

where z=ydata(5,1)=0.1424

Could someone help me out on this?

Thanks a lot in advance

Comments

  • I'm not sure that is not possible to use a simple [u]fit[/u] function. I try to explain my idea, tell me if I am wrong.

    You have to find a fit function that depends on two variables (beta1 and beta2). y= 0.1424*(1+beta(1)*x+beta(2)*x.^2)
    You have a constraint: the fitting curve have to pass through a specific point (let's name it (xc yc)). This give you a condition on beta1 and beta2:
    yc= 0.1424*(1+beta(1)*xc+beta(2)*xc.^2)

    0.1424=0.1424*(1+beta(1)*(-0.0076)+beta(2)*(-0.0076).^2)

    so you can write beta1 as function of beta2. I try to compute manually (but probably I'm wrong) and I found: beta1=-xc*beta2.

    So, now, the initial function become:
    y= 0.1424*(1+(0.0076*beta2)*x+beta2*x.^2)

    You can use the fittype function to tell to matlab what is the function for fitting:

    f = fittype('0.1424*(1+(0.0076*beta1)*x+(beta1)*x^2)');

    it recognize by itself that beta1 is the free parameter and x the independent variable.

    Now, let's fit...
    [c2,gof2] = fit(xdata,ydata,f);
    and plot...
    figure()
    plot(xdata,ydata,'*')
    hold on
    plot(c2)

    I resume all the code:

    [code]xdata= [-4.22; -3.64; -3.35; -2.79; -0.0076; 1.71];
    ydata= [0.34; 0.28; 0.26; 0.25; 0.1424; 0.137];
    f = fittype('0.1424*(1+(0.0076*beta1)*x+(beta1)*x^2)');
    [c2,gof2] = fit(xdata,ydata,f);
    figure()
    plot(xdata,ydata,'*')
    hold on
    plot(c2)[/code]

    It runs and it gives me a good result. I hope it is what you want.

  • hi!

    first of all, thanks a lot for the help!

    I tried your code and it gives me the following result:
    In fit>handlewarn at 714
    In fit at 328
    ??? Index exceeds matrix dimensions.

    Error in ==> cfit.plot at 67
    ydata = ydata(idx);

    moreover, i do not know whether I can use this approach because in fittype('0.1424*(1+(0.0076*beta1)*x+(beta1)*x^2)') the first 0.1424 is actually something that I have to take from a variable since this is all a part of a code which will go inside a loop: basically I will have to do 1000 iterations where everytime xdata and ydata will change and everytime the first part (the 0.1424) will be taken as a cell inside ydata.
    so basically instead of the number I will have to put the name of the variable (which will contain a different number for every loop) and I don't know whether it will recognize that.



  • I correct the error (in bold).
    For the 0.1424 you can substitute it with z and initialize z in another part of the code:


    xdata= [-4.22; -3.64; -3.35; -2.79; -0.0076; 1.71];
    ydata= [0.34; 0.28; 0.26; 0.25; 0.1424; 0.137];
    z=ydata(5);
    f = fittype('z*(1+(0.0076*beta1)*x+(beta1)*x^2)');
    [c2,gof2] = fit(xdata,ydata,f);
    figure()
    plot(xdata,ydata,'*')
    hold on
    plot(c2)
  • ok now it works, but the problem is that, as I expected, it "fit" treats z as a variable and it gives a value and confidence bounds for it as well!

    moreover, the 0.0076 you have in your formula should itself be linked to a cell in ydata, so the program should be this:

    xdata= [-4.22; -3.64; -3.35; -2.79; -0.0076; 1.71];
    ydata= [0.34; 0.28; 0.26; 0.25; 0.1424; 0.137];
    z=ydata(5);
    z2=xdata(5);
    f = fittype('z*(1-(z2*beta1)*x+(beta1)*x^2)');
    [c2,gof2] = fit(xdata,ydata,f);
    figure()
    plot(xdata,ydata,'*')
    hold on
    plot(c2)

    and yes, it works, but if you look at in "c2" you'll see that z and z2 are treated as variables for the fitting procedure and the number that is written in them is not "read" by the optimizer. Am I wrong?

    Indeed on the help for fittype you find: By default, the independent variable is assumed to be x and the dependent variable is assumed to be y. All other variables are assumed to be coefficients.

    and I also tried to input the fact that I want only beta1 as the coefficient in "f = fittype('z*(1-(z2*beta1)*x+(beta1)*x^2)', 'coefficients',{'beta1'});"

    but I got an error message

    Thanks again
  • Yes, you are right, sorry.

    You have to specify it in the fit type (for the example I use only z as parameter):

    xdata= [-4.22; -3.64; -3.35; -2.79; -0.0076; 1.71];
    ydata= [0.34; 0.28; 0.26; 0.25; 0.1424; 0.137];
    z=ydata(5);
    z2=xdata(5);
    f = fittype('z*(1-(0.0076*beta1)*x+(beta1)*x^2)',[b]'problem','z'[/b]);
    [c2,gof2] = fit(xdata,ydata,f,[b]'problem',z[/b]);
    figure()
    plot(xdata,ydata,'*')
    hold on
    plot(c2)

    Tell me if it's ok. if you want to use also z2 you have to use a cell array, as the matlab help says.
  • Ok this works now, thanks a lot!

    However I now see that in this way I will not get my desired result since putting that point passing through in the function as well is not entirely correct.

    Let me explain myself: I still want the fitting function to pass through the point (y(5),x(5)) BUT the function in general form needs to be

    ydata= beta(1)*(1+beta(2)*xdata+beta(3)*xdata.^2)

    so that even if I insert one point I won't be able to use the same technique since I now have 3 unknowns. I could use two points, solve manually the linear system and to the same (but actually I am trying to do this without success since b1 multiplies everything), but I would rather prefer the curve to be constrained to go just through 1 point, so that it has more freedom. Is there any other way or now it is really possible only with fmincon (which I am not so good at)?
  • In this case...
    I solved the function, imponing the passing through the point, so that beta1 now is function of beta2 and beta3.
    I found:
    beta1=(0.1424/(1-0.0076*beta2+0.0076^2*beta3))

    Then, I modified the code:
    xdata= [-4.22; -3.64; -3.35; -2.79; -0.0076; 1.71];
    ydata= [0.34; 0.28; 0.26; 0.25; 0.1424; 0.137];
    z=ydata(5);
    z2=xdata(5);
    f = fittype('[b](0.1424/(1-0.0076*beta2+0.0076^2*beta3))*(1-beta2*x+beta3*x^2)[/b]');
    c2 = fit(xdata,ydata,f)
    figure()
    plot(xdata,ydata,'*')
    hold on
    plot(c2)

    In this way it find the coefficients beta2 and beta3.

    Something wrong?
  • Yes now it works perfectly, thanks a lot!!

    Last question: is there a way (which would be above my expectations) with this method to minimize not the mean squared error but a "volume weighted mean squared error", where obviously I would have another vector jdata with all the volumes corresponding to the cells in xdata and ydata, let's say jdata= [100; 200; 300; 350; 1000; 900].

    where volume weighted mean squared error is given by:

    (SumOveri(Volume(i)*(ydata(i)-ymodel(i))^2))/(Sumoveri(Volume(i)))
  • I think this is no possible with the [u]fit[/u] function.
  • Ok i think to have solved it with the backslash method:

    xdata= [-4.22; -3.64; -3.35; -2.79; -0.0076; 1.71];
    ydata= [0.34; 0.28; 0.26; 0.25; 0.1424; 0.137];

    %Volume
    jdata= [100; 200; 300; 350; 1000; 900];
    sum = 0.0
    for i = 1:size(jdata,1)
    sum = sum + jdata(i,1)
    end
    for i = 1:size(jdata,1)
    coeff(i,1) = sqrt(jdata(i,1)/sum)
    end
    % Create Design matrix
    for i = 1:size(xdata,1)
    A(i,1) = coeff(i,1)*(xdata(i,1)-xdata(5,1))
    A(i,2) = coeff(i,1)*(xdata(i,1)^2-xdata(5,1)^2)
    end
    % Create vector of the right hand side
    for i = 1:size(ydata,1)
    b(i,1) = coeff(i,1)*(ydata(i,1)-ydata(5,1))
    end
    % Solve linear system A*beta = b to get regression
    % coefficients
    % !MATLAB command
    beta = A
    beta1=ydata(5,1)-beta(1,1)*xdata(5,1)-beta(2,1)*xdata(5,1)^2;
    beta2=beta(1,1)/beta1;
    beta3=beta(2,1)/beta1;

    however now i have a problem: in this way I do not get the values for

    sse
    r-squared
    adj r-squared
    rmse

    as i would get with the fit function. how can i calculate them manually?

  • I don't know if you read your mailbox. If you have a vector with weights you can use it in the fit function. So, have you a vector with the weights?
  • This is the code:

    xdata= [-4.22; -3.64; -3.35; -2.79; -0.0076; 1.71];
    ydata= [0.34; 0.28; 0.26; 0.25; 0.1424; 0.137];
    [b]jdata= [100; 200; 300; 350; 1000; 900];[/b]
    z=ydata(5);
    z2=xdata(5);
    [b]options = fitoptions('Weights',jdata);[/b]
    f = fittype('(0.1424/(1-0.0076*beta2+0.0076^2*beta3))*(1-beta2*x+beta3*x^2)');
    c2 = fit(xdata,ydata,f,[b]options[/b]);
    figure()
    plot(xdata,ydata,'*')
    hold on
    plot(c2)
  • that's perfect!!! thanks a lot!!
  • To specify the starting point you have to modify a bit the code:

    [b]options1[/b] = fitoptions('Weights',jdata);
    f = fittype('(0.1424/(1-0.0076*beta2+0.0076^2*beta3))*(1-beta2*x+beta3*x^2)');
    [b]options.StartPoint=[1 1];[/b]
    c2 = fit(xdata,ydata,f,[b]options1[/b]);

    I modified options in options1 because options is a reserved word and now that you have to use it give you some problems.
    The code above runs well.
Sign In or Register to comment.

Howdy, Stranger!

It looks like you're new here. If you want to get involved, click one of these buttons!

Categories

In this Discussion