［R］周期変動を含むデータの回帰モデルを作成する

周期変動を含むデータの回帰モデルを作成する場合、sin関数とcos関数を組み合わせたモデルで表現することができる。ここでは、以下の式を考える。

y = A + B * x + C * sin(2 * π * x) + D * cos(2 * π * x)

A,B,C,D:求めるモデルパラメータ

最小2乗法により、A～Dの4つのモデルパラメータ（定数）を決定する。なお、ここではその周期変動の周期は1とする。ここでsin関数、cos関数の引数を2πxとしているところが重要。sin関数、cos関数共に2πラジアンであり今は周期を1と仮定しているためにそのような指定をしている。xが1増えれば2π増えるが元に戻る、という考え方。

サンプルデータとその散布図として、以下のデータを扱う。横軸データがxベクトル、縦軸データがyベクトルである。データを散布図で示し、プロット点を青実線でつないでいる。

> x <- c(2014.0, 2014.2, 2014.5, 2014.7, 2015.2, 2015.6, 2016.1)
> y <- c(0, 0.2, 0, -0.2, 0.15, -0.05, 0.2)
> plot(x, y, type = "p")
> lines(x, y, col = "blue")

サンプルデータをプロットしたところで、さっそく周期変動を含む回帰モデルを作成し、それを線で描画してみる。

> r <- lm(y ~ x + sin(2 * pi * x) + cos(2 * pi * x))
> xet <- seq(min(x), max(x), by = 0.01)
> yet <- predict(r, data.frame(x = xet))
> lines(xet, yet, col = "red")

このようにして、近似曲線が得られた。

［サンプルデータ］書籍「統計解析のはなし」（石村貞夫著、東京図書）の表1.1.1（pp.2-8）

以下に統計解析の専門書である「統計解析のはなし」（石村貞夫著、東京図書）の、2～8ページに掲載されている表1.1.1のデータを掲載する。

本掲載は同書籍の利便性を高めるため私的利用を目的に掲載をしているものである。同書籍にて統計解析を学んだ者として同書籍に敬意を払っており、恣意的に著作権等を侵害しているわけではないことを明記しておく。

なお、あくまで私的利用のために手入力をしたものであり、本データについてはその掲載内容との整合性について保証するものではない。

以下にデータをテキスト形式のファイル（data.txt）と、その内容を掲載する。リンク先のファイルをダウンロードしてもよいし、当該範囲内をマウスで選択肢、メモ帳に貼り付ければよい。データの内容は以下のとおり。

T大学女子学生（1年生）200人の身長(cm)、体重(kg)、バスト(cm)、ウェスト(cm)、ヒップ(cm)、その女子学生が結婚相手に望む男性の身長(cm)

data.txt

