Brz početak rada sa R: Eksponencijalni modeli (5. Deo)

U 3. i 4. delu koristili smo komandu lm() da ocenimo regresione modele. Videli smo kako da proverimo nelinearnost u našim podacima prilagođavanjem polinomskih modela i proverom da li se ti modeli bolje prilagođavaju podacima no linearni model. Sada ćemo da vidimo kako prilagoditi eksponencijalni model u R.

Kao i ranije, koristićemo datoteku sa brojevima (vreme raspada koji se odvija unutar izvora zračenja), merenih Gajgerovim brojačem u jednoj nuklearnoj elektrani. Tačke su registrovane tokom perioda od 30 sekundi za kratkotrajno veštačko radioaktivno jedinjenje. Učitavamo podatke i oduzmemo pozadinski broj od 623.4 brojanja u sekundi kako bismo dobili brojke koje se odnose na radioaktivni izvor. Kao i ranije, u našem R radnom prostoru smo isekli i nalepili sledeće podatke.

gajger <- structure(list(Vreme = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30), Brojevi = c(750, 725.2, 695, 708.5, 725, 690.9, 691.5, 678.6, 686.3, 668.9, 665.3, 669.7, 657.5, 661.6, 665.1, 648.1, 664.9, 647.9, 660, 643, 646.2, 653, 646.9, 639.3, 638.7, 636.8, 650.2, 633.2, 642.2, 649.3, 642.7)), .Names = c("Vreme", "Brojevi"), row.names = c(NA, -31L), class = "data.frame")
gajger

Opet uzimamo u obzir pozadinski broj i oduzimamo ga od skupa podataka.

pozadinskibroj <- 623.4
gajger$Brojevi <- gajger$Brojevi - pozadinskibroj

Pridružićemo celu datoteku kako bismo se mogli pozivati na sve promenljive direktno po imenu.

attach(gajger)
names(gajger)

Budući da fizika atomskog raspada uključuje eksponencijalne veze, da vidimo da li eksponencijalna funkcija odgovara podacima čak i bolje od kvadratne.

Eksponencijalna funkcija se može tretirati kao model promenljive log(Brojevi) po promenljivoj Vreme.

eksponencijalni.model <- lm(log(Brojevi)~ Vreme)
summary(eksponencijalni.model)

R daje sledeći rezultat:

Ovaj model je prilično dobar, iako objašnjava oko 83% varijacija u poređenju sa 90% objašnjenih kvadratnim modelom. Prikažimo grafički podatke i ocenjeni eksponencijalni model na grafikonu sa horizontalnom osom na kojoj su vrednosti vremena od 0 do 30 sekundi u intervalima od 0.1 sekunde.

vrednostivremena <- seq(0, 30, 0.1)
Brojevi.eksponencijalni <- exp(predict(eksponencijalni.model, list(Vreme= vrednostivremena)))
plot(Vreme, Brojevi,pch=16)
lines(vrednostivremena, Brojevi.eksponencijalni, lwd=2, col = "blue", xlab = "Vreme (s)", ylab = "Brojevi po Sekundi")

Imajte na umu da smo koristili eksponencijalnu vednost predviđenih vrednosti.

Dakle, prilagodili smo naš eksponencijalni model. Iako je fizika atomskog raspada eksponencijalna, u ovom slučaju ocenjeni eksponencijalni model odgovara podacima u manjoj meri od kvadratnog modela.

To je sve za sada. Sledeći put ćemo pogledati neku osnovnu sintaksu vezanu za grafikone.

David

Dodatak: Korišćeni R kodovi

# Učitavanje gajger datoteke u R radni prostor
gajger <- structure(list(Vreme = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30), Brojevi = c(750, 725.2, 695, 708.5, 725, 690.9, 691.5, 678.6, 686.3, 668.9, 665.3, 669.7, 657.5, 661.6, 665.1, 648.1, 664.9, 647.9, 660, 643, 646.2, 653, 646.9, 639.3, 638.7, 636.8, 650.2, 633.2, 642.2, 649.3, 642.7)), .Names = c("Vreme", "Brojevi"), row.names = c(NA, -31L), class = "data.frame")
gajger

# Dobijanje brojeva koji se odnose na izvor radioaktivnosti
# tako što će se oduzeti pozadinski broj
pozadinskibroj <- 623.4 
gajger$Brojevi <- gajger$Brojevi - pozadinskibroj

# Pridruživanje datoteke i pozivanje promenljivih po imenu
attach(gajger)
names(gajger) 

# Ocenjivanje eksponencijalnog regresionog modela
eksponencijalni.model <- lm(log(Brojevi)~ Vreme)
summary(eksponencijalni.model)

# Dijagram rasturanja i eksponencijalni regresioni model 
vrednostivremena <- seq(0, 30, 0.1)
Brojevi.eksponencijalni <- exp(predict(eksponencijalni.model, list(Vreme= vrednostivremena)))
plot(Vreme, Brojevi,pch=16)
lines(vrednostivremena, Brojevi.eksponencijalni, lwd=2, col = "blue", xlab = "Vreme (s)", ylab = "Brojevi po Sekundi")