# Alle eksempler fra `Noter om R' version 6.0.
# (Dog ikke eksempler indeholdt i opgaver eller Appendiks G.)

# Der er indsat kommentarer med afsnits- og sidenumre for at man
# lettere kan finde et givent eksempel.

# Afsnit 1.2. Side 7.
2 + 2

# Side 8.
5 -
3

# Afsnit 1.3.
plot(sin, 0, 6*pi)

# Afsnit 1.4. Side 9.
cbrt <- function(x) { x^(1/3) }

# Afsnit 2. Side 12.
17^(1/3)

# Afsnit 2.1.
options(digits=12)
17^(1/3)

options(digits=7)

# Afsnit 3. Side 14.
z <- 17^(1/3)

z

z^6 - z^3

z <- 5
z
z <- z + 3
z

ls()

rm(z)
z

# Afsnit 4. Side 15.
cbrt <- function(x) { x^(1/3) }

cbrt(1000)

root <- function(x, n) { x^(1/n) }

root(625, 4)

logis <- function(t) { K / (1 + a*exp(-r*t)) }

logis(5)

a <- 50
r <- 1
K <- 1
logis(5)
K <- 10
logis(5)

# Afsnit 5. Side 16.
plot(sin, 0, 6*pi)

p <- function(x) { 10 - 5*x^2 +  x^3 }
plot(p, -2, 6)

# Afsnit 5.1.
q <- function(x) { 7*x - x^2 }
plot(q, -2, 6, add=TRUE)

# Afsnit 5.2.
plot(q, -2, 6)
plot(p, -2, 6, add=TRUE)

# Side 18.
plot(q, -2, 6, ylim=c(-16, 100))
plot(p, -2, 6, add=TRUE)

# Afsnit 5.3.
plot(p, -2, 6, ylab="y")
plot(q, -2, 6, add=TRUE)

plot(p, -2, 6, ylab="y", main="To funktioner",
     sub="Polynomierne p(x) og q(x)")
plot(q, -2, 6, add=TRUE)

# Side 19.
plot(p, -2, 6, ylab="")
plot(q, -2, 6, add=TRUE)
title(main="To funktioner", sub="Polynomierne p(x) og q(x)")
title(ylab="y")

# Afsnit 5.4.
plot(p, -2, 6, col="blue")
plot(q, -2, 6, add=TRUE, col="red")

# Afsnit 5.5.
g1 <- function(x) { x }
g2 <- function(x) { 2 * x }
g3 <- function(x) { 3 * x }
plot(g1, xlab="x", ylab="g1 g2 g3", main="Linjetyper",
     lty="dashed", las=0)
plot(g2, add=TRUE, lty="dotted")
plot(g3, add=TRUE, lty="twodash")

plot(g1, xlab="x", ylab="g1 g2 g3", main="Linjetykkelser",
     lwd=1, las=1)
plot(g2, add=TRUE, lwd=2)
plot(g3, add=TRUE, lwd=3)

# Afsnit 5.6. Side 20.
plot(sin, 0,10, ylab="")
plot(cos, 0,10, add=TRUE, lty="dashed")
text(9.8, 0.2, "sin(x)")
text(8.3, -0.9, "cos(x)")

# Side 22.
usr <- par("usr")
text(usr[2], usr[4], "Teksten", adj=c(1.05, 1.25))

# Afsnit 5.7.
plot(sin, 0, 10, ylab="")
plot(cos, 0, 10, add=TRUE, lty="dashed")
legend(2.6, 1, c("sin(x)","cos(x)"), lty=c("solid","dashed"))

legend("bottomleft", c("sin(x)","cos(x)"), lty=c("solid","dashed"))

# Afsnit 5.8. Side 23.
locator(1)

# Afsnit 5.9. Side 24.
f <- function(x) { x * sin(1/x) }
plot(f, 0.001, 0.1)
plot(f, 0.001, 0.1, n=5000)

# Afsnit 5.10. Side 25.
plot(p, -2, 6, axes=FALSE)
axis(1, pos=0)
axis(2, pos=0)

plot(p, -2, 6, axes=FALSE)
axis(1, at=seq(-2,6), pos=0)
axis(2, pos=0)

plot(p, -2, 6, axes=FALSE)
axis(1, at=seq(-2,6), pos=0)
axis(2, pos=0)
axis(2, pos=0, at=seq(-20,50), labels=FALSE, tcl=-0.25)

# Side 27.
plot(p, -2, 6, axes=FALSE)
axis(1, at=seq(-2,6))
axis(2)
axis(2, at=seq(-20,50), labels=FALSE, tcl=-0.25)
box()

# Afsnit 6.2. Side 28.
p <- function(x) { 10 - 5*x^2 + x^3 }
win.metafile("poly.emf")
plot(p, -2, 6)
dev.off()

postscript("poly.eps")
plot(p, -2, 6)
dev.off()
pdf("poly.pdf")
plot(p, -2, 6)
dev.off()

