# Gyakorlat 2019. 05. 08.
clear
close all
output_precision(7)

disp("Gyakorlat 2019. 05. 08.");

disp("\n\
Legkisebb negyzetek modszere (Least Squares)")

disp("\n\n\
Alap problema:\n\
- y = f(x; a0, a1, ..., ak) parameterfuggo fuggveny\n\
- adott n db meres, azaz (x1, y1) ... (xn, yn) (kulonbozo x-ek)\n\
  melyekre yi = f(xi)")

disp("\n\n\
Legkisebb negyzetek: az n meresi ponthoz (xi, yi) meghatarozzuk\n\
azokat az a0, a1, ..., ak parametereket amire\n\
gyok( (y1 - f(x1; a0, ... ak))^2 + ... + (yn - f(xn; a0, ... ak))^2 ) minimalis.")

disp("\n\n\
Nyilvan a negyzetgyok nem szamit, igy\n\
(y1 - f(x1; a0, ... ak))^2 + ... + (yn - f(xn; a0, ... ak))^2\n\
minimalizalasa a cel. Ez a hibaknak a negyzeteinek osszege.\n\
Ezt az osszeget minimalizaljuk, azaz a negyzetes hibak legkisebb\n\
osszeget keressuk, innen a modszer neve: legkisebb negyzetek.\n\n\
A modszer altalanosan hasznalhato\n\
parameterfuggo f-re.")

disp("\n\
Az y = f(x; a0, a1, ..., ak) osszefugges a megfigyelt jelenseg modellje.\n\
A celunk a modell parametereinek (a0, ..., ak) meghatarozasa\n\
a mereseink (x1, y1) ... (xn, yn) alapjan.\n\n\
Nyilvanvaloan ha a parameterek szamat (k + 1) nezzuk, akkor legalabb\n\
ennyi meresre van szukseg, hogy egyaltalan remeny legyen ezek meghatarozasara,\n\
hiszen egy meres yi ~ f(xi; a0, a1, ..., ak) egy egyenletet eredmenyez.")

disp("\n\
Az elobbi egyenletben ~ jel szerepel = helyett. Ez azt jelkepezni,\n\
hogy a mereseket a gyakorlatban nem tudjuk pontosan elvegezni,\n\
azaz a meresek zajosak. Matematikailag:\n\
yi = f(xi; a0, a1, ..., ak) + zaj\n\n\
A zaj milyensege, esetleges kapcsolat ket kulonbozo meresben eszlelt\n\
zaj kozott fontos jellemzoje a problemanak. Mi most azt tesszuk fel, hogy\n\
a zaj meresenkent egymastol fuggetlen, N(0, sigma^2) eloszlast kovet.\n\n\
Mivel nem tanultatok meg valoszinusegszamitast/statisztikat, igy\n\
roviden errol annyit, hogy a normalis eloszlas nagyon fontos szerepet\n\
jatszik a termeszetben, surusegfuggvenyet mar biztos lattatok, ezt\n\
szoktak haranggorbenek, Gauss gorbenek nevezni, 1. abra.")

# !!!!!!!!!!!!!!
# Kod az abrazolashoz, nem szukseges megerteni

figure (1)
x = -5 : 0.01 : 5;

# Normalis eloszlas N(mu, sigma^2) surusegfuggvenyenek abrazolasa
# mu: a veletlen zaj varhato erteke: 0
mu = 0;
# sigma: a szoras (~ mennyire ter el a varhato ertektol a zaj)
sigma = 1;

y = 1 / sqrt( 2 * pi * sigma^2 ) * exp( - (x - mu).^2 / (2 * sigma^2) );
plot(x, y)
hold on
grid on
line( [sigma, sigma], [0, 1], "color", "k")
line( [2*sigma, 2*sigma], [0, 1], "color", "r")
line( [3*sigma, 3*sigma], [0, 1], "color", "g")

line( [-sigma, -sigma], [0, 1], "color", "k")
line( [-2*sigma, -2*sigma], [0, 1], "color", "r")
line( [-3*sigma, -3*sigma], [0, 1], "color", "g")

title("N(mu, sigma) surusegfuggvenye, mu = 0, sigma = 1 esetben")
legend("1/sqrt(2 \pi sigma^2) * e^{( -(x-mu)^2 / 2sigma )}", "+/- sigma: 68%", "+/- 2*sigma: 95%", "+/- 3*sigma: 99.7%", "location", "northeastoutside")

