[vastaus aiempaan viestiin]
| Kirjoittaja: | Seppo Mustonen |
|---|---|
| Sähköposti: | - |
| Päiväys: | 13.3.2005 11:23 |
Kilpatehtävän 5 (Epälineaarisen regressiomallin estimointi)
ratkaisivat Jari Lipsanen, Joonas Kauppinen ja Kalle Suominen.
Olen pyytänyt ja saanut luvan kommentoida heidän ratkaisujaan.
Tehtävähän oli tavallaan mahdoton, sillä olennaisesti monotonisesti
muuttuvia tilanteita pystyy mallintamaan hyvin erilaisin
funktiomuodoin.
Jotta voisi löytää sen "oikean", pitäisi tietää mahdollisimman
paljon itse ilmiöstä, josta havainnot on saatu. Harkitsin ennen
tehtävänantoa pitkään, olisiko ollut syytä kertoa jotain aineiston
(kuvitellusta) taustasta, mutta päädyin lopulta hyvin niukalle linjalle.
Koko homman idea juontaa alkunsa yli kymmenen vuoden takaa.
Tällöin Survosta oli käytössä kouluissa supistettu SURVOS-versio ja
meillä oli yhteistoimintaa useiden aktiivisten matematiikan opettajien
kanssa. Tässä touhussa silmiini osui pieni aineisto CUP, jonka oli
kerännyt Rolf Dahlström oppilaineen Pohjoispuiston lukiossa
Hyvinkäällä.
Aineisto CUP esiintyy mm. Survon englanninkielisessä käyttöoppaassa
(sivulla 256) ja piirrosesimerkkinä myös nykyisessä Survossa
(DEMO -> BOOK -> /EX-P256).
Toistan sen tässä:
Temperature T (in C) of a coffee cup as a function of time t (in min)
DATA CUP:(t,T) 0,69.5 1,65.0 2,62.5 3,60.0 4,58.0 5,56.0 6,54.5 7,
53.0 8,51.5 9,50.0 10,48.5 11,47.5 12,46.5 13,45.5 14,44.5 15,43.5 16,
43.0 17,42.0 18,41.0 19,40.5 20,39.5 21,39.0 22,38.5 23,38.0 24,37.5 25,
36.5 END
Kokeessa on siis seurattu kuuman kahvikupin jäähtymistä 25 minuutin
ajan. Jossain opettajille tarkoitetussa monisteessa oli yritetty
kehitellä tilastollista mallia sille, millä tavalla lämpötila T
(asteina) muuttuu ajan t (minuutteina) mukaan mutta mielestäni hyvin
vaatimattomin tuloksin.
En hallitse fysiikkaa, mutta amatöörinä arvelin, että lämmön haihtuminen
jotensakin muistuttaa säteilyilmiöiden intensiteetin laskua, jolloin
oli luonnollista kokeilla eksponentiaalista mallia ja käyttää
Survon ESTIMATE-operaatiota seuraavasti:
MODEL MCUP1
T=T0+a*exp(-b*t)
T0=24 a=1 b=0.1
ESTIMATE CUP,MCUP1,CUR+1 / METHOD=M RESULTS=0
Estimated parameters of model MCUP1:
T0=32.004 (0.65557)
a=36.078 (0.553967)
b=0.0774714 (0.00322618)
n=26 rss=4.540517 R^2=0.99788 nf=244
Malli on suhteellisen hyvästä selitysasteestaan huolimatta kelvoton,
sillä esim. rajalämpötila kupin täydelleen jäähdyttyä olisi T0=32
astetta eli siis liikaa; tuskin luokka tai laboratorio on on ollut
mikään huonosti lämmitetty sauna.
Päätin yleistää mallia sellaiselta ajatuspohjalta, että lämpöä katoaa
kupista kahdella tapaa toisistaan riippumatta eli erikseen nesteen
pinnalta höyryävänä ja erikseen kupin seinämien ja pohjan kautta.
Tällöin olisikin kyseessä kahden eksponentiaalisen komponentin
summa eri parametrein. Kokeilin mallia
MODEL MCUP2
T=T0+a*exp(-b*t)+c*exp(-d*t)
T0=24 a=1 b=0.1 c=1 d=0.01
ESTIMATE CUP,MCUP2,CUR+1 / METHOD=M RESULTS=0
Estimated parameters of model MCUP2:
T0=29.2965 (0.664362)
a=3.4019 (0.439222)
b=0.875557 (0.222983)
c=36.7547 (0.347845)
d=0.0632298 (0.00282832)
n=26 rss=0.641183 R^2=0.99970 nf=1341
Mallin selitysaste on parantunut selvästi, mutta raja-lämpötila on
edelleen ehkä korkea, T0=29 astetta.
Olisi todella arvokasta tietää, mikä oli ympäristön lämpötila kokeen
aikana. Oletetaan, että se olisi ollut 24 astetta. Tällöin vakioimalla
T0 saadaan (Huom. risuaita T0:n edessä alkuarvossa)
.....................
#T0=24 a=1 b=0.1 c=1 d=0.01
ESTIMATE CUP,MCUP2,CUR+1 / METHOD=M RESULTS=0
Estimated parameters of model MCUP2:
a=6.82352 (0.62342)
b=0.337007 (0.0499705)
c=38.347 (0.675188)
d=0.0445664 (0.000970874)
n=26 rss=1.108772 R^2=0.99948 nf=560
Malli ei T0:n vakioinnista juuri huonontunut. Selitysaste on korkea
ja karkeat diagnostiset tarkastelut (residuaalien jakauma) osoittavat,
että mallintaminen on onnistunut kohtalaisen hyvin.
Oikean T0-arvon puuttuminen vain aiheuttaa tarpeetonta epämääräisyyttä
parametrien arvoissa.
* * *
Tämän esimerkin innoittamana loin Kilpatehtävä 5:n aineiston
seuraavasti:
FILE CREATE KILPA5,16,4
FIELDS:
1 N 4 Y (#.#####)
2 N 4 X (###)
END
FILE INIT KILPA5,100
VAR X,Y TO KILPA5
X=ORDER
Y=round(Y1,5)
Y1=2*exp(-0.03*X)+4*exp(-0.15*X)+N.G(0,0.02^2,rand(2005))
Tämä 100 havainnon aineisto noudattaa siis MCUP2-tyyppistä mallia
parametrein T0=0, a=2, b=0.03, c=4, d=0.15 ja jäännöstermi on
normaalinen odotusarvolla 0 ja ja varianssilla 0.02^2=0.0004 eli
teoreettinen jäännösneliösumma on 100*0.02^2=0.04.
Kun piirtää kuvan aineistosta, on vaikea välttyä ajatukselta, että
Y lähestyy nollaa, kun X kasvaa ja tämän uskoin hieman auttavan
tehtävän ratkaisijoita. Esim. tässä mallityypissä voi silloin
rauhassa asettaa T0=0 ja yleensäkin riittää rajoittua malleihin,
joissa Y lähestyy nollaa X:n kasvaessa.
...........................
Tällä "oikealla" mallilla saadaan tulokset:
MODEL M2
Y=a*exp(-b*X)+c*exp(-d*X)
a=2 b=0.1 c=1 d=0.1
ESTIMATE KILPA5,M2,CUR+1 / METHOD=M RESULTS=0
Estimated parameters of model M2:
a=2.03758 (0.0349776)
b=0.030543 (0.000373821)
c=3.97454 (0.0314643)
d=0.151668 (0.00210342)
n=100 rss=0.034732 R^2=0.99970 nf=260
eli parametriestimaatit ovat hyvin lähellä todellisia ja jäännösneliö-
summa 0.0347 on jopa alle teoreettisen arvon 0.04. Tämä ei ole
outoa, sillä estimoinnissa jää aina varaa mistä valita :)
Esittelen nyt kilpailijoiden ratkaisut. Jokainen käytti ESTIMATEa
mutta aika erilaisin mallimäärityksin. Yhteistä kaikille on kuitenkin
parametrien lukumäärä 4 kpl. kuten "oikeassakin".
Olen alla toistanut estimoinnit omalla tavallani mutta tulokset
ovat samat kuin kilpailijoiden alkuperäiset.
......................................................................
Jari Lipsanen:
MODEL MODLX
Y=d+a/((X^b)+c)
d=0 a=1 b=1 c=0
ESTIMATE KILPA5,MODLX,CUR+1 / METHOD=M RESULTS=0
Estimated parameters of model MODLX:
d=-0.20477 (0.0114716)
a=60.9412 (1.6263)
b=1.13812 (0.0115319)
c=9.83792 (0.331222)
n=100 rss=0.046322 R^2=0.99960 nf=460
Jari pääsee parhaaseen selitysasteeseen, mutta negatiivisen raja-arvon
d kustannuksella. Malli heikkenee jonkin verran, jos d pakotetaan
nollaksi.
.....................................................................
Joonas Kauppinen:
MODEL MODKX
Y=a+b*LOG(X)+c*d^X
a=2.5 b=-0.5 c=3 d=1
ESTIMATE KILPA5,MODKX,CUR+1 / RESULTS=0
Estimated parameters of model MODKX:
a=2.51591 (0.0456916)
b=-0.531731 (0.0108803)
c=3.31122 (0.0427112)
d=0.896585 (0.00127065)
n=100 rss=0.072271 R^2=0.99938 nf=430
Mallin muoto on todella erikoinen mutta tulos hyvä aineistossa nähdyllä
vaihtelualueella. Tässä kuitenkin valitettavasti Y -> -oo (joskin
hyvin hitaasti) b*log(X)-komponentin vaikutuksesta.
.....................................................................
Kalle Suominen:
MODEL MODSX
Y=a*LOG(X)+b+c*X+g*SIN(X)
a=-1 b=6 c=0 g=0
ESTIMATE KILPA5,MODSX,CUR+1 / METHOD=M RESULTS=0
Estimated parameters of model MODSX:
a=-1.61288 (0.0201797)
b=5.88216 (0.047276)
c=0.0168984 (0.000644573)
g=-0.0155013 (0.0116917)
n=100 rss=0.655222 R^2=0.99435 nf=16
Tässäkin mallin muoto on varsin eksoottinen ja mallia vaivaavat
samanlaiset asymptoottiset ongelmat kuin Joonaksen.
.....................................................................
Laskin jokaisesta mallista residuaalit ja niiden autokorrelaatiot
20 ensimmäisestä havainnosta, joissa mahdolliset systemaattiset
poikkeamat ovat todennäköisimpiä.
Viive "oikea" Jari Joonas Kalle
1 -0.0565 0.4078 0.3133 0.4254
2 0.0002 0.4295 0.2763 0.2321
3 0.0278 0.2192 0.0213 0.0752
Tämä osoittanee, että jäännöstermien käyttäytyminen kilpailijoiden
malleilla on hiukan liian systemaattista jo nähdyllä
vaihtelualueellakin.
* * *
Kahden tai useamman eksponentiaalisen termin summan muodostama malli on
yleensä varsin hankala estimoitava. Koska sitä tarvitaan mm.
tietyissä fysikaalisissa ja kemiallisissa ilmiöissä kts. esim.
http://www.statsci.org/other/prony.html
on sen estimointiin kehitetty erikoismenetelmiä, jotka ovat parempia
kuin raaka pns-estimointi, jota ESTIMATE käyttää.
-Seppo
| Vastaukset: |
|---|
Survo-keskustelupalstan (2001-2013) viestit arkistoitiin aika ajoin sukrolla, joka automaattisesti rakensi viesteistä (yli 1600 kpl) HTML-muotoisen sivukokonaisuuden. Vuoden 2013 alusta Survo-keskustelua on jatkettu entistäkin aktiivisemmin osoitteessa forum.survo.fi. Tervetuloa mukaan!