Could anyone help me to fit y=1/x on my data?

Hi,

could someone help me how I can fit y=1/x (y=constant/x) on my data with displaying equation?

clc
close all
clear all
hold on
grid on
set(gca,'fontweight','bold','fontsize',16);
xlim([0 3])
ylim([0 350])
x =[0.5,1,2]; 
y =[273.3333333,121.8287037,66.72297297]; 
e =[32.63764351,23.3035763,14.67708363]; 
errorbar(x,y,e ,'o b','MarkerSize',12, 'LineWidth',8);
x1 =[0.5,1,2]; 
y1 =[248.4821429,115.0833333,51.58333333]; 
e1 =[26.74351835,9.791376249,6.60976072]; 
errorbar(x1,y1,e1 , ' ^ b','MarkerSize',12, 'LineWidth',8);
x2 =[0.5,1,2]; 
y2 =[278.4895833,118.4151786,67.5]; 
e2 =[32.11137659,26.89490119,15.47805907]; 
errorbar(x2,y2,e2 , 'square b','MarkerSize',12, 'LineWidth',8);
x3 =[0.5,1,2]; 
y3 =[226.25,117.2678571,53.20833333]; 
e3 =[30.41523969,25.3971604,22.89566137]; 
errorbar(x3,y3,e3 , 'd b','MarkerSize',12, 'LineWidth',8);
x4 =[1,2]; 
y4 =[110.7954545,57.5]; 
e4 =[27.12006334,16.83560216]; 
errorbar(x4,y4,e4 , '* b','MarkerSize',12, 'LineWidth',8); 
hold off

Answer 1

[Mathieu NOE]

hello

here you are my friend

plot

code

clc
close all
clearvars

hold on
grid on
set(gca,'fontweight','bold','fontsize',16);
xlim([0 3])
ylim([0 350])
x =[0.5,1,2]; 
y =[273.3333333,121.8287037,66.72297297]; 
e =[32.63764351,23.3035763,14.67708363]; 
errorbar(x,y,e ,'o b','MarkerSize',12, 'LineWidth',8);
x1 =[0.5,1,2]; 
y1 =[248.4821429,115.0833333,51.58333333]; 
e1 =[26.74351835,9.791376249,6.60976072]; 
errorbar(x1,y1,e1 , ' ^ b','MarkerSize',12, 'LineWidth',8);
x2 =[0.5,1,2]; 
y2 =[278.4895833,118.4151786,67.5]; 
e2 =[32.11137659,26.89490119,15.47805907]; 
errorbar(x2,y2,e2 , 'square b','MarkerSize',12, 'LineWidth',8);
x3 =[0.5,1,2]; 
y3 =[226.25,117.2678571,53.20833333]; 
e3 =[30.41523969,25.3971604,22.89566137]; 
errorbar(x3,y3,e3 , 'd b','MarkerSize',12, 'LineWidth',8);
x4 =[1,2]; 
y4 =[110.7954545,57.5]; 
e4 =[27.12006334,16.83560216]; 
errorbar(x4,y4,e4 , '* b','MarkerSize',12, 'LineWidth',8); 

%% curve fit using fminsearch
% We would like to fit the function
% y = a + b./x
% to the data.
ym = mean([y;y1;y2;y3;[NaN y4]],'omitnan'); % create the mean y value for each x value

f = @(a,b,x) a + b./x;
obj_fun = @(params) norm(f(params(1), params(2),x)-ym);
sol = fminsearch(obj_fun, [mean(ym),-1]); % "good" initial guess => will give good fit 
a_sol = sol(1);
b_sol = sol(2);
y_fit = f(a_sol, b_sol, x);
Rsquared = my_Rsquared_coeff(ym,y_fit); % correlation coefficient

% plot(x, y_fit, '-',x,ym, 'r .', 'MarkerSize', 40)
plot(x, y_fit, 'r-')
title(['Fit a + b./x to data : R² =  ' num2str(Rsquared) ], 'FontSize', 20)
xlabel('x data', 'FontSize', 20)
ylabel('y data', 'FontSize', 20)
eqn = " y = "+a_sol + " + " +b_sol +" / x";
legend(eqn)
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Rsquared = my_Rsquared_coeff(data,data_fit)

    % R2 correlation coefficient computation
    sum_of_squares = sum((data-mean(data)).^2);% The total sum of squares
    sum_of_squares_of_residuals = sum((data-data_fit).^2);% The sum of squares of residuals, also called the residual sum of squares:
    Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;% definition of the coefficient of correlation is
    
end

[Walter Roberson]

Note that this is not pure constant/x : this is constant/x + constant . Which is often a better idea, but is not necessarily what you need.

[Walter Roberson]

You could look at the eqn variable, which describes it as a string, or you could ook at the sol variable, where sol(1) is the additive constant and sol(2) is the "c" value that x is to be divided into.

Answer 2

[Walter Roberson]

x =[0.5,1,2]; 
y =[273.3333333,121.8287037,66.72297297]; 
C = mean(x.*y)

C = 130.6471

predy = C./x;

syms c
residue = expand(sum((c./x - y).^2))

residue = 

bestc = solve(diff(residue,c),c)

bestc = 

vpa(bestc)

ans = 

133.68702033999999825720719638325

predy2 = double(bestc)./x;
plot(x, y, ':*', x, predy, '--.b', x, predy2, '-+r')
legend({'raw', 'mean', 'calculus'})

The calculus approach provides the theoretical optimum, and is valid for any number of points. If you do not have the symbolic toolbox, it is possible to create a formula that could be applied

The C value calculate from mean() should be a cheap approximation that is not too far off.

syms X [1 3]
syms Y [1 3]
syms c
residue = sum((c./X - Y).^2);
general_solution = solve(diff(residue,c),c)

general_solution = 

Obviously the 2 factors out, and the general formula becomes

c = sum(y./x)./sum(1./x.^2)