# Kod vege
# !!!!!!!!!!!!!!


disp("\n\
Zaraskent a normalis eloszlashoz: a kepen latjatok, hogy kb mit jelent a szoras.\n\
Ha sokszor merjuk csak magat a zajt es a kapott zajok erteket atlagoljuk,\n\
akkor egy 0 koruli szamot kapunk, a zaj varhato erteke 0, a zaj atlagos erteke 0.\n\n\
Viszont szinte sose 0 a zaj, hanem 0 korul valamilyen veletlen szam.\n\
Sok meres eseten azt fogjuk latni, hogy a zajok 68%-a +/- 1 szorasnyi\n\
tavolsagra van az atlagos ertektol. A zajok 95%-a +/- 2 szorasnyira van\n\
az atlagos ertektol. Vegul klasszikusan meg szoktak emliteni, hogy\n\
a zajok gyakorlatilag +/- 3 szorasnyira helyezkednek el, hiszen ez a\n\
tartomany mar 99.7%-ban tartalmazza oket.")

disp("\n\
Gyakorlatban igyekeznek minel tobb merest elvegezni ami\n\
meg gazdasagosan (meresi ido, kotlseg, human eroforras,\n\
a megfigyelt jelenseg idotartama) kivitelezheto. Azaz \n\
altalaban n joval nagyobb mint k:  n >> k.\n\n\
A sok (jo) meres altal a zaj hatasa a megoldasra csokken.")


disp("\n\
Megjegyzesek:\n\
- Ha rendelkezunk mogottes informaciokkal arrol, hogy\n\
  milyen a zaj, kevesebb meresi adatbol is lehet kovetkezteteseket levonni.\n\
- Ha a zaj fugg a merestol is pl nagyobb x-ekre nagyobb a zaj, kisebbekre kisebb,\n\
  akkor ezt fontos figyelembe venni kulonben teves eredmenyeket kaphatunk.\n\
  Ha ilyen problemaval szembesultok a gyakorlatban keressetek ra a\n\
  weighted least squares eljarasra.\n\
- Elofordul, hogy mar az eppen jelenlevo x ertekben sem vagyunk biztosak,\n\
  pl az is egy masik zajos meres eredmenye, akkor kicsit valtozik a\n\
  hibak negyzetosszege, ezzel most nem foglalkozunk.\n\n\
  A 2-3. abran szemleltetjuk, hogy milyen hibak negyzetet\n\
  minimalizaljuk, ha x nem zajos (ezzel dolgozunk ma, 2. abra),\n\
  illetve ha x zajos (3. abra).")
  
# !!!!!!!!!!!!!!
# Kod az abrazolashoz, nem szukseges megerteni

meresekSzama = 10;
x = linspace(0, 10, meresekSzama);
y = 1.1 * x + randn(1, meresekSzama) * 3;
egyutthatok = polyfit(x, y, 1);

xAbrazolas = -3 : 0.1 : 13;

figure(2)
plot(x, y, "b+")
axis("equal")
hold on 
plot(xAbrazolas, polyval(egyutthatok, xAbrazolas), "r")
for i = 1 : meresekSzama
  line([x(i) x(i)], [y(i) polyval(egyutthatok, x(i))], "color", "g")
endfor
title("Legkisebb negyzetek: a zold hibak negyzetenek osszegenek minimalizalasa")
legend("(x, y) meresek, x pontos, y zajos", "legjobb negyzetek egyenes modellje", "hibak", "location", "northeastoutside")

figure(3)
iranyVektor = [1 polyval(egyutthatok, 1)] - [0 polyval(egyutthatok, 0)];
iranyVektor = iranyVektor / norm(iranyVektor);
normalVektor = [iranyVektor(2) -iranyVektor(1)];
plot(x, y, "b+")
axis("equal")
hold on 
plot(xAbrazolas, polyval(egyutthatok, xAbrazolas), "r")
for i = 1 : meresekSzama
  # (xi, yi) tavolsaga az egyenestol
  tavolsag = dot([x(i) y(i)] - [0 polyval(egyutthatok, 0)], normalVektor); 
  # metszespont
  metszespont = [x(i) y(i)] - tavolsag * normalVektor;
  # vonal
  line([x(i) metszespont(1)], [y(i) metszespont(2)], "color", "g")