no, height, weight, bust, waist, hip, male
1, 151, 48, 81, 60, 86, 175
2, 154, 44, 75, 60, 78, 176
3, 160, 48, 80, 58, 85, 178
4, 160, 52, 86, 63, 91, 180
5, 163, 58, 90, 66, 92, 172
6, 156, 58, 80, 68, 93, 175
7, 158, 62, 83, 66, 91, 172
8, 156, 52, 85, 66, 90, 170
9, 154, 45, 80, 61, 88, 170
10, 160, 55, 80, 63, 90, 180
11, 154, 54, 85, 63, 90, 170
12, 162, 47, 80, 58, 90, 175
13, 156, 43, 80, 63, 87, 180
14, 162, 53, 78, 60, 85, 180
15, 157, 54, 83, 63, 90, 175
16, 162, 64, 94, 68, 92, 175
17, 162, 47, 80, 63, 88, 178
18, 169, 61, 80, 66, 93, 185
19, 150, 38, 75, 62, 83, 178
20, 162, 48, 83, 59, 85, 178
21, 154, 47, 77, 63, 87, 180
22, 152, 58, 85, 66, 93, 170
23, 161, 46, 80, 58, 84, 175
24, 160, 47, 78, 59, 80, 172
25, 160, 45, 80, 56, 80, 170
26, 153, 40, 80, 56, 81, 170
27, 155, 40, 80, 58, 83, 170
28, 163, 55, 80, 63, 87, 170
29, 160, 62, 84, 63, 95, 176
30, 159, 50, 85, 58, 80, 175
31, 164, 50, 81, 62, 88, 177
32, 158, 46, 78, 63, 83, 176
33, 150, 45, 85, 58, 80, 178
34, 155, 49, 82, 63, 82, 178
35, 157, 53, 82, 63, 88, 170
36, 161, 57, 85, 61, 92, 170
37, 168, 60, 80, 66, 93, 175
38, 162, 55, 80, 63, 88, 175
39, 153, 47, 80, 63, 88, 175
40, 154, 50, 80, 60, 85, 170
41, 158, 53, 79, 61, 90, 180
42, 151, 46, 82, 61, 85, 170
43, 155, 50, 82, 60, 84, 175
44, 155, 45, 80, 60, 85, 180
45, 165, 50, 80, 60, 85, 185
46, 165, 51, 82, 62, 90, 178
47, 154, 48, 82, 60, 90, 178
48, 148, 48, 78, 62, 90, 170
49, 169, 55, 82, 61, 87, 185
50, 158, 54, 87, 63, 88, 178
51, 146, 43, 82, 60, 87, 175
52, 166, 63, 93, 63, 98, 176
53, 161, 53, 82, 62, 90, 180
54, 143, 42, 80, 60, 80, 170
55, 156, 46, 80, 57, 92, 172
56, 156, 69, 98, 72, 98, 170
57, 149, 47, 75, 60, 82, 170
58, 162, 48, 76, 61, 90, 180
59, 159, 50, 80, 66, 85, 182
60, 164, 55, 80, 63, 90, 178
61, 162, 45, 78, 60, 84, 180
62, 167, 49, 80, 60, 84, 178
63, 159, 51, 83, 60, 86, 178
64, 153, 51, 82, 63, 85, 180
65, 146, 44, 83, 61, 91, 175
66, 156, 58, 95, 65, 90, 180
67, 160, 53, 84, 63, 94, 175
68, 158, 48, 85, 58, 85, 175
69, 151, 46, 83, 59, 85, 175
70, 157, 48, 78, 63, 83, 182
71, 151, 43, 79, 57, 83, 171
72, 156, 50, 82, 58, 85, 175
73, 166, 58, 86, 66, 90, 175
74, 159, 49, 80, 63, 83, 170
75, 157, 50, 79, 60, 84, 175
76, 156, 47, 80, 62, 84, 170
77, 159, 47, 80, 60, 85, 180
78, 156, 52, 82, 63, 85, 172
79, 156, 47, 80, 61, 87, 175
80, 161, 50, 80, 60, 88, 173
81, 151, 51, 82, 63, 89, 176
82, 162, 53, 80, 60, 90, 175
83, 153, 45, 76, 60, 80, 175
84, 157, 51, 80, 61, 89, 178
85, 153, 57, 90, 66, 93, 175
86, 159, 56, 85, 64, 98, 185
87, 157, 52, 82, 60, 88, 170
88, 158, 52, 85, 63, 88, 175
89, 159, 51, 83, 63, 89, 175
90, 159, 48, 80, 58, 87, 175
91, 159, 49, 83, 60, 88, 175
92, 153, 45, 78, 58, 80, 175
93, 153, 45, 80, 59, 83, 170
94, 164, 50, 80, 60, 82, 175
95, 157, 53, 80, 63, 93, 172
96, 157, 45, 78, 59, 80, 175
97, 155, 56, 86, 63, 90, 170
98, 149, 53, 86, 66, 93, 176
99, 160, 52, 82, 62, 92, 175
100, 150, 46, 80, 63, 85, 173
101, 161, 48, 80, 60, 85, 178
102, 155, 51, 82, 61, 87, 175
103, 156, 48, 80, 58, 83, 178
104, 156, 50, 82, 60, 90, 173
105, 151, 43, 80, 59, 82, 173
106, 153, 50, 79, 62, 89, 170
107, 161, 54, 85, 63, 85, 180
108, 165, 58, 88, 63, 90, 180
109, 155, 49, 78, 62, 85, 172
110, 156, 50, 80, 63, 88, 180
111, 165, 58, 88, 63, 90, 181
112, 156, 48, 82, 60, 80, 178
113, 154, 48, 80, 62, 90, 173
114, 160, 55, 81, 63, 92, 173
115, 159, 59, 85, 63, 94, 180
116, 155, 56, 82, 62, 93, 165
117, 165, 63, 83, 64, 85, 175
118, 151, 45, 83, 60, 85, 173
119, 156, 48, 80, 60, 80, 180
120, 158, 44, 78, 58, 80, 180
121, 155, 48, 81, 61, 85, 172
122, 162, 51, 86, 63, 90, 175
123, 163, 48, 85, 62, 89, 180
124, 153, 51, 78, 63, 82, 175
125, 157, 44, 80, 58, 83, 175
126, 153, 50, 82, 60, 88, 175
127, 151, 40, 79, 57, 80, 178
128, 153, 48, 83, 63, 85, 178
129, 153, 47, 82, 61, 86, 175
130, 153, 48, 83, 60, 84, 180
131, 162, 52, 80, 60, 90, 170
132, 156, 50, 80, 62, 85, 175
133, 156, 49, 80, 63, 90, 170
134, 153, 50, 85, 63, 90, 178
135, 161, 52, 78, 60, 85, 180
136, 156, 51, 78, 60, 88, 180
137, 164, 52, 96, 61, 90, 177
138, 155, 45, 79, 63, 85, 170
139, 159, 48, 78, 60, 84, 175
140, 168, 70, 88, 70, 90, 183
141, 159, 47, 75, 59, 88, 175
142, 170, 56, 82, 58, 88, 180
143, 164, 67, 90, 66, 92, 180
144, 168, 56, 80, 61, 90, 180
145, 169, 60, 81, 63, 93, 185
146, 156, 40, 80, 56, 82, 180
147, 156, 50, 79, 63, 90, 175
148, 158, 43, 78, 57, 82, 163
149, 155, 54, 80, 63, 87, 175
150, 155, 50, 83, 62, 88, 175
151, 152, 50, 80, 65, 85, 175
152, 157, 47, 75, 62, 83, 170
153, 153, 45, 80, 60, 88, 170
154, 166, 52, 80, 60, 85, 178
155, 155, 46, 80, 60, 86, 175
156, 163, 58, 85, 69, 90, 175
157, 150, 47, 80, 59, 86, 175
158, 160, 53, 80, 63, 90, 180
159, 153, 47, 83, 63, 86, 175
160, 168, 54, 85, 62, 88, 175
161, 165, 56, 85, 63, 88, 182
162, 159, 53, 84, 62, 87, 180
163, 153, 47, 82, 60, 85, 173
164, 151, 47, 88, 60, 88, 173
165, 156, 48, 80, 60, 86, 175
166, 158, 50, 86, 62, 88, 175
167, 155, 48, 83, 60, 89, 180
168, 165, 50, 85, 58, 80, 180
169, 157, 52, 82, 60, 89, 172
170, 160, 54, 82, 63, 88, 180
171, 163, 45, 83, 60, 92, 175
172, 163, 53, 86, 65, 90, 173
173, 162, 48, 85, 63, 85, 185
174, 167, 57, 83, 63, 85, 178
175, 160, 53, 83, 63, 85, 180
176, 162, 43, 80, 60, 85, 175
177, 158, 50, 85, 60, 85, 175
178, 159, 43, 78, 62, 84, 175
179, 159, 56, 71, 65, 87, 181
180, 160, 47, 85, 60, 85, 180
181, 162, 51, 80, 63, 87, 175
182, 162, 54, 80, 63, 87, 180
183, 158, 53, 84, 63, 87, 172
184, 147, 45, 80, 62, 90, 173
185, 164, 60, 87, 66, 94, 178
186, 147, 38, 85, 58, 80, 173
187, 156, 43, 80, 56, 84, 173
188, 161, 60, 92, 63, 95, 178
189, 160, 46, 78, 60, 84, 180
190, 154, 53, 84, 64, 85, 170
191, 163, 54, 83, 63, 92, 180
192, 158, 45, 83, 57, 89, 165
193, 157, 54, 85, 63, 90, 170
194, 164, 63, 85, 66, 85, 178
195, 152, 41, 79, 57, 86, 175
196, 157, 50, 82, 60, 90, 170
197, 155, 50, 80, 63, 88, 180
198, 163, 56, 82, 66, 90, 180
199, 155, 51, 83, 62, 90, 177
200, 152, 42, 78, 58, 86, 176