png("poly.png")
plot(p, -2, 6)
dev.off()

# Afsnit 7. Side 29.
f <- function(x) { x - cos(x) }
res <- uniroot(f, c(-10, 10))
res

res$root

uniroot(p, c(-2,0)) $ root
uniroot(p, c(0,2)) $ root
uniroot(p, c(4,6)) $ root

# Side 30.
polyroot(c(10, 0, -5, 1))

Re( polyroot(c(10, 0, -5, 1)) )

# Afsnit 8. Side 31.
u <- function(t) { 20 * t / (t+1) + 10 }
f <- function(t) { 4000 * u(t) - 5000 * t }
optimize(f, interval=c(0,100), maximum=TRUE)

# Afsnit 9. Side 32.
integrate(sin, 0, pi)

f <- function(x) { 7*x^2+x^4 }
integrate(f, 0, 3)

res <- integrate(f, 0, 3)
res$value * 3

f1 <- function(x) { exp(-2*x) }
integrate(f1, 0, Inf)
f2 <- function(x) { 1/x^2 }
integrate(f2, 1, Inf)
f3 <- function(x) { exp(-x^2) }
integrate(f3, -Inf, Inf)

f4 <- function(x) { sin(1/x) }
integrate(f4, 0, 1)
integrate(f4, 0, 1, subdivisions=10000)

# Afsnit 10.1.1. Side 34.
d <- read.table("Grise2Fast.txt", header=TRUE, dec=",")

# Side 35.
d

fejl <- read.table("Grise2Fast.txt", header=TRUE)
fejl

fejl2 <- read.table("Grise2Fast.txt", dec=",")
fejl2

# Afsnit 10.1.2. Side 36.
d <- read.csv2("Grise2KommaDK.csv")

d <- read.table("Grise2KommaDK.csv", sep=";", dec=",", header=TRUE)

d <- read.csv("Grise2KommaUS.csv")

d <- read.table("Grise2KommaUS.csv", sep=",", header=TRUE)

# Afsnit 10.2. Side 37.
d <- read.table("clipboard", header=TRUE, dec=",")

d <- read.table(pipe("pbpaste"), header=TRUE, dec=",")

d <- read.table("X11_clipboard", header=TRUE, dec=",")

d <- read.table("clipboard", header=TRUE)

# Afsnit 10.3.
library(RODBC)
channel <- odbcConnectExcel("Grise2.xls")
d <- sqlFetch(channel, "Ark1")
odbcClose(channel)

# Afsnit 10.5. Side 38.
write.table(d, "minedata.txt")

d <- read.table("minedata.txt")

write.table(d, "minedata.csv", sep=";", dec=",", quote=FALSE,
            row.names=FALSE)

write.csv2(d, "minedata.csv", row.names=FALSE)

write.table(d, "output.csv", sep=",", quote=FALSE, row.names=FALSE)

write.csv(d, "output.csv", row.names=FALSE)

# Afsnit 10.6. Side 39.
summary(d)

summary(fejl)

nrow(d)

# Side 40.
head(d)
tail(d, 3)

# Afsnit 10.7.
d$Kontrol

d$Kontrol[16]

min(d$Kontrol)
max(d$Kontrol)

log(d$Kontrol)

# Side 41.
d$LogK <- log(d$Kontrol)
d$LogV <- log(d$Vaekst)
d

K <- d$Kontrol
K
K[16]
min(K)

attach(d)
Kontrol

detach(d)
Kontrol

Kontrol <- 17
attach(d)
Kontrol

# Afsnit 11. Side 42.
mean(d$Kontrol)

mean(d)

max(d$Kontrol, d$Vaekst)

range(d$Kontrol, d$Vaekst)

# Afsnit 12. Side 44.
plot(d)

# Afsnit 12.1.
plot(d$Tid, d$Kontrol, xlab="min", ylab="pH")

# Side 45.
plot(d$Tid, d$Kontrol, ylim=c(5,7), xlab="min", ylab="pH")
points(d$Tid, d$Vaekst, pch=4)

plot(d$Tid, d$Kontrol, ylim=range(d$Kontrol, d$Vaekst),
     xlab="min", ylab="pH")
points(d$Tid, d$Vaekst, pch=4)

plot(d$Tid, d$Kontrol, xlim=c(0,500), ylim=range(d$Kontrol, d$Vaekst),
     xlab="min", ylab="pH")
points(d$Tid, d$Vaekst, pch=4)

# Afsnit 12.3. Side 49.
plot(d$Tid, d$Kontrol, ylim=c(5,7), xlab="min", ylab="pH")
points(d$Tid, d$Vaekst, pch=4)
legend("topright", c("Kontrol","Vaekst"), pch=c(1,4))

legend(x, y, c("Linje","Symbol","Begge"),
       lty=c("solid","blank","solid"), pch=c(-1,1,1))