endfor
title("Legkisebb negyzetek: a zold hibak negyzetenek osszegenek minimalizalasa")
legend("(x, y) meresek, x zajos, y zajos", "egyenes modell", "hibak", "location", "northeastoutside")

# Kod vege
# !!!!!!!!!!!!!!
  
disp("\n\
Mi csak erintolegesen foglalkozunk ezzel a temaval. Mint fentebb irtuk,\n\
altalanosan parameterfuggo modell eseten lehet alkalmazni a minimum keresest\n\
ahol a celfuggveny (amit minimalizalunk) az a hibak negyzetosszege.\n\n\
Specialisan erdekes amikor a f egy fuggvenyosztaly alapelemeinek, baziselemeinek\n\
egy kombinaciojakent all elo:\n\n\
f(x; a0, ..., ak) = a0 * F0(x) + ... + ak * Fk(x)\n\n\
igy tulajdonkeppen azt keressuk, hogy melyek a bazisfuggvenyek F0(x)...Fk(x)\n\
egyutthatoi ebben a kombinacioban amire a negyzetes hiba a legkisebb.\n\n\
Meg specialisabb eset (mi csak ilyennel foglalkozunk):\n\
a bazisfuggvenyek: 1, x, x^2, x^3, ..., x^k, tehat\n\
f(x; a0, ..., ak) = ak x^k + ... + a1 x + a0 egy k-adfoku polinom.\n\
A legkisebb negyzetek modszere pedig azt a k-adfoku polinomot keresi meg,\n\
ami az n db meresunket a lehet legjobban kozeliti.")

###########################################################################
disp("\n\n\
=============================================\n\n\
1. Pelda\n\n\
Mert adatok (ezeket mindig celszeru abrazolni, hogy lassuk mivel van dolgunk - 4. abra):")
x = [1 2 -1 0 3]
y = [2.01 3.98 -2.1 0.001 5.97]

figure (4)
plot(x, y, "b+") # _b_lue +, azaz kek keresztekkel csak a pontokat rajzolja, s nem koti oket ossze
title("1. Pelda: meresi adatok")
legend("(x, y) meresi pontok", "location", "northeastoutside")


disp("\n\
Latszik, hogy y ~ 2x, keressuk meg hat azt az f(x; a) = a * x-et\n\
amire a negyzetes hiba a legkisebb.\n\
Maskepp fogalmazva a megfigyelt jelensegre a modellunk:\n\
y = a * x + zaj\n\
Kerdes: az elvegzett 5 meres alapjan mennyi lehet az ""a"" szorzo erteke?")

disp("\n\
Erre az egyszeru peldara megmutatjuk, hogyan jon ki ez a celfuggveny\n\
derivalasabol (kulcsszavak, linearis legkisebb negyzetek, normalfuggveny):\n\n\
Egy adott a eseten a hibak negyzetosszege:\n\
E(a) = Szumma( (yi - a*xi)^2 )\n\n\
Ha E(a) minimalis, ott a derivaltja 0 (visszafele ugye nem igaz az implikacio!),\n\
tehat derivaljunk egyet (osszeg derivaltja, osszetett fgv derivaltja, negyzetfgv derivaltja):\n\
E'(a) = Szumma( 2 * (yi - a * xi) * xi ) = \n\
      = 2 * Szumma( xi * yi ) - 2 * a * Szumma( xi * xi ) = \n\
      = 2 * [ x1*y1 + ... + xn*yn - a * (x1^2 + ... + xn^2) ]\n\n\
Igy E'(a) = 0 akkor teljesul, ha\n\
a = (x1*y1 + ... + xn*yn) / (x1^2 + ... + xn^2) \n\n\
Ennek a tortnek a szamlaloja az x es y vektorok skalaris szorzata,\n\
a nevezoje x-nek onmagaval vett skalaris szorzata, vagy maskepp\n\
szolva x normajanak a negyzete.\n\n\
A skalaris szorzatot kiszamithatjuk dot(x,y) paranccsal\n\
(angolul: dot product), vagy tudva, hogy a skalaris szorzat nem mas,\n\
mint x es y mint matrixok szorzata\n\
(ahol x 1*n-es sorvektor es y n*1-es oszlopvetkor)\n\
akar ugy is mint x * transzponalt(y)\n\
(hiszen itt most x-et es y-t mint sorvektort vittuk be, igy\n\
 az y sorvektor transzponaltja az a megfelelo oszlopvetkor)")
 
