Matematiska institutionen Chalmers Tekniska Högskola och Göteborgs Universitet

Statistisk databehandling ht 95

Demonstration 1 : S^{+}

Programpaketet Splus (/S^{+}) är en utvidgning av programpaketet S. För att använda S^{+} krävs en S^{+}-manual [delas ut (gratis) i stencilform som kursmaterial], men ofta också en S-manual (ett begränsat antal finns att låna hos institutionssekreteraren).

S^{+} startas från Unix med kommandot Splus

besov> Splus
S-PLUS : Copyright (c) 1988, 1993 Statistical Sciences, Inc.
S : Copyright AT&T.
Version 3.2 Release 1 for Sun SPARC, SunOS 5.x :  1993
WWorking data will be in /users/math/palbin/.Data
>

S^{+} har två sorters grundläggande kommandon; expressions som evalueras direkt

> sqrt(2)
[1] 1.414214
>

och assignements med operatorn <- eller ->

> x <- 3**(0.5)
> x
[1] 1.732051
>

Förutom med Return kan kommandon även separeras med ;

> x <- log(2); exp(x) -> y
> x; y
[1] 0.6931472

> [1] 2

Om ett kommando ej är fullständigt vid ett radslut erhålles prompten + (istf. >) och S^{+} väntar på att kommandot skall avslutas

> sin(pi/4) *
+ 2**(-1/2)
[1] 0.5

Komprimerad hjälp beträffande en function fås mha. args(function)

> args(log)
function(x, base = 2.71828182845905)
NULL
> log(3, 3)
[1] 1

(Här betyder base = 2.71828182845905 att om log-funktionens andra argument base ej anges så får det “default”-värdet 2.71828182845905.)

Mycket mera omfattande hjälp fås med help(function) eller ? function

> help(gamma)
Gamma Function (and its Natural Logarithm)
DESCRIPTION:
Returns the gamma function or the log of the gamma
function.
USAGE:
gamma(x)
lgamma(x)
REQUIRED ARGUMENTS:
x:  numeric or complex object.  Missing values (NAs) are
allowed.
VALUE:
gamma or natural log of gamma function evaluated for each
value in x.
CLASSES:
This function will be used as the default method for
classes that do not inherit a specifc method for the
function or for the Math group of functions.  The result

will retain the class and the attributes.  If this
behavior is not appropriate, the designer of the class
should provide a method for the function or for the Math
group.
DETAILS:
NAs are returned when evaluation would cause numerical
problems, except that when an overflow occurs, an infinity
is returned.
These are members of the Math group of generic functions.
WARNING:
gamma(x) increases very rapidly with x.  Use lgamma to
avoid overflow.
BACKGROUND:
For positive integral values, gamma(x) is (x - 1)!.  For
negative integral values, gamma(x) is undefined.
SEE ALSO:
exp, sum, cumsum.
EXAMPLES:
gamma(6) # same as 5 factorial

Då help anropats blir man “kvar i” help: Man kommer ut med ett q

q
> gamma(6)
[1] 120

S^{+} är ett vektorspråk och dess funktioner opererar på vektorer som kan skapas med c(...)

> c(pi/2, pi/4) -> x
> y <- c(1, 1/sqrt(2))
> sin(x); y; sin(x)+y; cos(x)*y
[1] 1.0000000 0.7071068
> [1] 1.0000000 0.7071068
> [1] 2.000000 1.414214
> [1] 6.123234e-17 5.000000e-01

Man kan också tolka om vektorer till matriser och S^{+} arbetar även med matriser (men komponentvis; se multiplikationen nedan)

> x <- c(1, 2, 3, 4, 5, 6)

> matrix(x, 2, 3) -> y; y**2 -> z
> y; z
    [,1] [,2] [,3]
[1,]   1    3    5
[2,]   2    4    6
>      [,1] [,2] [,3]
[1,]    1    9   25
[2,]    4   16   36

Den ordning i vilken elementen från vektorn x hamnar i matrisen y är ej den ordning många förväntar sig: Detta är en vanlig källa till programmeringsfel!

Det är lätt att referera till enstaka element i vektorer/matriser

> x[4]; y[2,2]; z[2,1]
[1] 4
> [1] 4
> [1] 4