# Afsnit 12.4. Side 50.
t <- seq(0, 2*pi, len=100)

x <- cos(t)
y <- sin(t)
plot(x, y, type="l", asp=1, xlim=c(-2,2), ylim=c(-2,2))

points(0, 0, pch=3, col="red")

lines(x * t/3, y * t/3, col="blue")

# Side 51.
dx <- function(t) { cos(t)/3 - t*sin(t)/3 }
t0 <- uniroot(dx, c(1,5)) $ root
t0

x0 <- cos(t0)*t0/3
x0
lines(c(x0, x0), c(-2, 2), col="green")

# Afsnit 13. Side 52.
t <- c(1.0, 1.5, 2.0, 2.5, 3.0)
y <- c(1.4, 0.3, -1.5, -3.1, -4.8)
lreg <- lm(y ~ t)
lreg$coefficients

plot(t, y, ylim=c(-6,6), xlim=c(0,4), axes=FALSE)
axis(1, pos=0)
axis(2, pos=0)
abline(lreg, col="blue")

# Side 53.
model <- lm(Kontrol ~ Tid, data=d)
model$coefficients

plot(d$Tid, d$Kontrol, ylim=c(5, 7))
abline(model)

# Afsnit 13.1.
d$logKontrol <- log(d$Kontrol)
d$logTid <- log(d$Tid)

llmodel <- lm(logKontrol ~ logTid, data=d)
llmodel$coefficients

plot(d$logTid, d$logKontrol)
abline(llmodel)

# Side 54.
exp(llmodel$coefficients[1])

a <- llmodel$coefficients[2]
b <- exp(llmodel$coefficients[1])
pH <- function(t) { b*t^a }

plot(d$Tid, d$Kontrol, ylim=c(5, 7))
plot(pH, 0, 1500, add=TRUE)

# Afsnit 14.1. Side 55.
v <- c(3, 4)

v

c(v, v, 42)

seq(11,15)

seq(15,11)

11:15
15:11

seq(0, 1, by=0.1)
seq(20, 0, by=-4)

s <- seq(0, 2*pi, len=10)
s

# Side 56.
seq(0, by=0.1, len=14)

rep(1, 10)

rep(c(3, 4), 5)

rep(c(3, 4), each=5)
rep(c(3, 4), 5, each=2)

rep(c(3, 4), c(5, 2))

c(3, 8:2, 4, rep(5, 10))

# Afsnit 14.2.
u <- c(7, 9, 13, 107, 109, 113)
u
u[3]

u[c(3,4,4,1)]

u[1:4]

# Side 57.
u[-2]
u[c(-3,-5)]
u[-2:-4]

w <- 1:10
w
w[3] <- 17
w
w[c(2,6,8)] <- 101:103
w

# Afsnit 14.3.
sum(v)
mean(v)
min(v)
max(v)

summary(v)

c(1,2,3) + c(7,9,13)
c(1,2,3) * c(7,9,13)

sum(c(1,2,3) * c(7,9,13))

# Side 58.
sqrt(sum(v * v))

length(v)

c(3, 4) + c(11, 13, 17, 19)

c(3, 4, 3, 4) + c(11, 13, 17, 19)

c(2, 3, 5, 7, 11) * 2

c(2, 3, 5, 7, 11) * c(2, 2, 2, 2, 2)

sqrt(c(49, 81, 169))

cbrt <- function(x) { x^(1/3) }
cbrt(c(8, 27, 64))

# Side 59.
sqrt(c(1, 4, 9))
sapply(c(1, 4, 9), sqrt)
sum(c(1, 4, 9))
sapply(c(1, 4, 9), sum)

# Afsnit 14.4.
c("red", "green", "blue")
rep(c("red", "green", "blue"), 2)

# Afsnit 14.5.
v <- c(a=1,b=7)
v

2*v
c(v,5)

names(v)

names(v) <- c("x1","x2")
v
names(v)[1] <- "A"
v

# Side 60.
v["x2"]

names(v) <- NULL
v

# Afsnit 15. Side 61.
y <- 1
for (x in 2:6) { y <- y*x }
y
x

y <- 1
x <- 2
y <- y*x
x <- 3
y <- y*x
x <- 4
y <- y*x
x <- 5
y <- y*x
x <- 6
y <- y*x
y
x

# Afsnit 15.1. Side 62.
y <- 1
for (x in 2:6) { y <- y*x }

print(19)

y <- 1
for (x in 2:6) { y <- y*x ; print(x) ; print(y) }

# Side 63.
y <- 1
for (x in 2:6) { y <- y*x ; print(c(x,y)) }

# Afsnit 15.2.
x <- 0
for ( t in 1:100 ) { x <- x + runif(1, -1, 1) }
x

X <- rep(0, 100)
x <- 0
for ( t in 1:100 ) { x <- x + runif(1, -1, 1) ; X[t] <- x }