x
y

disp("\n\
y transzponaltja: transpose(y)")
transpose(y)
disp("vagy y'")
y'

disp("\n\
A skalaris szorzatok:\n\
1) x * y:")
dot(x,y)

disp("vagy mint matrixok szorzata:")
x * y'

disp("\n\
2) x * x:")
dot(x,x)

disp("vagy mint matrixok szorzata:")
x * x'

disp("\n\
Igy a legkisebb negyzetes kozelitessel kapjuk, hogy a legjobb a*x alaku\n\
fuggveny az 5 pontbol allo meresi adatunkra az")
a = dot(x,y) / dot(x,x)

disp("\n\
5. abra: a meresi pontjaink es a kalibralt a*x modellunk.")

figure (5)
plot(x, y, "b+")
hold on
plot(x, a * x, "r") # r_ed azaz piros vonal
title("1. Pelda: meresi adatok + kalibralt modell")
legend("(x, y) meresi pontok", "a * x", "location", "northeastoutside")

###########################################################################

disp("\n\n\
=============================================\n\n\
2. pelda\n\n\
Most megunknak generalunk egy a*x jellegu zajos merest, hogy jobban lassuk\n\
a modszert. Legyen a = 2 es vegezzunk el 100db merest az x = 30 ... 40\n\
tartomanyban. Azt feltelezzuk, hogy a meres zaja N(0,1) tipusu, azaz\n\
0 varhato erteku 1 szorasnegyzetu normalis eloszlast kovet (1. abra)")

meresekSzama = 100

disp("\n\
x tartomany: x = linspace(30, 40, meresekSzama)")
x = linspace(30, 40, meresekSzama); # 30-tol 40-ig meresekSzama db meresi pont egyenletes lepeskozzel

disp("\n\
Minden mereshez tartozik egy zaj: zaj = randn(1, meresekSzama)")
zaj = randn(1, meresekSzama); #rand -> random (veletlen) | n -> normalis | randn -> veletlen normalis eloszlasu szamokat general
# randn(n,m) egy n x m -es matrixot ad aminek az elemei mind veletlen szamok amik az N(0, 1) eloszlast kovetik
# probaljatok ki!

disp("\n\
Megjegyzes: mivel itt tenyleg veletlen zajt generaltok, ha tobbszor futtatjatok\n\
ezt az orat bemutato programot, akkor mindig kulonbozo eredmenyeket fogtok kapni!")

disp("\n\
zajos meresek: y = 2 * x + zaj")
y = 2 * x + zaj;

disp("\n\
Meresi adatokat abrazoljuk 6. abra")

figure (6)
plot(x, y, "b+")
title("2. Pelda: meresi adatok")
legend("(x, y) meresi pontok", "location", "northeastoutside")

disp("\n\
Legkisebb negyzetek modszere a legjobban illeszkedo a * x meghatarozasara:\n\
a = dot(x, y) / dot(x, x)  \n\n\
(ne lepodjetek meg, minden futtatasnal mas, hisz mindig uj veletlen merest kaptok!")
a = dot(x, y) / dot(x, x)

disp("\n\
7. abran nezzuk a mereseket, a kalibralt, azaz a legjobban illeszkedo modellt\n\
es az idealizalt (y = 2 * x) modellt.\n\n\
Vegyetek eszre, hogy a legjobban illeszkedo modell gyakrolatilag sosem fog\n\
egybeesni az idealizalt modellel. Ez nem problema, ez a zajos meres eredmenye.\n\
Nagyon sok merest elvegezve a zaj jelentosege csokken es a legjobban illeszkedo\n\
modell remelhetoleg kozel lesz az idealizalthoz.")

figure (7)
plot(x, y, "b+")
hold on 
plot(x, a * x, "r")
plot(x, 2 * x, "g")
title("2. Pelda: meresi adatok + modellek")
legend("(x, y) meresi pontok", "legjobb modell", "idealizalt modell", "location", "northeastoutside")

disp("\n\
Nagyitsatok ra a grafikonra az 7. abran, hogy lassatok elter a legjobb es\n\
az idealizalt modell (ha veletlen a = 2 adodna a legjobbra,\n\
futtasatok ujra, masodjara kizart, hogy ismet olyan legyen a veletlen zaj).")

