iris
, , Setosa
Sepal L. Sepal W. Petal L. Petal W.
[1,] 5.1 3.5 1.4 0.2
[2,] 4.9 3.0 1.4 0.2
[3,] 4.7 3.2 1.3 0.2
[4,] 4.6 3.1 1.5 0.2
.........
, , Versicolor
Sepal L. Sepal W. Petal L. Petal W.
[1,] 7.0 3.2 4.7 1.4
[2,] 6.4 3.2 4.5 1.5
[3,] 6.9 3.1 4.9 1.5
[4,] 5.5 2.3 4.0 1.3
.........
, , Virginica
Sepal L. Sepal W. Petal L. Petal W.
[1,] 6.3 3.3 6.0 2.5
[2,] 5.8 2.7 5.1 1.9
[3,] 7.1 3.0 5.9 2.1
[4,] 6.3 2.9 5.6 1.8
........
Crabs
> crabs sp sex index FL RW CL CW BD 1 B M 1 8.1 6.7 16.1 19.0 7.0 2 B M 2 8.8 7.7 18.1 20.8 7.4 3 B M 3 9.2 7.8 19.0 22.4 7.7 4 B M 4 9.6 7.9 20.1 23.1 8.2 5 B M 5 9.8 8.0 20.3 23.0 8.2Painters:
> painters
Composition Drawing Colour Expression School
Da Udine 10 8 16 3 A
Da Vinci 15 16 4 14 A
Del Piombo 8 13 16 7 A
Del Sarto 12 16 9 8 A
...............
The Carraci 15 17 13 13 E
...............
Durer 8 10 10 8 F
Holbein 9 10 16 13 F
Pourbus 4 15 6 6 F
Rubens 18 13 17 17 G
Teniers 15 12 13 6 G
Van Dyck 15 10 17 13 G
Linear Discriminant Analysis
lda<-function(x, grouping, prior, tol = 0.0001)
{
n <- dim(x)[1]
p <- dim(x)[2]
g <- as.factor(grouping)
counts <- tapply(rep(1, length(g)), g, sum)
> counts
c s v
25 25 25
prop <- counts/sum(counts) # allow for supplied prior
if(missing(prior))
prior <- prop
# means by group (rows) and variable (columns)
group.means <- tapply(x, list(rep(g, ncol(x)), col(x)), mean)
> group.means
1 2 3 4
c 5.904 2.784 4.160 1.316
s 4.996 3.440 1.448 0.208
v 6.424 2.988 5.452 2.064
f1 <- sqrt(diag(var(x - group.means[g, ])))
> f1
[1] 0.4943573 0.3179325 0.3894695 0.1906461
if(any(f1 < tol))
stop(paste(" variable(s)", paste(format((1:p)[f1 < tol]),
collapse = " "), "appear to be constant within groups"))
scaling <- diag(1/f1)
> scaling
[,1] [,2] [,3] [,4]
[1,] 2.022828 0.000000 0.000000 0.00000
[2,] 0.000000 3.145322 0.000000 0.00000
[3,] 0.000000 0.000000 2.567595 0.00000
[4,] 0.000000 0.000000 0.000000 5.24532
X <- sqrt(1/(n - length(prop))) * (x - group.means[g, ]) %*% scaling
> var(X)
[,1] [,2] [,3] [,4]
[1,] 0.013888889 0.007908634 0.010613375 0.005338647
[2,] 0.007908634 0.013888889 0.006883927 0.007176478
[3,] 0.010613375 0.006883927 0.013888889 0.007414400
[4,] 0.005338647 0.007176478 0.007414400 0.013888889
> 1/0.013888889
[1] 72
X.s <- svd(X, nu = 0)
> X.s
$d:
[1] 1.6476488 0.8258493 0.7190378 0.4442096
$v:
[,1] [,2] [,3] [,4]
[1,] -0.5235503 -0.5369555 0.04565989 -0.6599159
[2,] -0.4853020 0.2540247 0.80326968 0.2339044
[3,] -0.5366184 -0.3234145 -0.41417823 0.6602274
[4,] -0.4499087 0.7365822 -0.42559298 -0.2718446
rank <- sum(X.s$d > tol)
if(rank < p)
warning("variables are collinear")
scaling <- scaling %*% X.s$v[, 1:rank] %*% diag(1/X.s$d[1:rank])
> scaling
[,1] [,2] [,3] [,4]
[1,] -0.6427658 -1.3152144 0.1284524 -3.005105
[2,] -0.9264297 0.9674758 3.5137812 1.656211
[3,] -0.8362332 -1.0055073 -1.4789794 3.816209
[4,] -1.4322925 4.6783469 -3.1046651 -3.209999
xbar <- apply(prior %*% group.means, 2, sum)
> xbar
1 2 3 4
5.774667 3.070667 3.686667 1.196
X <- sqrt((n*prior)/(length(prop) - 1)) * scale(group.means, xbar, F) %*% scaling
> var(X)
[,1] [,2] [,3] [,4]
[1,] 133.65578 -46.29705 273.5546 -89.07820
[2,] -46.29705 18.36213 -88.7734 29.44335
[3,] 273.55465 -88.77340 575.2826 -185.95171
[4,] -89.07820 29.44335 -185.9517 60.22636
X.s <- svd(X, nu = 0)
rank <- sum(X.s$d > tol * X.s$d[1])
scaling <- scaling %*% X.s$v[, 1:rank]
structure(list(prior = prior, means = group.means, scaling = scaling,
lev = levels(g), svd = X.s$d[1:rank]), class = "lda")}}
Example of results from lda:
tr <- sample(1:50,25)
train <- rbind(iris[tr,,1],iris[tr,,2],iris[tr,,3])
test <- rbind(iris[-tr,,1],iris[-tr,,2],iris[-tr,,3])
cl <- factor(c(rep("s",25),rep("c",25), rep("v",25)))
z <- lda(train, cl)
z
$prior:
c s v
0.3333333 0.3333333 0.3333333
$means:
1 2 3 4
c 5.904 2.784 4.160 1.316
s 4.996 3.440 1.448 0.208
v 6.424 2.988 5.452 2.064
$scaling:
[,1] [,2]
[1,] 0.8592838 -0.6290820
[2,] 2.0422890 2.5686477
[3,] -2.5358040 -0.8402268
[4,] -2.9940143 3.1021325
$lev:
[1] "c" "s" "v"
$svd:
[1] 39.509837 3.745213
> predict(z,test)$class
[1] s s s s s s s s s s s s s s s s s s s s s s s s s
c c c c c c c c c c c v c c c c c c v c c c c c c
v v v v v v v v v v v v v v v v v v c v v v v v v
.......................................
.
.
.
.
.
.
.
.
10..
.
s .
e .
c .
o .
n .
d .
.3 3 3 c
l . vvvv c2 c2c2c2
i . vv vv33s s1 cccccc
n 0.. 333vvvv1 s1ss 1 ccc22
e .v v3sv ss1s1 1 ccc
a v v v vss sssssss c
r . v 1 1
. s
d .
i .
s .
c .
r .
. .
-10..
v .
a .
r .
i .
a .
b .
l .
e .
.......................................
-5 0 5 10
first linear discriminant var.
> a_lda(painters[,1:4],School) > a$svd [1] 3.1837508 1.9269049 1.6496730 0.7660973 > confusion(School,predict(a,painters[,1:4])$class) A B C D E F G H A 5 0 1 2 0 0 0 2 B 4 1 1 0 0 0 0 0 C 0 2 2 0 2 0 0 0 D 0 0 0 9 0 0 1 0 E 0 0 0 1 4 0 1 1 F 1 0 0 0 0 2 1 0 G 1 0 0 1 1 0 4 0 H 0 0 1 0 0 0 0 3 attr(, "error"): [1] 0.4444444 > confusion(School,predict(a,painters[,1:4],dimen=2)$class) A B D E F G H A 6 0 2 0 0 0 2 B 6 0 0 0 0 0 0 C 3 1 1 1 0 0 0 D 0 0 8 0 0 2 0 E 2 0 1 2 0 1 1 F 1 0 0 0 1 1 1 G 1 0 1 1 0 4 0 H 0 0 0 0 1 0 3 attr(, "error"): [1] 0.5555556
Smoothing
> plot(times,accel,main="Polynomial Regression")
> lines(times,fitted(lm(accel~poly(times,3))))
> lines(times,fitted(lm(accel~poly(times,6))),lty=3)
> legend(40,-100,c("degree 3","degree 6"),lty=c(1,3),bty="n")
> plot(times,accel,main="Natural Splines",type="n")
> lines(times,fitted(lm(accel~ns(times,df=50))),lty=2)
> text(locator(),'df=50')
> plot(times,accel,main="Lowess")
> lines(lowess(times,accel))
> text(locator(),'default span')