6. lecke: Egy mérési feladat elemzése
| Cél: MATLAB-os tanulmányainkat egy konkrét gyakorlati mérnöki probléma elemzésével és megoldásával zárjuk. A felhasznált adatok valós mérésből származnak, a megválaszolandó kérdések szintén a gyakorlatból jönnek. A MATLAB által nyújtott támogatást kimondottan változatos módon használjuk fel a kidolgozás közben. Legfontosabb célunk itt egy olyan minta prezentálása, amely a későbbiekben Önnek mankóként szolgálhat valamely saját problémájának a megoldásához. Bár ez a (későbbi) probléma a most bemutatottól nyilván jelentős mértékben különbözni fog, a megoldás általános technikája (tervezés, modellkészítés, apparátus összegyűjtése, szemléltetés stb.) hasonló lesz. |
Kiegészítő anyagrészként - a törzsanyagon túlmenően - a mostani probléma finomabb, mélyebb részleteibe is betekintünk, szintén azzal a célzattal, hogy a későbbi önálló munkához ötleteket és mankót biztosítsunk. |
Követelmények: Ön akkor sajátította el megfelelően a tananyagot, ha (a MATLAB segítségével) |
- be tud tölteni a MATLAB-ba szöveges fájlból nagyobb méretű mátrixokat, vektorokat, és ezek megfelelő (adott) részeit ki tudja vágni;
- végre tud hajtani nagyobb méretű adathalmazon alapvető statisztikai vizsgálatokat;
- saját szavaival be tudja mutatni a paraméterbecslési feladat két alaptípusát;
- saját szavaival be tudja mutatni az interpolációs és a regressziós feladatot;
- végre tudja hajtani a Basic Fitting opció felhasználásával a lineáris és a polinomiális regressziót, továbbá létre tudja hozni konkrét értékek előrejelzésére a regressziós polinomot;
- végre tudja hajtani a simítási feladatot a csúszóátlag módszerrel és a szűrési feladatot adott szűrőfüggvénnyel.
|
Időszükséglet: A tananyag elsajátításához (a feladatok megoldásával együtt) hozzávetőlegesen 3 órára lesz szüksége. |
Kulcsfogalmak |
- Alapstatisztikák, hisztogram;
- Paraméterbecslés, deriváltbecslés;
- Interpoláció és regresszió;
- Regressziós polinom, együtthatók, hibanorma;
- Csúszóátlag, simítás, szűrés.
|
Bevezetés |
A mostani lecke legnagyobb részében egy olyan komplex gyakorlati feladatot oldunk meg, amely a MATLAB modul zárásaként összegzi tudásunkat. A kidolgozás során feltételezzük azt, hogy a lépéseket követő hallgató megfelelő módon ismeri és alkalmazni is tudja a MATLAB modul előző leckéiben megtanult anyagrészeket/ismereteket. Ez a lecke alapvetően már nem az új anyagrészek átadásáról szól, hanem a tanulási folyamat lezárásaként - a kitűzött gyakorlati probléma megoldása kapcsán - bemutatja azt, hogy mire lehetünk képesek a MATLAB segítségével, hogyan segíthet ez a programrendszer a mérnöki munkában. Természetesen, egy gyakorló műszaki szakember számára a tanulás soha nem ér véget: ezt azzal is érzékeltetjük, hogy néhány szükséges új ismeretet most fogunk bevezetni. |
Feladatunk összetettségéből következik, hogy ha az alapszintnél igényesebb, és ezáltal a valóságot jól közelítő modellt szeretnénk felépíteni (és a hozzá tartozó megoldást is be szeretnénk mutatni), akkor ahhoz néhány olyan haladóbb ötletre is szükség van, amelyek nem tartoznak bele a tantárgy törzsanyagába. Hangsúlyozzuk, hogy ezek fakultatív ismeretek (* jellel megjelölve) és nem kerülnek számonkérésre sem. |
A feladatban egy rögzített rugalmas test nyúlásának mérési adatait kell feldolgoznunk, a hibás adatok kiszűrésével, statisztikákkal, regresszióval. A megoldás része a megfelelő modell felállítása is (ellenőrzéssel). |
A leírt lépéseket követve Ön is oldja meg az útmutató alapján az egyes részfeladatokat. |
Adatoszlopok kinyerése, ábrázolás, alapvető elemzés |
Egy egyszerű szövegszerkesztővel nyissuk meg az adatfájlt (gyak-fel-adatok.txt ), és tanulmányozzuk a struktúráját. Látható, hogy öt adatoszlopunk van, amelyek rendre a megnyúlás (l) és az erő (F) értékeit tartalmazzák (1. és 2., ill. 3. és 4. oszlopok), továbbá egy azonosítót (5. oszlop), összesen 999 adatsorral. |
A mérési adatokat tartalmazó szövegfájl részlete | |
A fejléc számunkra most felesleges, ezért töröljük. Mentsük el a fájlt (data.txt). A továbbiakban külön említés nélkül az új szövegfájllal dolgozunk. |
Ebben a feladatban a 3. és a 4. oszlop adatait fogjuk felhasználni. Olvastassuk be a fájlt a MATLAB rendszerrel és töltsük be az adatokat egy l és egy F vektorba! |
>> load data.txt, l = data(:,3); F = data(:,4); |
Ellenőrizze a Workspace ablak segítségével, hogy az l és az F vektorok tartalma valóban a megfelelő! |
Próbálja meg az eredeti mérési adatokat tartalmazó fájlt (gyak-fel-adatok.txt) betölteni a MATLAB-ba! Értelmezze a rendszer üzenetét! |
Következő lépésként ábrázoljuk az erőt a megnyúlás függvényében! |
>> plot(l, F) |
"Nyers ábra" (elemzéshez) az l és F vektorokról | |
Az F(l) függvény grafikonját vizsgálva azt találjuk, hogy három részre bontható fel. |
Az első szakasz - feltehetően a befogás és a mérés technikai korlátai miatt - még nem ad reális képet a fizikai jelenségről ("inicializálás", kb. az első 70 adat). Ezután egy exponenciális/logaritmikus jellegű rész következik, majd egy lineáris szakasz figyelhető meg (a Hooke-törvény szerint). |
Minket a feladat "standard" megoldása során elsősorban a lineáris rész érdekel majd, a fejlettebb (fakultatív) megoldásnál elemezzük a nemlineáris részt. |
Osztáspontok vizsgálata, alapstatisztikák |
Idézze fel az alap Excel modulban tanultakat a statisztikai függvényekről! |
Az alapstatisztikák (egy konkrét adatsokaságról) természetesen a MATLAB-ban is nagyon egyszerűen, beszédes nevű parancsokkal kikérhetők. |
Néhány példa (egyszerű próbák): |
>> max(l), mean(l) % l legnagyobb eleme és elemeinek átlaga ans = 9.1125 ans = 4.4925 >> median(l) % l elemeinek mediánja (a rendezett minta közepe) ans = 4.4688 >> mean(l)==max(l)/2 % megnézzük, hogy az átlag egyenlő-e a max. felével ans = 0 |
Összefoglalóan, külön részletezés nélkül közöljük a MATLAB alap statisztikai függvényeit (megjegyezzük, hogy a Toolboxok ennél lényegesen többet tudnak - ezeket a lehetőségeket azonban jegyzetünkben nem tudjuk tárgyalni). |
A MATLAB alap statisztikai függvényei (egyváltozós eset)Függvény | Jelentés | max | Maximális elem | min | Minimális elem | mean | Átlag (várható érték torzítatlan becslése) | median | Rendezett minta közepe | std | Tapasztalati szórás | sum | Összeg | cumsum | Kumulatív részösszegek | sort | Növekvő sorrendbe rendezés | diff | Szomszédos elemek differenciái | hist | Hisztogram (gyakoriság oszlopdiagram) |
1. táblázat |
Ha valamelyik statisztikai függvény (fogalom) jelentésében bizonytalan, akkor lapozza fel a megfelelő matematikai tankönyvet (Matematika 3., valószínűség számítás és statisztika), illetve nézze meg a MATLAB súgót a példákkal! |
A MATLAB alap statisztikai függvényei (többváltozós eset)Corrcoef | Korrelációs mátrix | Cov | Kovariancia mátrix |
2. táblázat |
Határozza meg a megnyúlás vektor elemeinek összegét és szórását! |
Néhány statisztikai függvény felhasználásával - a fenti egyszerű próbákon túlmenően - végrehajtunk a mért adatokon egy konkrét vizsgálatot. |
Az adatokat áttekintve láthatjuk, hogy jól érzékelhető szabályosság figyelhető meg a megnyúlási értékek sorozatában. Vizsgáljuk meg a MATLAB-bal, hogy vajon egyenlőközűnek vagy közel egyenlőközűnek tekinthetők-e az osztáspontok! |
Ehhez elkészítjük a szomszédos l értékek differenciáit (távolságát). |
>> dl = diff(l); |
Hasonlítsa össze az l és dl lista elemeinek a számát! Mivel magyarázza az eltérést? |
A dl vektor elemeit manuálisan is átnézhetjük (Workspace ablak) - legalábbis egy darabig -, de ez a módszer az adatok nagy száma miatt nem elég hatékony, és ráadásul kisebb eltérések így nem is feltétlenül érzékelhetők. |
A dl különbségvektor részlete a variable editor ablakban | |
Ezért egy hisztogramot (gyakoriság oszlopdiagramot) kérünk az adatokra, hogy jobban áttekinthessük a struktúrát. |
>> hist(dl) |
Mind a kézi átnézés, mind a grafikon azt mutatja, hogy az adatsorban csak 0,0063 és 0,0125 körüli értékek szerepelnek. Ha egy módosított hisztogramot is kikérünk (az összes eset darabszámával - hist(dl, 999)), akkor erről már jól látjuk, hogy valójában csak pontosan két érték fordul elő (0,00625 és 0,0125). |
Hisztogramok az l értékek differenciáiról | |
Ha még a második hisztogram sem "győzött meg" bennünket teljes mértékben erről, akkor meg is számoltathatjuk az elméleti értékek nagyon kis sugarú környezetébe eső adatok számát a pontos ellenőrzéshez (ehhez ciklusszervező utasítást és feltételvizsgálatot is bevetünk; a 10*eps sugarú környezetre a kis adattárolási pontatlanságok miatt van szükség): |
>> db = 0; for i = 1:998 if (dl(i) > 0.0125 - 10*eps & ... dl(i) < 0.0125 + 10*eps) db = db + 1; end; end; db db = 460 |
Teljesen hasonlóan eljárva, a 0,00625 nagyon pici környezetébe 538 adat fog esni, megvan tehát mind a 998. |
A grafikon részekre bontása |
A következőkben megkeressük az F(l) diagram egyenesnek tűnő szakaszát. A 2. ábrát elemezve azt a munkahipotézist fogalmazzuk meg, hogy nagyjából a 250. esettől már közel lineáris a kapcsolat. Ezt a feltevést deriváltbecsléssel jól ellenőrizhetjük: ha tényleg lineáris a függvényünk, akkor ettől kezdve a differenciahányadossal becsült derivált közel konstans lesz. |
Eltároljuk az eredeti adatsort, és lecsípjük belőle a vizsgálni kívánt részt (most: 100-tól): |
>> l_orig = l; F_orig = F; l = l_orig(100:end); F = F_orig(100:end); |
Elkészítjük a különbségsorozatokat és a differenciahányadosok sorozatát: |
>> dF = diff(F); dl = diff(l); deriv_F = dF./dl; % vigyázzunk a pontozott műveletre! |
Most ábrázolni szeretnénk az egyes l értékekhez tartozó differenciahányadosokat. Itt azonban - első próbálkozásra - kellemetlen meglepetésben van részünk: |
>> plot(l, deriv_F) ??? Error using ==> plot Vectors must be the same lengths. |
A hibát az okozza, hogy a két vektor hossza nem azonos (a különbségsorozat eggyel rövidebb). Ezt úgy tudjuk kijavítani, hogy még egy példányban odamásoljuk a lista végére az utolsó elemet. |
>> deriv_F(end+1) = deriv_F(end); |
Így már elvégezhető az ábrázolás: |
>> plot(l, deriv_F) |
A differenciahányados-értékek ábrája | |
Az ábrán jól látszik, hogy feltevésünkkel ellentétben az F'(l) deriváltfüggvény kb. 4,2-ig csökkenő tendenciát mutat, azaz eddig a függvény biztosan nem lineáris. |
Idézze fel az Excel számítások modul 4. leckéjében tanultakat a paraméterbecslési feladatról! |
A deriváltra a nem lineáris szakaszon a következő közelítő becslést adjuk (indoklás: a függvény alakja alapján látjuk, hogy negatív kitevős exponenciális típusú, megfelelő eltolással): |
|
ahol d az egyenes rész meredeksége az F = F(l) kapcsolatban. Mivel egy exponenciális függvény határozatlan integrálja szintén exponenciális függvény, ezért az F(l) függvény jó közelítéssel |
|
alakú lehet. (Itt kihasználtuk azt a becslésnél szokásos technikát, hogy könnyebb a deriváltból következni az eredeti függvényre, mint direkt módon becsülni az eredetit.) |
A paraméterbecslési feladatnak két típusa létezik: az első esetben maga a matematikai modell ismert, de az aktuális paraméterek értékei ismeretlenek (ilyen típusú feladattal találkoztunk korábban, az Excel számítások modul 4. leckéjében); a második esetben nem ismerjük a matematikai modellt, ekkor ismert alapfüggvények, például algebrai polinomok, ill. trigonometrikus vagy exponenciális függvények kombinációit tekintjük a modell egy közelítésének (a másodlagos cél ekkor ezek együtthatóinak a meghatározása). |
Ahogy már az Excel számítások modul 4. leckéjében is említettük, egy függvény paraméteres alakjának helyes megbecslése - a fent bemutatott módon - általában csak komoly tudás és nagy gyakorlati tapasztalat birtokában végezhető el! |
Folytatva a feladatunk megoldását, a továbblépésre ezután két lehetőségünk van. |
a) Csak az egyenesnek látszó részen dolgozunk (ez a rész tartozik az Informatika II. tárgy törzsanyagába), azaz az adataink második felén, kb. az 500. ponttól kezdve. |
*b) Az exponenciális rész további vizsgálatához a (2) összefüggés A, b, c, d, e paramétereit a négyzetes eltérések minimalizálásával keressük meg. (A mért adatsor és a becsült függvényértékek eltérései négyzetének összegét minimalizáljuk az fminsearch() függvény segítségével.) Erre a feladatra visszatérünk a lecke legvégén. |
Regresszió a lineáris részre, előrejelzés |
Az interpoláció és a regresszió alaphelyzete az, hogy adottak a jelenséghez tartozó (összetartozó) bemeneti-kimeneti értékpárok. Ha ezek az értékpárok hibátlannak tekinthetők és számuk egyenlő a keresett paraméterek (együtthatók) számával, akkor interpolációs feladatról beszélünk. Ha ezek az értékpárok hibával terheltek és számuk meghaladja a meghatározandó paraméterek (együtthatók) számát, akkor regressziós feladattal állunk szemben. |
A technikai megoldás legfontosabb lépése az, hogy meghatározzuk a jelenséghez tartozó modellfüggvényt (közelítő függvényt). A modellfüggvény ezután alkalmas lesz arra, hogy a bemenet/kimenet kapcsolatát olyan pontokban is megadjuk, amelyekben korábban nem ismertük. |
Az interpoláció tipikus technikái: polinomiális interpoláció, spline interpoláció, két- vagy többváltozós interpoláció (ezekkel mi az Informatika II. tananyagban a terjedelmi korlátok miatt nem tudunk foglalkozni). A regresszió tipikus technikái: lineáris regresszió, nemlineáris regresszió. |
A mi esetünkben (mostani mintafeladat és az Informatika II. tárgy keretein belül bemutatott, illetve számon kért feladatok) ismert a konkrét megvalósítás kimért adatsora és még a műszaki-fizikai jelenség modellje is. Ez általában ténylegesen regressziós feladat. |
Ezen belül a fontosabb feladattípusok (és célok) a következők: |
- A modell paramétereinek minél jobb becslése, mert ennek műszaki-fizikai jelentése van (hővezetési tényező, rugalmassági modulus stb.);
- A mért pontok közötti intervallumokban is szeretnénk minél jobb becslést adni a lehetséges értékekre;
- A teljes mért intervallumon kívül is szeretnénk előrejelzést adni a további függvényértékekre (trend).
|
A közelítés jóságának mérőszámai: |
- A független változó mért értékeinél a függő változó mért és becsült értékeinek korrelációs együtthatója;
- Ugyanennek a négyzete az R2 (determinációs együttható);
- Az eltérések szórása.
|
A MATLAB az ilyen típusú feladatok megoldásához változatos megoldási lehetőségeket nyújt, az alap lehetőségektől (amelyek már szintén meglehetősen "erősek") a Toolboxok változatos parancskészletéig (pl. nlinfit, cftool). |
A legegyszerűbb esetben adott x, y pontsorozathoz kérhető polinomos közelítés a plot rajzablak menüjében. |
Feladatunkban mi most a lineáris részre kérünk regressziót. Ehhez töröljük a változóinkat, újra betöltjük az adatokat, majd vesszük az adathalmaz 2. felét, és kirajzoltatjuk. |
>> clear all, load data.txt >> l_orig = data(:,3); F_orig = data(:,4); >> l = l_orig(500:end); F = F_orig(500:end); plot(l, F) |
Lépjünk be a Figure 1 ablak Tools/Basic Fitting ("alapszintű" illesztés) menüpontjába. Az ablakban válasszuk ki, ill. jelöljük be a következő opciókat: |
- Lineáris kapcsolat;
- Show equations (mutassa a rendszer az ábrán az egyenleteket);
- Plot residuals (rajzolja ki a rendszer a hibatagokat);
- Show norm of residuals (jelenítse meg a rendszer a hibák normáját).
|
A "Significant digits" (szignifikáns jegyek) értékét növeljük meg 4-re. |
A gomb megnyomásával nyissuk ki a "Numerical results" panelt is. Itt lehetőség van az illeszkedés adatainak mentésére, ill. direkt módon is kimásolhatjuk a parancsablakba a megfelelő változókat. |
Lineráris regresszió a Basic Fitting eszközzel | |
Láthatjuk, hogy az illeszkedés alapvetően pontos, ugyanakkor egyes szakaszokon a kirajzolódó "hullámzó" hibák azt sejtetik, hogy magasabb fokú görbe még pontosabban leírhatja a mért adatok viselkedését. (Ennek elemzésével részletesen nem foglalkozunk.) |
A Basic Fitting eszközzel végezze el a regressziós közelítést másod- és negyedfokú polinomokkal is. Milyen hibanormák adódnak így? |
Az elkészült regressziós függvényt felhasználjuk arra, hogy a mért tartományon kívül is megjósoljuk a várható - fizikai - viselkedést. Ehhez másoljuk be a parancsablakba a megfelelő paramétereket: |
>> p1 = 0.20365; p2 = 0.41372; |
Ezután legyártjuk néhány megfelelő x értékkel az képlettel meghatározott y értékeket, amelyek így jó becslést adnak az F függvény értékére (most három új pontot hozunk létre). |
>> x = [10 11 12]; y = p1*x + p2; ertektabla = [x; y] ertektabla = 10.0000 11.0000 12.0000 2.4502 2.6539 2.8575 |
Végül egy új ábrára kirajzoljuk újra a mérési adatsor pontjait (teljes adathalmaz), és ugyanide elhelyezzük az újonnan számolt pontokat is: |
>> figure(2), plot(l_orig, F_orig), hold on, plot(x, y, 'ok') |
Előrejelzés az adatokra | |
A regressziós polinom alapján adjon két intervallumon belüli pontra is becslést, és ellenőrizze a becslés pontosságát a valódi értékek alapján (pl. x = 4 és x = 6 helyek). |
Az adatok simítása és szűrése |
A mért adatsorozat feldolgozásának fontos lépése a mérnöki gyakorlatban a simítás (egyenetlenségek csökkentése) és ennek részeként egyes feltűnően kiugró adatok törlése (ezek feltehetően mérési hibából származnak). |
Ha a Zoom In eszközzel kellő mértékben ránagyítunk a mérési adatainkat bemutató ábra megfelelő részére, akkor tanulmányozhatjuk a finom struktúrát: egy ismétlődő apró "hullámzás" figyelhető meg a grafikonon. A Pan eszközzel "mászkálva" találhatunk olyan helyeket, ahol a lokális eltérések nagyobbak. |
Az adatok finomszerkezete | |
Rátérve a simítási feladat megoldására, a MATLAB lehetőségeit itt ketté kell választanunk. Ha elérjük a Signal Processing Toolbox szolgáltatásait, akkor a probléma kezelésére nagyon változatos eszköztár áll rendelkezésünkre (a Toolboxok apparátusát jegyzetünkben nem tárgyaljuk). |
Ha csak a MATLAB alapszolgáltatásai használhatóak, akkor egyszerű simításként alkalmazhatjuk például az 5-pontos csúszóátlag módszert, amely öt egymás melletti érték átlagával (bal oldali két szomszéd, jobb oldali két szomszéd és maga az alappont) helyettesíti a középsőt. Erre készítettünk egy saját függvényt (shape5.m ). |
Saját simító függvény (5 pontos csúszóátlag) | |
A megvalósításnál egy MATLAB-os for ciklust használtunk, a vizsgálatból értelemszerűen ki kell hagyni az adatsor első kettő és utolsó kettő elemét. |
A feladatunk megoldása a függvényt használva: |
>> l = l_orig; F = F_orig; >> Z = shape5(F); |
Kirajzoltatjuk az eredeti és a simított adatsorozat eltérését: |
>> plot(l, F-Z) |
Az eredeti és a simított adatsorozat eltérése | |
Elemezze az ábrát nagyobb nagyítással is. Alap statisztikai függvényekkel írassa ki a legnagyobb eltérést! A fent megismert módszerrel (osztáspontok vizsgálata) számoltassa meg, hogy hány esetben nagyobb az eltérés 0,0004-nél, 0,00025-nél. |
Láthatjuk, hogy az eltérések zöme - abszolút értékben - 0,00025-nél kisebb. |
Ha a mért adatsor és a simított adatsor közötti különbség jelentős (ill. jelentősnek ítéljük), akkor a kiugró eseteket érdemes törölni (szűrési feladat), mint extra zajból, vagy mérési hibából származó értékeket. (Mi ezt a határt most 0,0003-nél húzzuk meg.) |
Erre a feladatra szintén egy saját függvényt készítettünk (outlier_filter), amely paraméterként megkapja az eredeti adatsor (x, y) mellett a simított értékek sorozatát is (ys), és az x(k), y(k), ys(k) adathármasokból csak azokat hagyja meg, ahol az eltérés a szintén megadott limitnél kisebb. Az eredmények rendre az xf, yf, ysf sorozatokba kerülnek, amelyeket üres sorozatokból kiindulva készítünk el. (Fontos megjegyezni, hogy az új sorozatok hossza általában eltér az eredetitől!) |
Saját szűrő függvény | |
A függvény segítségével most tehát lecsípjük a 3E-4-nél nagyobb hibájú adatokat: |
>> limit = 3E-4; >> [ls, Fs, Zs] = outlier_filter(l, F, Z, limit); |
Kirajzoltatjuk az eredeti és a simított-szűrt adatsorozat eltérését: |
>> plot(ls, Fs-Zs) |
Az eredeti és a simított-szűrt adatsorozat eltérése | |
Megnézzük, hogy az eljárás hány adatot törölt. |
>> length(ls) ans = 939 |
Látjuk, hogy nagyjából az adatok 6%-át "vesztettük el". |
Ezután ábrázoljuk a redukált adatsorozat második felét, és az ábrára a megszokott módon lineáris regressziót kérünk, együtthatókkal, hibanormával (Basic Fitting; 6.13. ábra). |
>> l = ls(470:end); F = Fs(470:end); plot(l, F) |
Az együtthatók új értékei 0,20332 és 0,4162, míg a régiek 0,20365 és 0,41372; hasonlóan az új hibanorma 0,07968, a régi pedig 0,09042. |
Ezeket az eredményeket úgy összegezhetjük, hogy jól láthatóan a simítás és a pár kiugró pont törlése kevéssé módosította a korábbi értékeket, hiszen a sorozat eléggé lineáris volt. (A Signal Processing Toolbox parancsait használva is nagyon hasonló eredmény adódik.) |
Lineráris regresszió a simított-szűrt adatsorozatra | |
A fent megismert módon (megfelelő paraméterek a parancsablakba, y értékek legyártása stb.) adjon előrebecslést most is az F(l) függvény értékeire a 10, 11 és 12 helyen. Hasonlítsa össze az előrejelzett új értékeket a régiekkel! Mekkora eltérés tapasztalható? |
Készítsen új ábrát, amelyen bemutatja a szűrt adathalmaz mellett a most kapott három előrejelzett pontot! |
Az eredeti adatok felhasználásával oldja meg újra elölről a teljes simítási és szűrési feladatot oly módon, hogy 3 pontos csúszóátlag módszert használ (ehhez a shape5 függvényt át kell írni, új neve legyen shape3), majd a szűrés a fentivel egyezően a 0,0003-es limittel történjen. Hány pontot távolított el most a szűrés? |
Ezután ábrázolja a redukált adatsorozat második felét, és az ábrára a megszokott módon kérjen lineáris regressziót, együtthatókkal, hibanormával (Basic Fitting). Hasonlítsa össze az a kapott új értékeket a korábbiakkal! Mekkora eltérés tapasztalható? |
*Általános regresszió, paraméterbecslés a legkisebb négyzetek módszerével |
Ebben a - hangsúlyozottan fakultatív (!) - leckerészben egy olyan általános regressziós eljárást mutatunk be, amely a mostani konkrét feladat mellett sok más probléma megoldására is használható. Csak annyi a teendőnk, hogy a mintamegoldás megfelelő sorait átírjuk (ezek most piros kommentezésűek). |
Célunk a feladatnak megfelelően az exponenciális rész regresszióval történő becslése. |
Az eredeti adatsor 100. pontjánál kezdjük az illesztést. |
>> clear all >> load data.txt >> l = data(:,3); F = data(:,4); xdata=l(100:end); ydata=F(100:end); % adat-előkészítés |
Következik az típusú illesztés végrehajtása, a mellékelt (lásd lent) fitc függvény hívja meg az fminsearch parancsot. |
A fitc függvény részlete (a kódot bemutatjuk lent) | |
>> [estimates, model] = fitc(xdata,ydata) % a kapott paraméterek estimates = -0.3234 0.5801 1.1298 0.1902 0.5311 model = @fitc/expfunc >> est_fun = [num2str(estimates(1),4), '*exp(-', ... num2str(estimates(2),4), '*(x-', num2str(estimates(3),4), '))+ ', ... num2str(estimates(4),4), '*x+', num2str(estimates(5),4)] % a képletünk karakterlánca szövegdarabkákból összeragasztva est_fun = -0.3234*exp(-0.5801*(x-1.13))+ 0.1902*x+0.5311 >> [sse, FittedCurve] = model(estimates); % FittedCurve a becsléssel kapott értéksorozat >> sse % minimalizált négyzetes eltérés sse = 0.0091 >> R_square = min(min(corrcoef(ydata,FittedCurve).^2)) % determinációs együttható R_square = 1.0000 |
Ezután kirajzoljuk az adatokat és a regressziós függvényt, majd előrejelzést is készítünk. |
>> plot(xdata, ydata, '.'), hold on >> plot(xdata, FittedCurve, 'r') % az első rajz szándékosan vastag, a második (piros) vékony >> xlabel('x'), ylabel('f(estimates,x,forecast)') % x és y tengely címkéje >> title(['Fitting to function ', est_fun]); % diagramcím >> x = [9:0.2:10]; y = inline(est_fun), plot(x, y(x), 'k') % előrejelzés a [9, 10] tartományban % a vonal színe fekete y = Inline function: y(x) = -0.3234*exp(-0.5801*(x-1.13))+ 0.1902*x+0.5311 >> legend('data', est_fun, 'forecast', 2) % jelmagyarázat >> text(6,0.5, [{['R^2: ', num2str(R_square,4)]}; {['SSE: ', num2str(sse,4)]}]) % R-négyzet és SSE érték megjelenítése az ábrán (ellenőrzött, "semleges" helyen) % R-négyzet pozícionálását kell átírni megfelelő helyre >> hold off |
A parancssorozat kiadása után kapott eredmény (függvényrajz, előrejelzéssel, jelmagyarázattal) a 6.15. ábrán látható. |
Általános regressziós feladat megoldása | |
Az eredményt elemezve látható, hogy az illesztésünk R-négyzet értéke (determinációs együtthatója) gyakorlatilag 1, azaz az illesztés szinte tökéletes. Így a lineáris tag most kapott együtthatója az anyag fizikai viselkedését leíró rugalmassági együttható: 0.1902. |
Megjegyezzük, hogy ez pontosabb a lineáris regresszióval kapott értéknél! |
Végül közöljük a fitc() függvény definícióját is. Itt látunk példát függvénydefiníciók egymásba ágyazására és tetszőleges darabszámú paraméter átvételére (params). |
function [estimates, model] = fitc(xdata, ydata) % Call fminsearch with starting point. start_point = [-1, 0.5, 1, 0.2, 0.4]; % saját alkalmas kezdőpontok A, b, c, d, e-re model = @expfunc; % itt tesszük hozzáférhetővé a modellfüggvényt estimates = fminsearch(model, start_point); % expfunc accepts curve parameters as inputs, and outputs sse, % the sum of squares error for % A*exp(-b*(xdata-c) + d*xdata + e % and the FittedCurve. FMINSEARCH only needs sse, but we want % to plot the FittedCurve at the end. function [sse, FittedCurve] = expfunc(params) A = params(1); % 1. paraméter b = params(2); % 2. paraméter c = params(3); % 3. paraméter d = params(4); % 4. paraméter e = params(5); % 5. paraméter FittedCurve = A.*exp(-(b * xdata-c))+d*xdata+e; % a függvény definíciója ErrorVector = FittedCurve - ydata; sse = sum(ErrorVector .^ 2); end end |