###########################################################################

disp("\n\n\
=============================================\n\n\
3. pelda\n\n\
Terjunk vissza az 1. Pelda adataihoz")
x = [1 2 -1 0 3]
y = [2.01 3.98 -2.1 0.001 5.97]

disp("\n\
Bar ""ranezesre"" y = 2*x latszik, s igy iment y ~ a * x formaban\n\
kerestuk a megoldast kezzel, keressuk most y ~ a1 * x + a0 = f(x; a0, a1)\n\
formaban, azaz keressuk a legjobb elsofoku polinomot, azaz a legjobban\n\
illeszkedo egyenest!")

disp("\n\
Hasonloan mint fentebb derivalassal ki lehet hozni a normalegyenleteket,\n\
errol szakirodalomban reszletes leirast talaltok. Elokerulnek megint a\n\
Vandermonde matrixok, a rovidebb formula eseten a pszeudoinverz fogalma\n\
es meg sok minden mas.\n\n\
Mi ezt atugorjuk. Az Octave-nak van beepitett polinom illeszto\n\
modulja, ezt fogjuk hasznalni. Ez a hatterben a ma targyalt legkisebb\n\
negyzetes illesztest alkalmazza.\n\n\
egyutthatok = polyfit(x, y, 1)\n\n\
poly -> polinom\n\
fit  -> illesztes\n\
1    -> elsofoku, azaz egyenes: a1 * x + a0")
egyutthatok = polyfit(x, y, 1)

disp("\n\
Azaz a legjobban illeszkedo egyenes egyenlete:\n\
y = 2.0119 * x - 0.0397")

disp("\n\
Ha megvannak az egyutthatok, akkor hasznalhatjuk akar az Octave\n\
polyval(egyutthatok, x) parancsat, hogy az egyutthatok altal\n\
meghatarozott polinomot kiertekeljuk az x helyen:\n\n\
polyval(egyutthatok, 1)")
polyval(egyutthatok, 1)

disp("\n\
polyval(egyutthatok, 1.3)")
polyval(egyutthatok, 1.3)

disp("\n\
Hasonloan az eddig elmondottakhoz, ne lepodjunk meg, hogy nem az\n\
y ~ 2 * x + 0 adodik.\n\n\
Azon se lepodjunk meg, hogy nem is\n\
y ~ 1.998667 * x + 0 adodik (a = 1.998667 jott ki az 1. Peldaban)\n\
ez termeszetes, hiszen most ket parametert valaszthatunk meg, a1-et es a0-t\n\
igy tobb lehetosegunk van, igy minden bizonnyal a negyzetes hibaosszeg is kisebb.")

a = dot(x,y) / dot(x,x);

figure (8)
plot(x, y, "b+")
hold on 
plot(x, polyval(egyutthatok, x), "r")
plot(x, a * x, "g")
plot(x, 2 * x, "k")
title("3. Pelda: meresi adatok + modellek")
legend("(x, y) meresi pontok", "legjobb elsofoku modell a1*x + a0", "legjobb a*x modell", "idealizalt modell", "location", "northeastoutside")

disp("\n\
Vizualisan ezt nehez erzekelni (8. abra) a nagyjabol 2*x tendenciat koveto\n\
keves meresi pont kozott, szamoljuk ki mik a hibak!")

disp("\n\
Az (x, y) meresi pontokra az y = 2 * x modell hibainak negyzetosszege:")
# Szumma( (yi - 2 * xi)^2 )
hiba = 0;
# a hiba valtozoba gyujtjuk a hibainak negyzeteinek osszeget
for i = 1 : length(x)
  hiba = hiba + ( y(i) - 2 * x(i) )^2;
endfor
# kiirjuk az eredmenyt
hiba

disp("\n\
Az (x, y) meresi pontokra az y = 1.99867 * x modell (a * x) hibainak negyzetosszege:")
# Szumma( (yi - 1.99867 * xi)^2 )
hiba = 0;
# a hiba valtozoba gyujtjuk a hibainak negyzeteinek osszeget
for i = 1 : length(x)
  hiba = hiba + ( y(i) - 1.99867 * x(i) )^2;
endfor
# kiirjuk az eredmenyt
hiba

