echo on
% == Least Squares case uitwerken ==

% De huis prijs grafiek
%     prijs woonoppervlakte, perceel
D = [ 375, 135, 156;
      399, 190,  82;
      349, 145,  85;
      390, 190,  77;
      399, 150, 141;
      349, 160, 132;
    ];
% De formule p = x1 * wo + x2 * po toepassen
A = D(:,2:3)
p = D(:,1)
Ap = [ A p ]
rref(Ap)
% Niets bruikbaars

% De kleinste kwadraten aanpak
AT = A'
ATA = AT * A
ATp = AT * p
ATAp = [ ATA ATp ]
T = rref(ATAp)

%Bepaal nu de geschatte prijzen
x = T(:,3)
pschati = A * x

% In dit model wordt de vraagprijs (x1) zwaarder gewogen en wel
percent = x(2) / x(1) * 100

%pschati prijs wo, po 
[ pschati D ]
% Huizen met een groot po worden fout gemeten
% Voor mijn eigen modelprijs formule ga ik krijgen naar 2de hands
% auto advertenties (van vanderplas.nl in autoplaza leiderdorps dagblad
% woensdag 22 mei 2007)
% Merksterkte heb ik de volgende waardes gegeven
% Mini = 11;
% BMW = 10;
% Audi = 5;
% Fiat = 1;
% prijs, Bouwjaar, Aantal kilometer x1000, aantal deuren, CD, airco
% merksterkte
% speler
Autos = [
  22950, 2003, 113, 5, 1, 1,  5;
  23950, 2004,  64, 5, 1, 1, 10;
  13950, 2001, 116, 3, 0, 1, 10;
  24950, 2003, 102, 5, 0, 1, 10;
  24950, 2002,  47, 2, 1, 1, 10;
   4750, 2001,  54, 3, 1, 0, 1;
  14950, 2005,  39, 3, 1, 1, 1;
  13950, 2002,  86, 3, 1, 1, 11;
  17750, 2003,  55, 3, 1, 1, 11;
  ];

% Ik laat alles even zwaar wegen
A = Autos(:,2:end -1)
p = Autos(:,1)

% De kleinste kwadraten aanpak
AT = A';
ATA = AT * A;
ATp = AT * p;
ATAp = [ ATA ATp ];
T = rref(ATAp);

%Bepaal nu de geschatte prijzen
x = T(:,end);
pschati = A * x;

tp = sum(p);
tsp = sum(pschati);
prijs = [ p pschati ];
for i=1:size(prijs,1)
    prijs(i,3) = prijs(i,2) / prijs(i,1) * 100;
end
% prijs, geschatte prijs
int32(prijs)
total = int32(sum(prijs(:,:)))
%De excessen
min(prijs(:,end)), max(prijs(:,end))
int32( sum(prijs(:,3) / size(prijs,1) ) )

% Dit presteert redelijk, als ik echter de excessen weghaal dan gaat hij
% een stuk beter. Ik zal de auto nader moeten onderzoeken waarom deze
% eruit springt

%
% Olympische 400 meter recordtijden
%

% het script lijn400.m maakt lijnbenadering van de recordtijden als
% functie van olympisch jaar de designmatrix X en waarneemvector y.
% 
% Ik heb een matrix gemaakt van de orginele data en daar een mapping mee
% gemaakt
% zie lijn400.m
lijn400

% Ik neem aan dat de symbolen XT betekenen matrix X getranponeerd.
XT = X';
XTX = XT * X;
XTy = XT * y;
XTXy = [ XTX XTy ];

% inverse
beta = inv(XTX) * XTy
% Reductie
beta = rref(XTXy)
% Matlab \ operator
beta = X\y
% Alle geven dezefde antwoorden, echter wel anders geformateerd

% voor de functie rt400 zie het bestand rt400.m
for i=1:size(X,2)
    rtwaarden_lin(i) = rt400(X(i,2), beta);
end

% plot van de functie
figure;
grid on;
hold on;
plot(X(:,2), y, '*');
plot(X(:,2), rtwaarden_lin, '-r');
hold off;
title('Olympische spelen, mannen 400 meter - lineare plot');
pause
close all hidden;

% voor het maken van de K designmatrix zie kwa400.m
% of gebruik iets als
% K = X;
% K(:,3) = K(:,2) .^ 2;
K = kwa400;
KT = K';
KTK = KT * K;
KTy = KT * y;
KTKy = [ KTK KTy ];

rref(KTKy)
beta = K \ y
% zelfde (maar niet in notatie)

t = 0:38;
rtwaarden_kwa = beta(1) + beta(2) * t + beta(3) * (t .^ 2);

% plot van de functie
figure;
grid on;
hold on;
plot(X(:,2), y, '*');
plot(X(:,2), rtwaarden_lin, '-r');
plot(t, rtwaarden_kwa, '-g');
hold off;
title('Olympische spelen, mannen 400 meter - lineare and kwadratische plot');
pause
close all hidden;

%
% Temperatuurverloop
%

% data invul
clear X;
t_vul;

% rtwaarden bepaling
beta_l = X \ y_l;
beta_h = X \ y_h;

t=0:365;
func_l = beta_l(1) + beta_l(2) * cos(( 2 * pi * t) / 365) + ...
                + beta_l(3) * sin(( 2 * pi * t) / 365);
func_h = beta_h(1) + beta_h(2) * cos(( 2 * pi * t) / 365) + ...
                + beta_h(3) * sin(( 2 * pi * t) / 365);

% plot van de functie
figure;
grid on;
hold on;
plot(orig(:,1), y_l, 'r*');
plot(orig(:,1), y_h, 'g*');
plot(func_l, '-r');
plot(func_h, '-g');
hold off;
axis ([0, 365, 0, 100]);
title('Gemiddelde temperature - kwadratische plot');
pause
close all hidden;