# Side 64.
X <- rep(0, 100)
x <- 0
x <- x + runif(1, -1, 1) ; X[1] <- x
x <- x + runif(1, -1, 1) ; X[2] <- x
 ... (i alt 100 linjer af denne type) ...
x <- x + runif(1, -1, 1) ; X[99] <- x
x <- x + runif(1, -1, 1) ; X[100] <- x

X

plot(0:100, c(0,X), type="l", xlab="t", ylab="x")

# Afsnit 15.3. Side 65.
fib <- rep(1,30)
for(i in 3:30) { fib[i] <- fib[i-1] + fib[i-2] }
fib

# Afsnit 16.1. Side 66.
A <- matrix(c(1,2,3,4), 2)

A

A2 <- matrix(c(1,2,3,4), 2, byrow=TRUE)
A2

M1 <- matrix(1:12, 3)
M1
M2 <- matrix(1:8, 4)
M2

# Side 67.
E <- diag(2)
E

t(M2)

# Afsnit 16.2.
A[1,2]
A[,1]
A[1,]

M1[c(3,1),]

M1[,c(2,4)]

M1[1, 2:4]
M1[1:2, 3]
M1[1:2, 3:4]

# Side 68.
M1[,-2]
M1[-1,]
M1[2, -3]

# Afsnit 16.3.
A <- matrix(1:4, 2)
v <- c(5,6)
w <- c(7,8)
cbind(v, w)
cbind(A, w)

rbind(v, w)
rbind(A,w)

c(v, w)

# Afsnit 16.4. Side 69.
A <- matrix(1:4, 2)
A
B <- matrix(9:6, 2)
B

A + B

A %*% B

# Side 70.
7 * A

M1 %*% M2

A %*% c(5,6)
c(5,6) %*% A
c(1,2,3) %*% c(7,9,13)

# Afsnit 16.4.1.
R <- A
R <- R %*% A
R <- R %*% A
R <- R %*% A
R

# Side 71.
R <- A
for (n in 2:4) { R <- R %*% A }
R

matpow(A, 4)

# Afsnit 16.5.
det(A)

solve(A)

A %*% solve(A)
solve(A) %*% A

A <- matrix(c(1,2,3,4), 2)
x <- solve(A, c(3,5))
x

# Side 72.
A %*% x

# Afsnit 16.6.
evA <- eigen(A)
evA
ev1 <- evA$vectors[,1]
ev2 <- evA$vectors[,2]

A %*% ev1
evA$values[1] * ev1

# Side 73.
M <-matrix(c(2.0,0.7,0,1.6,0,0.6,1.3,0,0.3),3)
M
evM <- eigen(M)
evM$values
evM$vectors

abs(evM$values)

# Afsnit 16.7.
M<-matrix(1:12,3)
M
nrow(M)
ncol(M)

dim(M)

# Afsnit 16.8. Side 74.
rownames(A)=c("r1","r2")
A
rownames(A)[2]="rk2"
A
rownames(A)

colnames(A)=c("s1","s2")
A

A["r1","s2"]

colnames(A) <- NULL
rownames(A) <- NULL
A

a<-1:2
b<-3:4
cbind(a,b)
rbind(a,b)

names(a)<-c("x1","x2")
cbind(a,b)

# Afsnit 16.9. Side 75.
M <- matrix(c(2.0, 0.7, 1.5, 0.4), 2)
M

v0 <- c(100,0)
v <- v0

for(t in 1:10) { v <- M %*% v }

v

M <- matrix(c(2.0, 0.7, 1.5, 0.4), 2)
v0 <- c(100,0)
v <- v0
V <- matrix(v0,2)

# Side 76.
for(t in 1:10) { v <- M %*% v ; V <- cbind(V,v) }

v
V

plot(0:10, V[1,]/1000, type="o", ylim=c(0, max(V)/1000),
     col="red", xlab="år", ylab="antal kaniner (tusinder)",
     main="Udvikling af kaninpopulation")
points(0:10, V[2,]/1000, type="o", pch=2, col="blue")
legend("topleft", c("Unge","Gamle"), pch=c(1,2), lty="solid",
       col=c("red", "blue"))

# Afsnit 17.1.1. Side 77.
y <- x <- 2

x <- 2

(x <- 2)

(plot(sin,0,6*pi))

# Afsnit 17.2. Side 78.
1+2; 2^5; sqrt(2)

1+2;

# Afsnit 17.2.1.
(1+2
-3)

1+2-
3

# Afsnit 17.3.
{ 1+2; 2^5; sqrt(2) }

# Side 79.
{ x<-2; y<-3; x+y }
x

{ x<-2
  y<-3
  x+y }

# Afsnit 17.4.
plot(function(x) { x^3 - 5*x^2 }, -2, 6)

f <- function(x) x^3 - 5*x^2
f <- function(x) ( x^3 - 5*x^2 )
f <- function(x) { x^3 - 5*x^2 }