disp("\n\
Az (x, y) meresi pontokra az y = 2.0119 * x - 0.0397 modell hibainak negyzetosszege:")
# Szumma( (yi - (2.0119 * xi - 0.0397))^2 )
hiba = 0;
# a hiba valtozoba gyujtjuk a hibainak negyzeteinek osszeget
for i = 1 : length(x)
  hiba = hiba + ( y(i) - (2.0119 * x(i) - 0.0397) )^2;
endfor
# kiirjuk az eredmenyt
hiba

disp("\n\
Latszik, hogy 1.99867 * x jobban fekszik a meresi pontokra mint 2 * x\n\
illetve a 3 opcio kozul egyertelmuen 2.0119 * xi - 0.0397 illeszkedik\n\
a legjobban.")

disp("\n\n\
Megjegyzes: ha noveljuk az illeszkedo polinom fokszamat, elobb utobb\n\
a hibak negyzetosszege 0 lesz! Ez nem kivanatos jelenseg!!\n\n\
Emlekezzetek a polinomialis interpolaciora (Lagrange interpolaciora)\n\
ha k = n - 1, azaz (n - 1)-ed foku polinomot illesztunk az n db \n\
meresi pontra, akkor a legjobban illeszkedo polinom minden meresi ponton\n\
at fog menni (mar ha az xi-k mind kulonbozoek).\n\n\n\
Gyakorlati pelda mondjuk egy kilott lovedek, elrugott labda magassaganak a\n\
kovetese lehetne s ezaltal a gravitacios egyutthato becslese.\n\
Tudjuk, hogy a mozgasegyenlet masodofoku polinom, igy ha vegzunk 100 merest mig\n\
a targy a levegoben halad arra masodfoku polinomot illesztunk.\n\
Egy 99-ed foku polinom illesztese ertelmetlen volna es a \n\
celunkhoz sem jutnank kozelebb, azaz a gravitacios konstans\n\
meghatarozasahoz.")

disp("\n\n\
Demonstralas a jelen peldara:\n\n\
Most 5 darab meresi pontunk van, tegyuk fel, hogy tudjuk, hogy\n\
a fizikai/biologiai/kemiai stb modell ami ezt zajosan generalta\n\
az y = a * x + zaj alaku.\n\
Erre kiszamoltuk, hogy a meresi pontjaink alapjan y ~ 1.998667 * x \n\
a leheto legjobb modell.\n\n\
5 meresi pontra egy 4-ed foku polinom mar tokeletesen illeszkedik:\n\
egyutthatok4 = polyfit(x, y, 4)")
egyutthatok4 = polyfit(x, y, 4)

disp("\n\
Az (x, y) meresi pontokra illesztett 4-edfoku polinom hibainak negyzetosszege:")
# Szumma( (yi - polyval(egyutthatok4, xi)^2 )
hiba = 0;
# a hiba valtozoba gyujtjuk a hibainak negyzeteinek osszeget
for i = 1 : length(x)
  hiba = hiba + ( y(i) - polyval(egyutthatok4, x(i)) )^2;
endfor
# kiirjuk az eredmenyt
hiba

disp("\n\
1.18 * 10^-30 az gyakorlatilag 0...")

disp("\n\
Abrazoljuk a meresi pontokat, a legjobban illeszkedo a*x modellt\n\
illetve a tokeletesen illeszkedo 4-edfoku polinomot: 9. abra.\n\
Hogy lassuk a pontos tendenciat abrazoljuk az illeszteseket\n\
-3 es 8 kozott.")

xAbrazolas = -3 : 0.01 : 8;

figure (9)
plot(x, y, "b+")
hold on 
plot(xAbrazolas, a * xAbrazolas, "r")
plot(xAbrazolas, polyval(egyutthatok4, xAbrazolas), "k")
title("3. Pelda: meresi adatok + modellek")
legend("(x, y) meresi pontok", "legjobb a*x modell", "legjobb 4-edfoku modell (Lagrange interpolacio)", "location", "northeastoutside")

disp("\n\
Latszik, hogy hiaba illeszkedik a 4-edfoku polinom pontosan\n\
minden (zajos) meresi pontra, ha nekunk az y ~ a * x\n\
jelenseget kell megfigyelni, a 4-edfoku altal adott eredmenyek\n\
nem segitenek, sot kifejezetten rosszak tudnak lenni.")

