［R］Silverman (1985)によるオートバイ衝突事故の実験データ

MASSパッケージにデータフレーム形式で含まれている。

> install.packages("MASS")
> library(MASS)
> head(mcycle)
  times accel
1   2.4   0.0
2   2.6  -1.3
3   3.2  -2.7
4   3.6   0.0
5   4.0  -2.7
6   6.2  -2.7
> tail(mcycle)
    times accel
128  52.0  10.7
129  53.2 -14.7
130  55.0  -2.7
131  55.0  10.7
132  55.4  -2.7
133  57.6  10.7
> nrow(mcycle)
[1] 133
> plot(mcycle$times, mcycle$accel)

引用元の論文は以下のとおり。

Silverman, B.W. (1985), Some Aspects of the Spline Smoothing Approach to Non‐Parametric Regression Curve Fitting. Journal of the Royal Statistical Society: Series B (Methodological), 47: 1-21.
doi:10.1111/j.2517-6161.1985.tb01327.x

［R］行列のトレース（固有和）を求める

Rには行列のトレースを直接計算する関数は標準では搭載されていない。

sum関数とdiag関数を組み合わせると簡単に求めることができる。

> d <- 1:9
> mx <- matrix(d, 3, 3)
> mx
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9
> sum(diag(mx))
[1] 15
> d <- 1:16
> mx <- matrix(d, 4, 4)
> mx
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16
> sum(diag(mx))
[1] 34

［R］同じ乱数を発生させる

同じ乱数を発生させるにはset.seed関数を使う。引数の値（種という）に応じて乱数の発生を設定するため、同じ値を設定すれば、同じ乱数が得られるようになる。

> runif(3)
[1] 0.3814675 0.2523344 0.4245695
> runif(3)
[1] 0.1408280 0.5990673 0.1032862
> runif(3)
[1] 0.1752793 0.1997167 0.5213036
> set.seed(3)
> runif(3)
[1] 0.1680415 0.8075164 0.3849424
> set.seed(3)
> runif(3)
[1] 0.1680415 0.8075164 0.3849424
> set.seed(4)
> runif(3)
[1] 0.585800305 0.008945796 0.293739612

R起動時には、自動で任意の種が設定される。

この関数を使用すると種に基づいて乱数が発生されるようになるが、この種は、Rを再起動しても同じ種を与えれば同じ乱数が発生される。機種が異なっても（WindowsとLinuxなど）、種さえ同じであれば全く同じ乱数が発生される。

［R］自然スプラインによるスプライン補間

splinefun関数を使う。Rの初期状態で使うことができる。以下、実行例。

> x <- c(1, 2, 4, 6, 8, 16)
> y <- c(0, 1, 2, 4, 3, 5)
> rfun <- splinefun(x, y, method = "natural")
> estx <- seq(min(x), max(x), length.out = 80)
> esty <- rfun(estx)
> plot(x, y)
> lines(estx, esty, col = "red")

splinefun関数のmethodオプションにnaturalを必ず指定すること。splinefun関数の戻り値は関数であり、スプライン補間による推定値の計算に使う。

［R］順列の総数を求める

Rには順列の総数nPrを求める関数が標準で搭載されていない。順列は以下のように書きかえることができるため、階乗を求めるfactorial関数を使用して計算する。

以下は、順列の総数 nPr を求める関数 permuを作成して計算した例。なお、n >= r、0! = 1であることに注意。n = r の場合は階乗そのものである。

₅P₂ = 20
₅P₃ = 60
₁₀P₃ = 720

> permu <- function(n, r) {factorial(n) / factorial(n - r)}
> permu(5, 2)
[1] 20
> permu(5, 3)
[1] 60
> permu(10, 3)
[1] 720
> permu(3, 3)
[1] 6

［R］組合せの総数を求める

choose関数を使う。以下、計算例。

₅C₂ = 10

₇C₃ = 35

> choose(5, 2)
[1] 10
> choose(7, 3)
[1] 35

［R］階乗を計算する

factorial関数を使う。以下、実行例。なお、4の階乗（4!）は24である。

> factorial(4)
[1] 24

ガンマ関数の値を計算するgamma関数を使っても計算できる。ガンマ関数の以下の性質を利用する。

Γ(n + 1) = n! (nは整数かつn≧1)

> gamma(4 + 1)
[1] 24

factorial関数とgamma関数はその実行速度が異なる。gamma関数のほうが断然早いので、gamma関数の使用を推奨する。

> system.time(for (i in 1:1000000) factorial(4))
   ユーザ   システム       経過  
      0.39       0.00       0.39 
> system.time(for (i in 1:1000000) gamma(4 + 1))
   ユーザ   システム       経過  
      0.11       0.00       0.11

［R］スプライン関数によるデータの平滑化

smooth.spline関数を使う。Rの初期状態で使うことができる。以下、実行例。

> x <- c(1, 2, 4, 6, 8, 16)
> y <- c(0, 1, 2, 4, 3, 5)
> r <- smooth.spline(y ~ x, all.knots = TRUE)
> estx <- seq(min(x), max(x), length.out = 80)
> esty <- predict(r, estx)$y
> plot(x, y)
> lines(estx, esty, col = "red")
> print(r)
Call:
smooth.spline(x = y ~ x, all.knots = TRUE)
Smoothing Parameter  spar= 0.5378882  lambda= 0.006212978 (14 iterations)
Equivalent Degrees of Freedom (Df): 2.943406
Penalized Criterion (RSS): 1.55995
GCV: 1.001812

 

［R］標準正規分布に従う確率変数（乱数）を得る

rnorm関数を使う。引数には、発生させる乱数の数を与える。

> rnorm(1)
[1] -1.35796
> rnorm(3)
[1] -0.5810260 -2.1525327  0.2157626

標準正規分布とは、平均が0、標準偏差が1の正規分布である。0付近の数値が多く得られ、正負の実数が得られる。

［R］行列の階数を得る

QR分解を行う関数qrの戻り値に階数がある（rank）。これを使用する。

> d1 <- rep(1, 9)
> ma1 <- matrix(d1, nrow = 3, byrow = TRUE)
> ma1
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    1    1    1
[3,]    1    1    1
> qr(ma1)$rank
[1] 1
> d2 <- c(1,1,1,0,1,0,0,0,1)
> ma2 <- matrix(d2, nrow = 3, byrow = TRUE)
> ma2
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    0    1    0
[3,]    0    0    1
> qr(ma2)$rank
[1] 3
> d3 <- c(1,1,1,0,0,1,0,0,0)
> ma3 <- matrix(d3, nrow = 3, byrow = TRUE)
> ma3
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    0    0    1
[3,]    0    0    0
> qr(ma3)$rank
[1] 2