# Afsnit 17.5.
1:40   # giver mig tallene fra 1 til 40
# Denne linje er kun en kommentar.

# Afsnit 18. Side 80.
M <- matrix(c(2.0, 0.7, 1.5, 0.4), 2)
v0 <- c(100,0)
v <- v0
V <- matrix(v0,2)
for(t in 1:20) { v <- M %*% v ; V <- cbind(V,v) }

fremskriv <- function(M, v0, antal.trin) {
   v <- v0
   V <- matrix(v0,2)
   for(t in 1:antal.trin) { v <- M %*% v ; V <- cbind(V,v) }
   V
}

# Side 81.
M <- matrix(c(2.0, 0.7, 1.5, 0.4), 2)
v0 <- c(100,0)
V1 <- fremskriv(M, v0, 20)

M <- matrix(c(2.0, 0.7, 1.5, 0.4), 2)
V1 <- fremskriv(M, c(100,0), 10)
V2 <- fremskriv(M, c(50,50), 10)
V3 <- fremskriv(M, c(0,100), 10)

plot.udvikling <- function(V) {
   årnumre <- 0:(ncol(V)-1)
   plot(årnumre, V[1,]/1000, type="o", col="red",
        main="Udvikling af kaninpopulation", xlab="år",
        ylab="Antal kaniner (tusinder)")
   points(årnumre, V[2,]/1000, type="o", col="blue")
   legend("topleft", c("Unge", "Gamle"),
          lty="solid", pch=1, col=c("red", "blue"))
}

plot.udvikling(V1)
plot.udvikling(V2)
plot.udvikling(V3)

# Side 82.
pct.unge <- function(V) { 100*V[1,]/(V[1,]+V[2,]) }

plot(0:10, pct.unge(V1), type="o", col="red", ylim=c(0,100),
     main="Udvikling af kaninpopulationer", xlab="år",
     ylab="Andel af unge i procent")
points(0:10, pct.unge(V2), type="o", col="blue")
points(0:10, pct.unge(V3), type="o", col="green")

kaningraf <- function(M, antal.trin) {
   V1 <- fremskriv(M, c(100,0), antal.trin)
   V2 <- fremskriv(M, c(50,50), antal.trin)
   V3 <- fremskriv(M, c(0,100), antal.trin)
   plot(0:antal.trin, pct.unge(V1), type="o", col="red", ylim=c(0,100),
        main="Udvikling af kaninpopulationer", xlab="år",
        ylab="Andel af unge i procent")
   points(0:antal.trin, pct.unge(V2), type="o", col="blue")
   points(0:antal.trin, pct.unge(V3), type="o", col="green")
}

M1 <- matrix(c(2.0, 0.7, 1.5, 0.4), 2)
kaningraf(M1, 10)
M2 <- matrix(c(4.0, 0.7, 2.5, 0.4), 2)
kaningraf(M2, 10)
M3 <- matrix(c(2.0, 0.9, 1.5, 0.6), 2)
kaningraf(M3, 10)

# Afsnit 19.1. Side 83.
x <- 3
x > 5
x < 5

x <- 1:10
x > 5
x == 3

# Side 84.
x[x>6]

sum(x>6)

x > 2 & x < 9
!(x > 2 & x < 9)
x <= 3 | x >= 7
x[!(x>2 & x<9) | x==5]

x > 2 && x < 9
x <= 3 || x >= 7

# Afsnit 19.1.1. Side 85.
heltal <- function(x) { floor(x) == x }

> heltal(42)
[1] TRUE
> heltal(4.2)
[1] FALSE

# Afsnit 19.2.
x <- 5
if (x<2) 2 else x
if (x>2) 2 else x

# Side 86.
if (x<2) 2
if (x>2) 2

if (x<2) 2
else x
if (x<2) 2 else
x

if (x<2) {
  2
} else {
  x
}

{ if (x<2) 2
  else x }

# Afsnit 19.2.1.
f <- function(x) { if (x<1) 0 else (x-1)/3 }

g <- function(x) { if (x<1) 0 else if (x<4) (x-1)/3 else 1 }

g <- function(x) {
   if (x<1) 0
   else if (x<4) (x-1)/3
   else 1
}

# Side 87.
g <- function(x) {
   if (x<0) 0
   else {
      if (x<4) (x-1)/3
      else 1
   }
}

# Afsnit 19.2.2.
primtal <- function(n) {
   er.primtal <- TRUE
   max.m <- floor(sqrt(n))
   for(m in 2:max.m) {
      if ( heltal(n/m) ) er.primtal <- FALSE
   }
   er.primtal
}

# Side 88.
primtal(21)
primtal(23)

primtal <- function(n) {
   er.primtal <- TRUE
   if (n > 3) {
      max.m <- floor(sqrt(n))
      for(m in 2:max.m) {
         if ( heltal(n/m) ) er.primtal <- FALSE
      }
   }
   er.primtal
}

