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

I'm not sure that is not possible to use a simple fit 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:

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)

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)','problem','z');
[c2,gof2] = fit(xdata,ydata,f,'problem',z);
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('(0.1424/(1-0.0076*beta2+0.0076^2*beta3))*(1-beta2*x+beta3*x^2)');
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 fit 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\b
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];
jdata= [100; 200; 300; 350; 1000; 900];
z=ydata(5);
z2=xdata(5);
options = fitoptions('Weights',jdata);
f = fittype('(0.1424/(1-0.0076*beta2+0.0076^2*beta3))*(1-beta2*x+beta3*x^2)');
c2 = fit(xdata,ydata,f,options);
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:

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

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.