Noga räknat har S^{+} två sorters vektorer: Först har vi “vanliga” sk. atomiska (eng. atomic) vektorer där varje element är ett tal (reellt eller komplext), en teckensträng eller ett logiskt värde (TRUE eller FALSE). Den andra sortens vektorer kallas rekursiva (eng. recursive) och har något eller några element som i sin tur är (atomiska eller rekursiva) vektorer eller matriser. Rekursiva vektorer benämnes också listor.

Listor kan skapas mha. kommandot list och man kan testa om en vektor är rekursiv (dvs. en lista) eller atomisk.

> vekt1 <- c(10, TRUE); vekt2 <- c(5, FALSE)
> lista <- list(vekt1,vekt2)
> is.atomic(vekt1); is.atomic(lista)
[1] T
> [1] F
> is.recursive(lista)
[1] T

Element i listor kan refereras till med dubbla “bracket”-paranteser

> lista[[1]]; lista[[1]][1]
[1] 10 1
> [1] 10

(Observera att det logiska värdet TRUE här skrivs ut som en etta.) Ofta ger man elementen i listor namn och då kan man referera till listans element mha dessa namn

> nylista <- list(elem1=vekt1,elem2=vekt2)
> vekt <- nylista$elem2
> vekt[1]; vekt[2]
[1] 5
> [1] 0

Observera att vissa funktioner i S^{+} (som tex. funktionen lsfit för linjär regression) har en lista som funktionsvärde. Element (dvs. vektorer eller matriser) hos dessa listor kan var mycket långa. Mha. $ kan man referera till de delar av funktionens värde-lista man är intresserad av (tex. nås regressionskoefficienterna från lsfit mha. lsfit$coef) och man slipper då utskrift av hela listan. Filer kan läsas in till S^{+} mha. funktionen scan; tex. kan filen data.s med utseendet

1 4
9 16 25

läsas in med

> x <- scan("data.s")
> sqrt(x)
[1] 1 2 3 4 5

Observera "..." runt filnamnet. För filer som ej finns i det directory S^{+} startades från måste man ange directory före filnamnet. Tex. kan filen data2.s med utseendet

1 4
9 16

belägen i directory /users/math/palbin/user läsas in med

> y <- scan("/users/math/palbin/user/data2.s")
> sqrt(y)
[1] 1 2 3 4

(Det fungerar tyvärr ej med scan("˜ /user/data2.s").)

Funktionen runif(n,min,max) ger n st. oberoende slumptal med likformig fördelning över intervallet [min, max].

> runif(5,1,2)
[1] 1.377098 1.496360 1.174278 1.553202 1.222433

Man kan definiera egna funktioner i S^{+}

> my.func <- function(x,y,z) sin(x)*(y**2)*exp(-z)
> my.func(pi/2,1,0)
[1] 1

Om U_{1}, …, U_{n} är likformigt fördelade slumptal över [0, 1] och F en kontinuerlig fördelningsfunktion med invers F^{−1}, så är X_{1} = F^{−1}(U_{1}), …, X_{n} = F^{−1}(U_{n}) slumptal med fördelning F.

Tex. har Weibull(a, b)-fördelningen F (x) = 1−exp{−(ax)^{b}} invers F^{−1}(y) = (− ln(1− y))^{1/b}/a och vi skapar en vektor X med 100 st. Weibull(1, 2)-slumptal med kommandona

> Finv <- function(a,b,x) (-log(1-x))**(1/b)/a
> X <- Finv(1,2,runif(100,0,1))
> X
 [1] 0.52001833 0.19623383 0.53592498 1.10367285 0.14149955 0.16924995
 [7] 1.02710354 1.24994595 1.21247929 1.18556703 0.90047396 0.75665895
[13] 0.93281487 0.29455009 1.46431833 0.24884127 0.25850963 1.05261643
[19] 0.92792423 0.08847906 1.39457337 0.46879738 0.50197726 1.01928398
[25] 1.39565971 1.29879476 1.11979418 0.65056252 1.36409325 0.41278285
[31] 0.27459245 1.66527698 0.36725148 0.26930565 0.66945274 1.36103108
[37] 0.87773169 1.09268324 0.55348753 0.45254122 1.33691553 1.07740933
[43] 1.68681872 1.97153569 0.60634804 1.54671789 0.85536057 1.40520623
[49] 0.88406436 0.98593807 0.76230990 0.36471573 0.71770006 0.73169387
[55] 1.52935740 0.28117906 2.56611372 1.13520558 1.03277283 0.75403254
[61] 0.90017817 2.20934149 0.75507928 0.86350916 0.70798122 1.09016944
[67] 0.06855055 0.34149005 0.33377808 1.37707546 1.10396063 0.88172748
[73] 0.69343386 0.38276787 1.03960206 0.22478093 0.81594809 0.67868702
[79] 1.16401369 1.19375474 1.35777449 0.72666033 0.34248738 0.82506722
[85] 1.73747984 0.64890786 0.58855173 0.24514754 1.09874622 0.29502881
[91] 0.94474390 1.02251647 0.76238478 0.93971260 1.84041605 1.62627373
[97] 1.51179879 0.19341062 0.30974725 0.23781206