# Afsnit 19.2.3. Side 89.
fibonacci <- function(n) {  # Fibonacci-tal for n>0.
  max.n <- max(n)           # Hvor mange fibonacci-tal skal bruges?
  fib <- rep(1, max.n)      # Lav en vektor med plads til max.n tal.
  if ( max.n > 2 ) {
    for ( i in 3:max.n ) {
       fib[i] <- fib[i-1] + fib[i-2]   # Udregn Fibonacci-tal nummer i.
    }
  }
  fib[n]                    # Returner de(t) n'te Fibonacci-tal.
}

fibonacci(6)
fibonacci(1:10)

# Afsnit 19.3.
f <- function(x) { if (x<1) 0 else (x-1)/3 }

g <- function(x) {
   if (x<1) 0
   else if (x<4) (x-1)/3
   else 1
}

# Side 90.
x <- 1:10
x
10-x
ifelse(x>5, x, 10-x)

f <- function(x) { ifelse(x<1, 0, (x-1)/3) }
g <- function(x) { ifelse(x<1, 0, ifelse(x<4, (x-1)/3, 1)) }

f(0:2)
g(0:5)

integrate(g,0,5)
plot(g, 0, 5)

# Afsnit 19.4.1. Side 91.
fak <- function(x) { if (x==1) 1 else x*fak(x-1) }
fak(4)

# Side 92.
fak <- function(x) {
   print(x)
   if (x==1) {
      print("if-gren valgt")
      1
   } else {
      print("else-gren valgt")
      x*fak(x-1)
   }
}

fak(4)

fak <- function(x) { ifelse(x==1, 1, x*fak(x-1)) }

# Side 93.
fak <- function(x) { ifelse(x<=1, 1, x*fak(x-1)) }

fak(1:5)

# Afsnit 19.5.
n <- 1
while(log(n)<5) { n <- n+1 }
n

# Side 94.
log(148)
log(149)

s <- 0 ; n <- 0
while(s<1000) { n <- n+1; s <- s+n }
n

# Side 95.
sum(1:44)
sum(1:45)

# Afsnit 19.5.1.
s <- 0 ; n <- 0
while(s<10) { n <- n+1; s <- s+n }

s <- 0 ; n <- 0
while(s<10) { n <- n+1; s <- s+n ; print(c(n,s)) }

# Side 96.
s <- 0 ; n <- 0
while(s<10) { n <- n+1; s <- s+n ; print(c(n,s)) ; flush.console() }

# Afsnit 19.5.2.
M <- matrix(c(2.0, 0.7, 1.5, 0.4), 2)
v0 <- c(2,0)
v <- v0

antal <- 0
while(sum(v)<1000) { antal <- antal+1 ; v <- M%*%v  }

antal
v
sum(v)

# Afsnit 19.5.3. Side 97.
primtal <- function(n) {
   er.primtal <- TRUE
   max.m <- floor(sqrt(n))
   for(m in 2:max.m) {
      if ( heltal(n/m) ) er.primtal <- FALSE
   }
   er.primtal
}

primtal <- function(n) {
   er.primtal <- TRUE
   max.m <- floor(sqrt(n))
   m <- 2
   while(m <= max.m) {
      if ( heltal(n/m) ) er.primtal <- FALSE
      m <- m + 1
   }
   er.primtal
}

primtal <- function(n) {
   er.primtal <- TRUE
   max.m <- floor(sqrt(n))
   m <- 2
   while(m <= max.m && er.primtal) {
      if ( heltal(n/m) ) er.primtal <- FALSE
      m <- m + 1
   }
   er.primtal
}

# Side 98.
primtal <- function(n) {
   max.m <- floor(sqrt(n))
   m <- 2
   while(m <= max.m && ! heltal(n/m) ) { m <- m + 1 }
   m > max.m
}

# Afsnit 19.5.4.
gæt.et.tal <- function() {
   tal <- floor(runif(1, 1, 101))
   antal.gæt <- 0
   gæt <- 0
   while(gæt!=tal) {
      indtastning <- readline("Indtast et gæt mellem 1 og 100: ")
      gæt <- as.numeric(indtastning)
      if(gæt < tal) { cat("Nej, tallet er større end", gæt, ".\n") }
      if(gæt > tal) { cat("Nej, tallet er mindre end", gæt, ".\n") }
      antal.gæt <- antal.gæt+1
   }
   cat("Du brugte", antal.gæt, "gæt.\n")
}

# Side 99.
gæt.et.tal()

# Afsnit 19.5.5.
matpow <- function(A, n) { # matrix-potensopløftning for heltal n>=0
  res <- diag(ncol(A))     # start: enhedsmatrix
  while (n > 0) {
    res <- res %*% A       # gang A på resultatet, n gange
    n <- n-1
  }
  res                      # det færdige resultat
}

# Side 100.
A <- matrix(c(1,2,3,4),2)
matpow(A,0)
matpow(A,1)
matpow(A,2)

