-
Notifications
You must be signed in to change notification settings - Fork 1
/
lynxsunspots.R
58 lines (44 loc) · 1.27 KB
/
lynxsunspots.R
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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
# totally stolen from
# http://pbil.univ-lyon1.fr/members/lobry/corcau/
# -- Jesse, 2013
data(sunspots)
#library(ts)
data(lynx)
# setwd("")
spots <- window(sunspots, freq = 1, start = 1880, end = 1900)
lnx <- window(lynx, start = 1880, end = 1900)
ratio <- max(lnx)/max(spots)
par(mai = rep(1, 4))
png(filename="lynxsunspots.png",
res=300,
bg="white",
type="quartz",
pointsize=12,
width=6, height=5, units="in"
)
plot(lnx,
main = "Sunspot Count\nAnd Lynx Population Density",
t = "b",
ylab = ""
)
lines(ratio*spots, col = "red", t = "b" )
axis(side = 4, col = "red", col.axis = "red", at = ratio*pretty(spots),
lab = pretty(spots))
legend(1887, 4500, col = c("red", "black"), c("spots", "lynx"), pch = 21)
dev.off()
png(filename="lynxsunspotstotal.png",
res=300, bg="white",
type="quartz",
pointsize=12,
width=6, height=5, units="in"
)
spots <- window(sunspots, freq = 1, start = 1821, end = 1934)
ratio <- max(lynx)/max(spots)
plot(lynx, main = "Sunspot Count\nAnd Lynx Population Density",
ylab = "")
lines(ratio*spots, col = "red")
axis(side = 4, col = "red", col.axis = "red", at = ratio*pretty(spots),
lab = pretty(spots))
mtext("Sun spots", adj = 1, col = "red")
mtext("Lynx", adj = 0)
dev.off()