書籍の奥付によれば著者は鶴見大学の助教授（掲載時）とあるため、著者が講義で鶴見大学の女子学生（1年生）にアンケートを採った結果と思われる。

［R］年月日の数値を日付（Dateオブジェクト）に簡単に変換する

以下の2通りの方法がある。いずれの方法もベクトルを扱うことができる。なお、ここではDateオブジェクトに変換することを考えており、時刻やタイムゾーンの情報は扱わない。

（1）as.Date関数とsprintf関数を組み合わせる使う
（2）as.Date関数とISOdate関数を組み合わせて使う

■（1）の方法

> year <- c(2014, 2015)
> month <- c(1, 2)
> day <- c(11, 12)
> as.Date(sprintf("%d-%d-%d", year, month, day))
[1] "2014-01-11" "2015-02-12"

■（2）の方法

> year <- c(2014, 2015)
> month <- c(1, 2)
> day <- c(11, 12)
> as.Date(ISOdate(year, month, day))
[1] "2014-01-11" "2015-02-12"

ISOdate関数の戻り値はPOSIXctオブジェクト（年月日、時刻、タイムゾーンの情報を持つ）であることから、年月日しか扱わないDateオブジェクトにするには、as.Date関数でDateオブジェクトにする必要がある。

> class(ISOdate(2015, 1, 1))
[1] "POSIXct" "POSIXt"
> class(as.Date(ISOdate(2015, 1, 1)))
[1] "Date"