matpow <- function(A, n) { # matrix-potensopløftning for heltal n
  if (n < 0) {
    n <- -n                # for negativ potens, vend n's fortegn
    A <- solve(A)          # og sæt A til den inverse af A
  }
  res <- diag(ncol(A))     # start: enhedsmatrix
  while (n > 0) {
    res <- res %*% A       # gang A på resultatet, n gange
    n <- n-1
  }
  res                      # det færdige resultat
}

matpow(A,-1)
matpow(A,-2)

# Afsnit 19.5.6. Side 101.
funpow <- function(f, n, x) {
  res <- x
  while (n > 0) {
    res <- f(res)
    n <- n-1
  }
  res
}

f <- function(s) { s + 0.04*s + 41000 }
funpow(f, 35, 0)

sapply(0:7, function(n) { funpow(f, n, 0) })

sapply(0:7, function(n) { c(n, funpow(f, n, 0)) } )

t(sapply(0:35, function(n) { c(n, funpow(f, n, 0)) } ))

A <- matrix(c(0.8, -0.1, 0.8, 0.4), 2)
b <- c(4, 2)
f <- function(v) { A %*% v + b }
f(c(15, 0))

# Side 102.
sapply(c(0:6), function(i) { funpow(f, i, c(15, 0)) } )

funpow <- function(f, n, x) {
  if (n==0) x else funpow(f, n-1, f(x))
}

# Afsnit 20.1. Side 103.
xx <- list(t=c(1.0, 1.5, 2.0, 2.5, 3.0),
           y=c(1.4, 0.3, -1.5, -3.1, -4.8))
xx

blandet <- list(A=matrix(c(1, 2, 3, 4), 2), d=42.2)

xx$t
xx$y
names(xx)

xx["y"]
xx[2]

xx$t <- 1:8
xx$ny <- 10
xx

# Afsnit 20.2. Side 104.
d <- data.frame(t=c(1.0, 1.5, 2.0, 2.5, 3.0),
                y=c(1.4, 0.3, -1.5, -3.1, -4.8))
d

d$ny <- 1:8

# Afsnit 20.3.
d[1,2]
d[2:4,2]
d[3,]
d[1:3,]

d1 <- d[ d$t > d$y ,]
d1

# Side 105.
d2 <- subset(d, t>y)
d2

# Afsnit 21. Side 106.
f <- function(x,y) { cos(x)*sin(2*y) }
x <- seq(0, 2*pi, len=50)
y <- seq(-pi, pi, len=50)
z <- outer(x, y, f)

# Afsnit 21.1.
persp(x, y, z)

persp(x, y, z, phi=35, theta=30)

persp(x, y, z, phi=35, theta=30, r=10)

persp(x, y, z, phi=35, theta=30, scale=FALSE)

persp(x, y, z, phi=35, theta=30, scale=FALSE, ticktype="detailed")

persp(x, y, z, phi=35, theta=30, col="red")
persp(x, y, z, phi=35, theta=30, col=rgb(1,0.5,0))

# Side 107.
persp(x, y, z, phi=35, theta=30, col="NA")
persp(x, y, z, phi=35, theta=30, col=c("red","green","blue"))
persp(x, y, z, phi=35, theta=30, col=rainbow(50))

persp(x, y, z, phi=35, theta=30, shade=0.7)
persp(x, y, z, phi=35, theta=30, shade=0.7, ltheta=110)

# Side 108.
persp(x, y, z, phi=35, theta=30, border="red")
persp(x, y, z, phi=35, theta=30, shade=0.7, border=NA)

# Afsnit 21.1.1.
zv <- z[1:49,1:49]
zv.skaleret <- (zv - min(zv)) / (max(zv) - min(zv))
zcol <- gray(zv.skaleret)
persp(x, y, z, phi=35, theta=30, col=zcol)

# Afsnit 21.2.
f <- function(x,y) { cos(x)*sin(2*y) }
x <- seq(0, 2*pi, len=50)
y <- seq(-pi, pi, len=50)
z <- outer(x, y, f)

image(x, y, z)

image(x, y, z, col=heat.colors(256))

# Side 109.
image(x, y, z, col=gray(seq(0,1,len=256)))

# Afsnit 21.3.
f <- function(x,y) { cos(x)*sin(2*y) }
x <- seq(0, 2*pi, len=50)
y <- seq(-pi, pi, len=50)
z <- outer(x, y, f)

contour(x, y, z)

# Side 110.
contour(x, y, z, nlevels=15)

contour(x, y, z, levels=c(-0.95,-0.75,-0.5,0,0.5,0.75,0.95))

contour(x, y, z, method="edge", col="red")

image(x, y, z)
contour(x, y, z, add=TRUE)

# Afsnit 22. Side 111.
x <- seq(0, 2, len=50)        # 50 støttepunkter i intervallet [0,2]
y <- seq(-1, 1, len=50)       # 50 støttepunkter i intervallet [-1,1]
z <- outer(x, y, f)           # Funktionsværdier i støttepunkterne
persp(x, y, z)                # Tegn grafen

