KURZUS: Számítási módszerek

MODUL: Matematikai számítások MATLAB-bal

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
1. ábra

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
2. ábra

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ényJelentés
maxMaximális elem
minMinimális elem
meanÁtlag (várható érték torzítatlan becslése)
medianRendezett minta közepe
stdTapasztalati szórás
sumÖsszeg
cumsumKumulatív részösszegek
sortNövekvő sorrendbe rendezés
diffSzomszédos elemek differenciái
histHisztogram (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)
CorrcoefKorrelációs mátrix
CovKovariancia 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
3. ábra

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
4. ábra

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
5. ábra

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):

F " ( l )=a e b( lc ) +d

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

F( l )=A e b( lc ) +dl+e

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
6. ábra

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 y=p1x+p2  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
7. ábra

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
8. ábra

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)
9. ábra

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
10. ábra

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
11. ábra

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
12. ábra

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
13. ábra

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 A e b( lc ) +dl+e 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)
14. ábra
>> [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
15. ábra

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
Önellenőrző kérdések

1. Ebben a feladatban egy rögzített rugalmas test nyúlásának mérési adatait dolgozzuk fel.

Törölje a mellékelt gyak-fel-adatok.txt adatfájlból a fejléc sorokat, mentse el így a fájlt, majd töltse be a MATLAB-ba.

Olvassa be az 1. és a 2. oszlop adatait egy l (megnyúlás) és egy F (erő) vektorba!

Törölje az l és F adatvektorok első 450 elemét!

A továbbiakban csak az új megnyúlás és erő adatokkal dolgozzon!

Adja meg az új l vektor méretét (sor×oszlop) a Workspace ablakból leolvasva!

Válasz.:

A választ 4 tizedesjegy pontosan adja meg.
2. Mennyi az F értékek átlaga?

Válasz:

A választ 2 tizedesjegy pontosan adja meg.
3. Mennyi az F értékek összege?

Válasz:

4. Ábrázolja az erőt a megnyúlás függvényében!

A Basic Fitting opció segítségével illesszen másodfokú regressziós görbét a mérési adatokhoz! Jelenítse meg az ábrán az egyenletet 4 szignifikáns jeggyel és kérje ki az eltérési normát is!

Mennyi a másodfokú tag együtthatója?

Válasz:

5. A megfelelő panel segítségével másolja be a MATLAB-ba a regressziós polinom együtthatóit és ezzel számoljon!

A választ 4 tizedesjegy pontosan adja meg.
A regressziós polinom alapján hány N/mekkora becsült erő tartozna a 10 mm-es megnyúláshoz?

Válasz:

6. A megadott simítófüggvény (shape5.m ) felhasználásával (és esetleges módosításával) végezzen simítást az 5-pontos csúszóátlag módszerrel az erő adatokon!

A megadott szűrőfüggvény (outlier_filter.m ) felhasználásával (és esetleges módosításával) törölje azokat az adatokat a megfelelő adatvektorokból (l, F és simított F), amelyekre az eredeti és a simított erő értékek eltérése nagyobb, mint 3,2×10-4 N!

Hány adatpont törlődött az l vektorból az eljárás során?

Válasz:

Mekkora az új, leszűrt és megrövidült l vektor hossza?

Válasz: