% REGRESIÓN LINEAL - Ecuación de Riedel
clear
% log(P) = A + B/T + C log(T) + D T^2
% supongamos que tenemos los siguientes datos experimentales
T = [290.1 302.4 311.2 318.7 325.1 330.5 334.9 338.9 343.0 346.2 349.9 353.5 356.2 358.9 362.3 365.2 367.9 370.5 373.2 375.8 378.5 381.3];
P = [8.6 15.4 22.5 30.5 39.0 47.6 55.5 63.8 73.0 81.3 91.3 102.0 110.9 120.1 132.8 144.5 155.8 167.6 180.0 193.5 207.9 223.4];

% transformación  y = A + B*x1 + C*x2 + D*x3
y = log10(P)'; x1 = 1./T; x2 = log10(T); x3 = T.^2;

X = [ones(length(T),1) x1' x2' x3'];

disp('coeficientes:  ')
b = inv(X'*X)*X'*y; b'

Pcalc = 10.^(X*b)';
% varianza
disp('Varianza =  ')
V = (P-Pcalc)*(P-Pcalc)'/(length(P)-4)

figure(1)
plot(T,Pcalc,T,P,'o')
xlabel('T (K)');ylabel('P (kPa)');title('Ajuste según ecuación de Riedel')

figure(2)
plot(T,P-Pcalc,'o',T,0)
xlabel('T (K)'); ylabel('P (kPa)');title('residuos')

% también se puede hacer un ajuste no lineal
function e = riedel(T,P,p)
    A = p(1);B = p(2);C = p(3); D = p(4);
    y = A + B./T + C*log10(T) + (D*T.^2);
    Pcalc2 = 10.^y;
    e = (P-Pcalc2)*(P-Pcalc2)';
endfunction
p0 = [53.5 -3540 -16.6 .0000077]; % es muy sensible al punto inicial
disp('coefs. con fminsearch')
p = fminsearch(@(p) riedel(T,P,p),p0)
    A = p(1);B = p(2);C = p(3); D = p(4);
    y = A + B./T + C*log10(T) + (D*T.^2); Pcalc2 = 10.^y;
figure(3)
    plot(T,P,T,Pcalc2,'r*')
disp('Varianza2 = ')
V2 = (P-Pcalc2)*(P-Pcalc2)'/(length(P)-4)