overflade <- function(f, x.min, x.max, y.min, y.max) {
   x <- seq(x.min, x.max, len=50)  # 50 x-støttepunkter
   y <- seq(y.min, y.max, len=50)  # 50 y-støttepunkter
   z <- outer(x, y, f)             # Funktionsværdier i støttepunkterne
   persp(x, y, z)                  # Tegn grafen
}

f <- function(x,y) { cos(x)*sin(2*y) }
overflade(f, 0, 2, -1, 1)
overflade(f, 0, 10, -5, 5)
g <- function(x,y) { x^2+cos(x*y) }
overflade(g, 0, 2, 0, 2*pi)

overflade <- function(f, x.min, x.max, y.min, y.max, n) {
   x <- seq(x.min, x.max, len=n)   # n x-støttepunkter
   y <- seq(y.min, y.max, len=n)   # n y-støttepunkter
   z <- outer(x, y, f)             # Funktionsværdier i støttepunkterne
   persp(x, y, z)                  # Tegn grafen
}

overflade(f, 0, 2, -1, 1, 50)
overflade(f, 0, 10, -5, 5, 200)
overflade(g, 0, 2, 0, 2*pi, 100)

# Afsnit 22.1. Side 112.
overflade <- function(f, x.min=0, x.max=1, y.min=0, y.max=1, n=50) {
   x <- seq(x.min, x.max, len=n)   # n x-støttepunkter
   y <- seq(y.min, y.max, len=n)   # n y-støttepunkter
   z <- outer(x, y, f)             # Funktionsværdier i støttepunkterne
   persp(x, y, z)                  # Tegn grafen
}

overflade(f)
overflade(f, 0, 10, -5, 5)

# Afsnit 22.1.1.
overflade(f, n=100)

overflade(f, x.max=10, y.max=10)

overflade(n=100, x.max=5, y.max=5, f=g)

seq(1, 7, length.out=3)
seq(1, 7, len=3)
seq(1, 7, l=3)

# Afsnit 22.1.2. Side 113.
overflade <- function(f, x.min=min(x.interval), x.max=max(x.interval),
                      y.min=min(y.interval), y.max=max(y.interval),
                      n=50, x.n=n, y.n=n,
                      x.interval=c(0,1), y.interval=c(0,1) ) {
   x <- seq(x.min, x.max, len=x.n)
   y <- seq(y.min, y.max, len=y.n)
   z <- outer(x, y, f)
   persp(x, y, z)
}

overflade(f, x.max=10, y.interval=c(-5,5) )
overflade(g, 0, 2, 0, 2*pi, y.n=100)

overflade <- function(f, x.min=min(x.interval), x.max=max(x.interval),
                      y.min=min(y.interval), y.max=max(y.interval),
                      n=50, x.n=n, y.n=n,
                      x.interval=c(0,1), y.interval=c(0,1),
                      x=seq(x.min, x.max, len=x.n),
                      y=seq(y.min, y.max, len=y.n) ) {
   z <- outer(x, y, f)
   persp(x, y, z)
}

overflade(g, y.max=2*pi, x=c((0:10)/5,3:5))

# Afsnit 22.2. Side 114.
overflade <- function(f, x.min=min(x.interval), x.max=max(x.interval),
                      y.min=min(y.interval), y.max=max(y.interval),
                      n=50, x.n=n, y.n=n,
                      x.interval=c(0,1), y.interval=c(0,1),
                      x=seq(x.min, x.max, len=x.n),
                      y=seq(y.min, y.max, len=y.n),
                      ... ) {
   z <- outer(x, y, f)
   persp(x, y, z, ...)
}

overflade(g, x.max=4, y.max=2*pi, col="red")
overflade(f, 0, 10, 0, 10, col="blue", theta=30, phi=35)

# Afsnit C. Side 139.
preg <- lm(Kontrol ~ Tid + I(Tid^2), data=d)
preg

m <- seq(0, 1500, len=101)
plot(d$Tid, d$Kontrol, ylim=c(5, 6.75))
lines(m, predict(preg, data.frame(Tid=m)), col="red")

preg3 <- lm(Kontrol ~ Tid + I(Tid^2) + I(Tid^3), data=d)
lines(m, predict(preg3, data.frame(Tid=m)), col="green")
logreg <- lm(Kontrol ~ I(log(Tid)), data=d)
lines(m, predict(logreg, data.frame(Tid=m)), col="blue")

# Afsnit D. Side 140.
library("scatterplot3d")

x <- runif(100,0,10)
y <- runif(100,0,10)

z <- x + 2*y + 3 + rnorm(100)

scatterplot3d(x, y, z)

t <- seq(0, 8*pi, len=300)
x <- cos(t)
y <- sin(t)
z <- t
scatterplot3d(x, y, z, type="l")