なお、処理は（2）の方法が断然早い。

> year <- seq(0, 9999, 1)
> month <- 1
> day <- 2
> system.time(as.Date(sprintf("%d-%d-%d", year, month, day)))
   ユーザ   システム       経過 
      0.11       0.00       0.11
> system.time(as.Date(ISOdate(year, month, day)))
   ユーザ   システム       経過 
      0.06       0.00       0.06

［R］コマンドの実行時間を計る

system.time関数を使えばよい。実行時間を計りたいコマンドをsystem.time関数の「( )」に挟んで記述して関数を実行すればよい。

以下は、正規分布をする乱数を1万個作成することを10万回繰り返したときの（長さが1万×10万＝10億のベクトルを作る）、実行時間を計測した例。

> system.time(d <- rep(rnorm(1e4), 1e5))
   ユーザ   システム       経過 
      1.33       2.41       3.86

表示される時間の単位は秒（secs）。精度は1/100秒（10ミリ秒）。表示されるそれぞれの時間については以下のとおり。

ユーザ：記述をしたコマンドの処理時間
システム：システム(OSなど)に依る処理時間
経過：コマンドの起動から終了までに要した時間

あくまでその時に計測をしつつコマンドを実行したときの計測結果を表示しているに過ぎない。そのため、実行するごとに出力は当然変わる。「ユーザ＋システム＝経過」となるはずだが、必ずしもそうはならない。

［R］プロット点の縁と中を別々の色で塗りつぶす

オプションpchに21～25の値を与えて描画すればよい。21～25はそれぞれ丸（○）、四角（□）、菱形（◇）、上向き三角（△）、下向き三角（▽）である。その際の縁の色はcolオプションで、中を塗りつぶす色はbgオプションで指定をする。以下、描画の例。

xl <- c(0, 5)
yl <- c(0, 5)
plot(0, 0, xlim = xl, ylim = yl, type = "n", xlab = "", ylab = "")
# まず直線を描画
segments(0, 0, 5, 5)
# 黒実線の円、中は透過
points(1, 1, pch = 1, cex = 3)
# 赤く塗りつぶされた丸
points(2, 2, pch = 16, cex = 3, col = "red")
#　縁は赤実線、中は白塗りつぶしの丸
points(3, 3, pch = 21, cex = 3, col = "red", bg = "white")
# 縁の線を太くして、中は青で塗りつぶしてみる
points(4, 4, pch = 21, cex = 3, lwd = 4, col = "red", bg = "blue")

実行結果は以下のとおり。 一般に近似曲線などはプロット点の下に描画したほうが見栄えは良く、図中の一番左下のプロット点は美しくない。左下から2番目や3番目のように描画するべきである。4番目はあくまで描画例であり、実際の図でこのような色の組み合わせは使ってはいけない。

［R］尖度の計算

資料の分布のとがり具合の度合いを表すものに尖度というものがある。

Rに標準で組み込まれている組み込み関数に尖度を計算させるものはないため、計算するには組み込み関数を組み合わせて計算する必要がある。

以下は自作をした尖度を計算する関数。第1仮引数valに尖度を計算させたい分布を数値ベクトルで与えると、尖度を関数の戻り値としている。

kurtosis <- function(val) {
    n <- length(val)
    var1 <- var(val) * (n - 1) / n  # 不偏分散から標本分散を求める
    sd1 <- sqrt(var1)               # 標本分散の平方根
    r <- sum((val - mean(val)) ^ 4) / (n * sd1 ^ 4)
    return(r)
}

なお、求めた尖度をa4とすると、a4の値により、資料の分布について以下のことが言える。

a4 > 3 →分布は正規分布よりとがっている
a4 = 3 →分布は正規分布と同じような形をしている
a4 < 3 →分布は正規分布よりなだらか