###########################################################################

disp("\n\n\
=============================================\n\n\
4. pelda\n\n\
Az eddigi peldakat kombinaljuk. Vegyunk fel 10 db zajos merest\n\
az x = 30 ... 40 intervallumban az y = 2 * x + zaj modellre.")
meresekSzama = 10

disp("\n\
x tartomany: x = linspace(30, 40, meresekSzama)")
x = linspace(30, 40, meresekSzama); # 30-tol 40-ig meresekSzama db meresi pont egyenletes lepeskozzel

disp("\n\
Minden mereshez tartozik egy zaj: zaj = randn(1, meresekSzama)")
zaj = randn(1, meresekSzama); 

disp("\n\
zajos meresek: y = 2 * x + zaj")
y = 2 * x + zaj;

disp("\n\
10. abra: meresi pontok")

figure (10)
plot(x, y, "b+")
title("4. Pelda: meresi adatok")
legend("(x, y) meresi pontok", "location", "northeastoutside")

disp("\n\
Tegyuk fel, hogy nem tudjuk, hogy a modell a * x alaku,\n\
a pontokbol akkor is latszik, hogy nagyjabol linearis tendenciat\n\
kovetnek a meresek, igy a legjobban illeszkedo egyenest keressuk:\n\n\
egyutthatok = polyfit(x, y, 1)")
egyutthatok = polyfit(x, y, 1)

disp("\n\
10 meresi pontra egy 10 - 1 = 9 foku polinom teljesen illeszkedne!\n\
(itt az x-ek kulonbozoek igy nem okozhat problemat, hogy egy x-et\n\
ketszer is megmerunk zajosan es kulonbozo erteket kapunk)\n\n\
egyutthatok9 = polyfit(x, y, meresekSzama - 1)\n\n\
(ezt nem iratjuk ki mert a 10db egyutthato hosszu lenne.\n\n\
- jo esellyel ha lefuttatjatok itt egy figyelmeztetes (warning) fog\n\
  megjelenni, hogy a polinom illesztesnel numerikus bizonytalansag lepett fel.\n\
  Ezen nem kell most fennakadni.")
egyutthatok9 = polyfit(x, y, meresekSzama - 1);

disp("\n\
A 11. abran latjuk, hogy hiaba illeszkedik tokeletesen a 9-edfoku\n\
modell a meresi pontokon, ha a jelenseg fizikai modellje valojaban linearis,\n\
akkor a 9-edfoku polinom ertekei mas helyeken teljesen felrevezetoek.")

xAbrazolas = 29.5 : 0.01 : 40.5;

figure (11)
plot(x, y, "b+")
hold on 
plot(xAbrazolas, polyval(egyutthatok, xAbrazolas), "r")
plot(xAbrazolas, polyval(egyutthatok9, xAbrazolas), "k")
title("4. Pelda: meresi adatok + modellek")
legend("(x, y) meresi pontok", "legjobb egyenes modell", "legjobb 9-edfoku modell (Lagrange interpolacio)", "location", "northeastoutside")

###########################################################################

disp("\n\n\
=============================================\n\n\
5. pelda\n\n\
Tekintsuk az alabbi 300 db zajos mereset egy jelensegnek: x, y\n\n\
(nem iratjuk ki mert hosszu volna)")

meresekSzama = 300

x = linspace(-4, 4, meresekSzama);
y = 3 * x.^2 + 2 * x + 1 + randn(1, meresekSzama); # y = [3 x^2 + 2 x + 1] + zaj

disp("\n\
Tegyuk fel, hogy nem ismerjuk a jelenseg modelljet, nezzuk a meresi pontok\n\
grafikonjat, 12. abra")

figure (12)
plot(x, y, "b+")
title("5. Pelda: meresi pontok")
legend("(x, y) meresi pontok", "location", "northeastoutside")

disp("\n\
Egyertelmuen parabolat latunk! Tehat jo esellyel egy masodfoku polinom a \n\
jelenseg fizikai modellje amit megfigyelunk. Igy egy masodfoku polinomot\n\
illesztunk legkisebb negyzetek modszerevel a mintara:\n\n\
egyutthatok = polyfit(x, y, 2)")
egyutthatok = polyfit(x, y, 2)

