リンク

無料ブログはココログ

R (数値計算)

2020年5月16日 (土)

[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)

Silverman1985

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

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

2020年5月14日 (木)

[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

2020年5月 3日 (日)

[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など)、種さえ同じであれば全く同じ乱数が発生される。

2019年9月19日 (木)

[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")

Nsp

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

2019年8月28日 (水)

[R]順列の総数を求める

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

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

5P2 = 20
5P3 = 60
10P3 = 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

2019年8月20日 (火)

[R]組合せの総数を求める

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

5C2 = 10

7C3 = 35

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

 

2019年7月 1日 (月)

[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

2019年6月25日 (火)

[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

Figure1 

[R]標準正規分布に従う確率変数(乱数)を得る

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

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

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

2019年3月13日 (水)

[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

より以前の記事一覧