Vissa fördelningar finns inprogrammerade i S^{+}. Tex. fås 25 st. N(5, 2)-slumptal med

> N <- rnorm(25, 5 ,2); N
> [1] 7.5020368 6.9515621 6.1662037 1.3967953 5.9963663 0.7638452 3.6892214
  [8] 6.2716771 4.9303204 5.7730262 5.4781295 1.4003307 4.4297668 4.1349703
 [15] 6.4611820 3.8329255 4.6703808 2.1754967 3.7131461 5.2736179 8.9416210
 [22] 6.2194380 7.5743943 3.2613323 5.3978069

[rweibull(n, shape, scale) ger n st. Weibull(a, b)-slumptal med (OBS!) a=1/shape och b=scale).]

Det finns ett flertal vanliga deskriptiva statistikor (strunta i varningstexten!)

> mean(N); var(N); median(N)
[1] 4.896224
Warning messages:
Looking for object "mean" of mode "function", ignored one of mode "numeric"
> [1] 4.226016
> [1] 5.273618
> min(N); quantile(N, 0.25); quantile(N, 0.5); quantile(N, 0.75); max(N)
[1] 0.7638452
> 253.713146
> 505.273618
> 756.219438
> [1] 8.941621

För tillgång till grafik i S^{+} måste funktionen openlook() anropas. Man får då ett grafiskt fönster (som kan ikonifieras eller behållas). Nu kan tex. vektorn x plottas mot (lika “långa”) vektorn y med plot(x, y) kommandot.

Den empiriska fördelningsfunktionen F_{n} för ett stickprov x_{1}, …, x_{n} definieras

= 0, −∞<x<x_{(1)},
F_{n}(x) = i/n, x_{(i)} ≤x<x_{(i+1)}, i ∈ {1, …, n−1},
= 1, x_{(n)} ≤x<∞.

(där x_{(1)}, …, x_{(n)} är det ordnade stickprovet). Man kan sådeles åskådliggöra den empiriska fördelningsfunktionen för stickprovet N skapat ovan med

> openlook(); plot(sort(N),1:25/25)

> plot(sort(N),1:25/25,type="l")

Grafisk utskrift flyttas till default skrivare mha. postscript(). Därefter kommer plottar att sparas i filen ps.out.0001.ps och (så småningom) skrivas ut. Grafiken flyttas tillbak till skärmen med dev.off(). Först då [eller då S^{+}-körningen avslutas] skrives ps.out.0001.ps ut på skrivaren. Nästa gång postscript() anropas sparas grafik i ps.out.0002.ps osv.

> hist(N)

> boxplot(N)

En S^{+}-körning avslutas med q().

Avslutningsvis nämner vi att S^{+} har en funktion qqplot(...) med vars hjälp man grafiskt kan studera hur väl ett stickprov x_{1}, …, x_{n} ansluter till en viss fördelning (-sfunktion) F:

Om y är en vektor som innehåller det sorterade stickprovet (så att y_{i} =x_{(i)}) och om z är en vektor vars element ges av z_{i} =F^{−1}((i−1/2)/n) (dvs. z_{i} är F ’s (i−1/2)/n-fraktil) så ger kommandot qqplot(y,z) en plot av y mot z där god anslutning till en rät linje indikerar att stickprovet ansluter väl till F. (Vektorn ((1−1/2)/n, …, (n−1/2)/n) kan skapas med kommandot ppoints(n,0.5).)

Fallet då F = Φ (normalfördelning) är så viktigt att det har det egna kommandot qqnorm(x). Här behöver man varken skapa en vektor y med det sorterade stickprovet eller en vektor z med Φ’s fraktiler.