disp("\n\
A 13. abran latjuk, hogy az illesztett modell hogyan nez ki.")
figure (13)
plot(x, y, "b+")
hold on 
plot(x, polyval(egyutthatok, x), "r")
title("5. Pelda: meresi pontok + modell")
legend("(x, y) meresi pontok", "illesztett modell", "location", "northeastoutside")

###########################################################################

disp("\n\n\
=============================================\n\n\
6. pelda\n\n\
Tekintsuk az alabbi 100 db zajos mereset egy jelensegnek: x, y\n\n\
(nem iratjuk ki mert hosszu volna)")

meresekSzama = 100

x = linspace(4, 6, meresekSzama);
y = exp(2 * x + 3) + randn(1, meresekSzama); # y = e^(2x + 3) + zaj

disp("\n\
Tegyuk fel, hogy nem ismerjuk a jelenseg modelljet, nezzuk a meresi pontok\n\
grafikonjat, 14. abra")

figure (14)
plot(x, y, "b+")
title("6. Pelda: meresi pontok")
legend("(x, y) meresi pontok", "location", "northeastoutside")

disp("\n\
Latszik, hogy az x tengelyen alig valtoznak az x ertekek, mig az y tengelyen\n\
hatalmasat valtozik a mert ertek. Ez tipikusan exponencialis jelensegre utal.\n\
Igy logikusnak tunik egy y ~ e^(a1 * x + a0) modell felallitasa, illetve\n\
ezen a1 es a0 egyutthatok keresese legkisebb negyzetek modszerevel.\n\n\
Ezt direktben is el lehet vegezni, de mi most csak a polinom illesztest\n\
tanultuk, igy kicsit trukkozni kell.\n\n\
Bar ez most ""trukkozes"", a valodi alkalmazasokban is ha \n\
exponencialis / logaritmikus tendenciat eszlelnek a meresekben, akkor\n\
szokas lineariazlni a merest -- atterni logaritmikus koordinata-rendszerre\n\
( log-log, lin-log, log-lin koordinata-rendszerek ), melyekben a mert adatok\n\
mar linearis tendenciat mutatnak.\n\n\
Lehet, hogy ez valami furcsa dolognak hangzik, pedig eleg gyakran\n\
talalkoztok vele. Logaritmikus skalan merjuk pl a hangerosseget (decibel, dB),\n\
a foldrengest (Richter).\n\n\
Pl a 40 dB az 10 000-szer olyan eros mint az 1 dB,\n\
a 100 dB az 10 000 000 000-szor olyan eros mint az 1 dB.")

disp("\n\n\
Tehat a 14. abra alapjan azt sejtjuk, hogy y ~ exp( a1 * x + a0 ),\n\
ebbol kovetkezoen log(y) ~ a1 * x + a0   (ahol log(.) a termeszetes alapu logaritmus)\n\n\
Ezt ellenorizzuk a 15. abran.")

figure (15)
plot(x, log(y), "b+")
title("6. Pelda: meresi pontok lin-log logaritmikus koordinata-rendszerben")
legend("(x, log(y)) meresi pontok", "location", "northeastoutside")

disp("\n\n\
Itt mar alkalmazhatjuk a tanult legkisebb negyzetes egyenes illesztest:\n\
egyutthatok = polyfit(x, log(y), 1)")
egyutthatok = polyfit(x, log(y), 1)

disp("\n\
16. abra az illesztett egyenes es a meresi pontok.")

figure (16)
plot(x, log(y), "b+")
hold on 
plot(x, polyval(egyutthatok, x), "r")
title("6. Pelda: meresi pontok lin-log logaritmikus koordinata-rendszerben + modell")
legend("(x, log(y)) meresi pontok", "illesztett modell", "location", "northeastoutside")


disp("\n\
17. abra az illesztett egyeneshez tartozo exponencialis fgv es a meresi pontok.")
figure (17)
plot(x, y, "b+")
hold on 
plot(x, exp(polyval(egyutthatok, x)), "r")
title("6. Pelda: meresi pontok  + modell")
legend("(x, y) meresi pontok", "illesztett modell", "location", "northeastoutside")

disp("\n\
Az illeszkedes gyakorlatilag tokeletes akarhanyszor is futtatjatok le\n\
a programot ennek az oka, hogy a jelhez kepest a zaj nagyon kicsi.\n\
Ennek ellenkezojerol majd latunk peldat az orai feladatokban.")
