R – plot()

2016-09-08 / tau / コメントする

概要

Rの高水準作図関数で、タイプ指定で散布図、折れ線グラフなど様々なグラフを描くことができる。

基本形

数字の対を与えて散布図を描く。

x <- seq(0, 1, 0.1)
y <- x ^ 2
plot(x, y)

x <- seq(0, 1, 0.1)

y <- x ^ 2

plot(x, y)

r-plot-basic

タイプ指定

plot()のタイプ

引数	機能
`type="p"`	点プロット(デフォルト)
`type="l"`	折れ線グラフ
`type="b"`	点と線のプロット
`type="c"`	“b”で点を描かないプロット
`type="o"`	点プロットと線プロットの重ね描き
`type="h"`	各点からx軸まで垂線を描く
`type="s"`	左側の値から右にステップ表示
`type="S"`	右側の値から左にステップ表示
`type="n"`	軸だけ描く(その後低水準関数でプロット)

折れ線グラフ

例えばtype指定によって折れ線グラフを描ける。

x <- seq(0, 1, 0.1)
y <- x ^ 2
plot(x, y, type="l")

x <- seq(0, 1, 0.1)

y <- x ^ 2

plot(x, y, type="l")

r-plot-type-l

関数を直接与える場合

関数と定義域を与えて、直接その形を描ける。

plot(sin, -pi, pi)

1	plot(sin, -pi, pi)

r-plot-function-sin

定義した関数を渡すこともできる。

para <- function(x) x^2
plot(para, -1, 1, asp=1, pos=0)

1 2	para <- function(x) x^2 plot(para, -1, 1, asp=1, pos=0)

r-plot-function-user

asp=1は縦軸と横軸のスケール比を1:1に設定する引数
pos=0はゼロで軸が交わるようにする引数

R – 整数化関数

2016-09-08 / tau / コメントする

実数を整数化する関数には以下の通り。

`floor(x)`	ｘ以上でない最大の整数(正の数なら切り捨て)
`ceiling(x)`	ｘ未満でない最小の整数(正の数なら切り上げ)
`trunc(x)`	0に近い方へ整数化(切り捨て)
`round(x, digits = 0)`	IEEE式で丸める
`signif(x, digits = 6)`	digitsで指定された有効桁数に丸める

> floor(3.2); floor(-3.2)
[1] 3
[1] -4
> ceiling(3.2); ceiling(-3.2)
[1] 4
[1] -3
> trunc(3.2); trunc(-3.2)
[1] 3
[1] -3
> round(3.2); round(-3.2); round(3.7); round(-3.7)
[1] 3
[1] -3
[1] 4
[1] -4

> floor(3.2); floor(-3.2)

[1] 3

[1] -4

> ceiling(3.2); ceiling(-3.2)

[1] 4

[1] -3

> trunc(3.2); trunc(-3.2)

[1] 3

[1] -3

> round(3.2); round(-3.2); round(3.7); round(-3.7)

[1] 3

[1] -3

[1] 4

[1] -4

R – ベクトルのソート

2016-09-08 / tau / コメントする

sort()関数

sort()関数で引数にベクトルを指定して、その要素をソートできる。デフォルトで昇順、decreasing属性で降順にもできる。

> x <- c(5, 1, 4, 2, 3)
> sort(x)
[1] 1 2 3 4 5
> sort(x, decreasing=TRUE)
[1] 5 4 3 2 1

> x <- c(5, 1, 4, 2, 3)

> sort(x)

[1] 1 2 3 4 5

> sort(x, decreasing=TRUE)

[1] 5 4 3 2 1

partial属性は、指定した位置の要素より小さい要素をその要素より前に、大きい要素を後ろに移動させる(完全なソートはしない)。

以下の例では、3番目の要素(5)より小さな要素を前へ、大きな要素を後ろへ移動させている。5の前の要素の列、後ろの列は、それぞれ必ずしもソートされない。

> x <- c(3, 7, 5, 4, 2, 1, 6)
> sort(x, partial=3)
[1] 2 1 3 4 5 7 6

> x <- c(3, 7, 5, 4, 2, 1, 6)

> sort(x, partial=3)

[1] 2 1 3 4 5 7 6

partial属性をベクトルで与えた場合、複数の位置の要素の前後に各要素が集められる。以下の例では、5番目の要素(5)と10番目の要素(11)が基準となり、5より小さい要素群、5と11の間の要素群、11より大きい要素群に分けられるが、それぞれは必ずしもソートされない。

> x <- c(9,3,8,4,5,15,2,7,10,11,6,14,13,1,12)
> sort(x, partial=c(5, 10))
 [1]  1  3  2  4  5  6  8  7  9 10 11 14 13 15 12

> x <- c(9,3,8,4,5,15,2,7,10,11,6,14,13,1,12)

> sort(x, partial=c(5, 10))

[1] 1 3 2 4 5 6 8 7 9 10 11 14 13 15 12

ちなみにこの指定の仕方で、基準となる要素の大小関係が逆転していると、よくわからない結果になる。

> x <- c(9,3,8,4,5,15,2,7,10,11,6,14,13,1,12)
> sort(x, partial=c(1, 11))
 [1]  1  3  8  4  5  6  2  7  9 10 11 14 13 12 15

> x <- c(9,3,8,4,5,15,2,7,10,11,6,14,13,1,12)

> sort(x, partial=c(1, 11))

[1] 1 3 8 4 5 6 2 7 9 10 11 14 13 12 15

rank()関数

ベクトルの要素を昇順でソートしたときの順序を要素とするベクトルを返す。

> x <- c(8, 6, 10, 2, 4)
> rank(x)
[1] 4 3 5 1 2

> x <- c(8, 6, 10, 2, 4)

> rank(x)

[1] 4 3 5 1 2

R – ベクトルの演算

2016-09-08 / tau / コメントする

数値演算

ベクトル同士の数値演算

同じ長さのベクトル同士の演算は、要素ごとにその演算が実行される。

> c(1, 2, 3) + c(4, 5, 6)
[1] 5 7 9
> c(1, 2, 3) * c(4, 5, 6)
[1]  4 10 18

> c(1, 2, 3) + c(4, 5, 6)

[1] 5 7 9

> c(1, 2, 3) * c(4, 5, 6)

[1] 4 10 18

ベクトルの長さが違う場合は、短い方の要素が循環して利用される。ただし長い方の要素数が短い方の要素数の整数倍になっていないと警告メッセージが出る。

> c(10, 20, 30, 40) + c(1, 2)
[1] 11 22 31 42
> c(10, 20, 30, 40, 50) + c(1, 2)
[1] 11 22 31 42 51
 警告メッセージ: 
 c(10, 20, 30, 40, 50) + c(1, 2) で: 
   長いオブジェクトの長さが短いオブジェクトの長さの倍数になっていません

> c(10, 20, 30, 40) + c(1, 2)

[1] 11 22 31 42

> c(10, 20, 30, 40, 50) + c(1, 2)

[1] 11 22 31 42 51

警告メッセージ:

c(10, 20, 30, 40, 50) + c(1, 2) で:

長いオブジェクトの長さが短いオブジェクトの長さの倍数になっていません

ベクトルの各要素への同じ演算の適用

要素数1のベクトルを用いると、ベクトルの全要素に同じ演算を施すことができる。

> c(1, 2, 3, 4) * 2
[1] 2 4 6 8
> 1 / c(2, 3, 4, 5)
[1] 0.5000000 0.3333333 0.2500000 0.2000000

> c(1, 2, 3, 4) * 2

[1] 2 4 6 8

> 1 / c(2, 3, 4, 5)

[1] 0.5000000 0.3333333 0.2500000 0.2000000

比較演算

数値演算と同じく、要素ごとの比較演算も可能。要素数が違う場合も同じように循環適用される。

結果は論理型ベクトルとして返される。

> c(1, 2, 3, 4, 5) == c(1, 1, 4, 4, 5)
[1]  TRUE FALSE FALSE  TRUE  TRUE
> c(1, 1, 2, 2, 3, 3) == c(1, 2)
[1]  TRUE FALSE FALSE  TRUE FALSE FALSE
> -2:2 > 0
[1] FALSE FALSE FALSE  TRUE  TRUE

> c(1, 2, 3, 4, 5) == c(1, 1, 4, 4, 5)

[1] TRUE FALSE FALSE TRUE TRUE

> c(1, 1, 2, 2, 3, 3) == c(1, 2)

[1] TRUE FALSE FALSE TRUE FALSE FALSE

> -2:2 > 0

[1] FALSE FALSE FALSE TRUE TRUE

数学関数

数学関数の引数にベクトルを使うと、各要素にその関数が適用される。

> abs(-2:2)
[1] 2 1 0 1 2
> sqrt(1:4)
[1] 1.000000 1.414214 1.732051 2.000000

> abs(-2:2)

[1] 2 1 0 1 2

> sqrt(1:4)

[1] 1.000000 1.414214 1.732051 2.000000

集合演算

ベクトルを集合とみなして、集合演算を行うことができる。

> x <- c(-2, -1, 0, 1, 2)
> y <- c(0, 1, 2, 3, 4)
> union(x, y)      # 和集合
[1] -2 -1  0  1  2  3  4
> intersect(x, y)  # 積集合
[1] 0 1 2
> setdiff(x, y)    # 差集合
[1] -2 -1
> setequal(x, y)   # 集合として等しいか
[1] FALSE
> is.element(x, y) # xの各要素がyに含まれるか(x %in% y でもよい)
[1] FALSE FALSE  TRUE  TRUE  TRUE

> x <- c(-2, -1, 0, 1, 2)

> y <- c(0, 1, 2, 3, 4)

> union(x, y) # 和集合

[1] -2 -1 0 1 2 3 4

> intersect(x, y) # 積集合

[1] 0 1 2

> setdiff(x, y) # 差集合

[1] -2 -1

> setequal(x, y) # 集合として等しいか

[1] FALSE

> is.element(x, y) # xの各要素がyに含まれるか(x %in% y でもよい)

[1] FALSE FALSE TRUE TRUE TRUE

R – ベクトルの操作

2016-09-08 / tau / コメントする

要素の取得・部分ベクトルの抽出

単一要素の取得

ベクトルの要素は[]で添字を指定して参照

先頭要素は(0番目ではなく)1番目
要素数より大きい数を指定した場合はNAが返される
0を指定した場合はnumeric(0)～0を0個並べたベクトルが返される

> x <- c(11, 12, 13, 14, 15)
> x[1]
[1] 11
> x[2]
[1] 12
> x[5]
[1] 15
> x[6]
[1] NA
> x[0]
numeric(0)

> x <- c(11, 12, 13, 14, 15)

> x[1]

[1] 11

> x[2]

[1] 12

> x[5]

[1] 15

> x[6]

[1] NA

> x[0]

numeric(0)

複数要素からなるベクトルを抽出

複数の要素をベクトル化して指定した場合は、それらの位置の要素からなるベクトルが返される。返り値のベクトルの要素の順番は引数のベクトルで指定した順番になる。

> x <- c(11, 12, 13, 14, 15)
> x[c(5, 3, 1)]
[1] 15 13 11

> x <- c(11, 12, 13, 14, 15)

> x[c(5, 3, 1)]

[1] 15 13 11

要素を除外したベクトルを取得

負の数を指定した場合は、その絶対値の位置の要素を除いたベクトルが返される。

複数の要素を除外したい場合は、除外したい位置の添字をベクトル化して指定
除外要素の指定は順不同
返り値のベクトルの要素の順番は元のベクトルに準ずる
負数を使う場合は正数は混在できない
0は使えるが意味を持たない

> x[-1]
[1] 12 13 14 15
> x[c(-4, -2)]
[1] 11 13 15
> x[c(2, -4)]
 x[c(2, -4)] でエラー:  負の添字と混在できるのは 0 という添字だけです 
> x[c(0, -4)]
[1] 11 12 13 15

> x[-1]

[1] 12 13 14 15

> x[c(-4, -2)]

[1] 11 13 15

> x[c(2, -4)]

x[c(2, -4)] でエラー: 負の添字と混在できるのは 0 という添字だけです

> x[c(0, -4)]

[1] 11 12 13 15

条件を指定した要素・ベクトルの抽出

要素の代わりに条件式を指定すると、その条件に合った要素からなるベクトルが返される。

条件式を変数として書くこともできて、

条件を満たす要素位置がTRUE、満たさない要素位置がFALSEという要素からなるベクトルが返される
そのベクトルを元のベクトルの引数として指定すると、TRUEの位置の要素からなるベクトルが返される

> x <- c(-3, -2, -1, 0, 1, 2, 3)
> x[x >= 0]
[1] 0 1 2 3
> condition <- (x < 0)
> condition
[1]  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE
> x[condition]
[1] -3 -2 -1

> x <- c(-3, -2, -1, 0, 1, 2, 3)

> x[x >= 0]

[1] 0 1 2 3

> condition <- (x < 0)

> condition

[1] TRUE TRUE TRUE FALSE FALSE FALSE FALSE

> x[condition]

[1] -3 -2 -1

要素の追加・ベクトルの結合

既存のベクトルの長さ以上の要素を指定して値を代入すると、ベクトルが拡張されて値が代入される。

間が空いていればNAが代入される。

> a <- c(1)
> a
[1] 1
> a[2] <- 4
> a[3] <- 9
> a
[1] 1 4 9
> a[6] <- 36
> a
[1]  1  4  9 NA NA 36

> a <- c(1)

> a

[1] 1

> a[2] <- 4

> a[3] <- 9

> a

[1] 1 4 9

> a[6] <- 36

> a

[1] 1 4 9 NA NA 36

c()関数の中でベクトルと要素を列記して結合できる。複数の要素とベクトルを混在させることが可能。

> x <- c(-4, -3, -2)
> y <- c(2, 3, 4)
> c(-6, -5, x, -1, 0, 1, y, 5, 6)
 [1] -6 -5 -4 -3 -2 -1  0  1  2  3  4  5  6

> x <- c(-4, -3, -2)

> y <- c(2, 3, 4)

> c(-6, -5, x, -1, 0, 1, y, 5, 6)

[1] -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6

append(v₁, v₂, after=n)関数は、v₁のn番目の要素の後ろにv₂を挿入。

> x <- c(1, 2, 3, 4, 5)
> y <- c(-1, -2, -3)
> append(x, y)
[1]  1  2  3  4  5 -1 -2 -3
> append(x, y, after=3)
[1]  1  2  3 -1 -2 -3  4  5

> x <- c(1, 2, 3, 4, 5)

> y <- c(-1, -2, -3)

> append(x, y)

[1] 1 2 3 4 5 -1 -2 -3

> append(x, y, after=3)

[1] 1 2 3 -1 -2 -3 4 5

また、以下のように現在の要素数 + 1の要素に値をセットすると要素が追加される。

> (a <- c(1, 2, 3))
[1] 1 2 3
> a[length(a)+1] <- 4
> a
[1] 1 2 3 4

> (a <- c(1, 2, 3))

[1] 1 2 3

> a[length(a)+1] <- 4

> a

[1] 1 2 3 4

要素・ベクトルの置換

replace()関数で指定位置の要素を指定した要素で置き換える。

複数要素を置き換える場合は、指定位置と置き換え要素をベクトルで指定
- 指定位置の数と置き換え要素の数が等しいこと
指定位置の要素をベクトルで置き換えることはできない模様

> x <- c(10, 20, 30, 40, 50)
> replace(x, 3, 99)
[1] 10 20 99 40 50
> replace(x, c(2, 4), c(22, 44))
[1] 10 22 30 44 50
> replace(x, 3, c(-1, -2))
[1] 10 20 -1 40 50
 警告メッセージ: 
 replace(x, 3, c(-1, -2)) で: 
   置き換えるべき項目数が、置き換える数の倍数ではありませんでした

> x <- c(10, 20, 30, 40, 50)

> replace(x, 3, 99)

[1] 10 20 99 40 50

> replace(x, c(2, 4), c(22, 44))

[1] 10 22 30 44 50

> replace(x, 3, c(-1, -2))

[1] 10 20 -1 40 50

警告メッセージ:

replace(x, 3, c(-1, -2)) で:

置き換えるべき項目数が、置き換える数の倍数ではありませんでした

条件を指定して、その条件を満足する要素のみ置換することも可能。

> x <- c(-2, -1, 0, 1, 2)
> replace(x, (x<0), -9)
[1] -9 -9  0  1  2

> x <- c(-2, -1, 0, 1, 2)

> replace(x, (x<0), -9)

[1] -9 -9 0 1 2

R – ベクトルの生成

2016-09-08 / tau / コメントする

値を直接指定して生成する

c(...)の関数の中に要素を列挙する。

> c(1, 2, 3, 4, 5)
[1] 1 2 3 4 5

1 2	> c(1, 2, 3, 4, 5) [1] 1 2 3 4 5

同じ値が繰り返されるベクトル

要素が値の場合。

> rep(3,5)
[1] 3 3 3 3 3
> rep(3,times=5)
[1] 3 3 3 3 3

> rep(3,5)

[1] 3 3 3 3 3

> rep(3,times=5)

[1] 3 3 3 3 3

ベクトルを与えると、そのベクトルの要素が繰り返される。

> rep(c(1,2,3),3)
[1] 1 2 3 1 2 3 1 2 3
> rep(c(1,2,3),times=3)
[1] 1 2 3 1 2 3 1 2 3

> rep(c(1,2,3),3)

[1] 1 2 3 1 2 3 1 2 3

> rep(c(1,2,3),times=3)

[1] 1 2 3 1 2 3 1 2 3

ベクトルの要素を要素ごとに繰り返す。

> rep(c(1,2,3),each=4)
 [1] 1 1 1 1 2 2 2 2 3 3 3 3

1 2	> rep(c(1,2,3),each=4) [1] 1 1 1 1 2 2 2 2 3 3 3 3

要素ごとに繰り返しながら全体を繰り返す。timesとeachの指定順は関係しない。

> rep(c(1,2,3),times=3,each=2)
 [1] 1 1 2 2 3 3 1 1 2 2 3 3 1 1 2 2 3 3

1 2	> rep(c(1,2,3),times=3,each=2) [1] 1 1 2 2 3 3 1 1 2 2 3 3 1 1 2 2 3 3

一定値で増減するベクトル

1ずつ増減

[開始値]:[終了値]で指定する。

> 1:5
[1] 1 2 3 4 5
> 2:-2
[1]  2  1  0 -1 -2

> 1:5

[1] 1 2 3 4 5

> 2:-2

[1] 2 1 0 -1 -2

seq()関数を使っても同じ。

> seq(1:5)
[1] 1 2 3 4 5
> seq(2:-2)
[1] 1 2 3 4 5

> seq(1:5)

[1] 1 2 3 4 5

> seq(2:-2)

[1] 1 2 3 4 5

seq()関数の本来の表現。

> seq(1,5)
[1] 1 2 3 4 5
> seq(2,-2)
[1]  2  1  0 -1 -2

> seq(1,5)

[1] 1 2 3 4 5

> seq(2,-2)

[1] 2 1 0 -1 -2

1以外のステップ値で増減するベクトル

seq()関数の”by=”引数でステップ値を指定。

> seq(0,8,by=2)
[1] 0 2 4 6 8
> seq(8,0,by=-2)
[1] 8 6 4 2 0

> seq(0,8,by=2)

[1] 0 2 4 6 8

> seq(8,0,by=-2)

[1] 8 6 4 2 0

単に引数を3つ指定した場合は、3つ目の値が”by=”の値とみなされる。

> seq(0,12,3)
[1]  0  3  6  9 12

1 2	> seq(0,12,3) [1] 0 3 6 9 12

終了値がジャストでない場合は手前まで。

> seq(0,10,by=3)
[1] 0 3 6 9
> seq(0,-10,by=-3)
[1]  0 -3 -6 -9

> seq(0,10,by=3)

[1] 0 3 6 9

> seq(0,-10,by=-3)

[1] 0 -3 -6 -9

ただし要素の変化方向と増減値の符号が逆の場合はエラー。

> seq(0,5,by=-2)
 seq.default(0, 5, by = -2) でエラー:  'by' 引数中に誤った符号があります 
> seq(5,0,by=2)
 seq.default(5, 0, by = 2) でエラー:  'by' 引数中に誤った符号があります

> seq(0,5,by=-2)

seq.default(0, 5, by = -2) でエラー: 'by' 引数中に誤った符号があります

> seq(5,0,by=2)

seq.default(5, 0, by = 2) でエラー: 'by' 引数中に誤った符号があります

ベクトルの長さ(要素数)を指定

length.out属性で指定。

> seq(0,18,length.out=3)
[1]  0  9 18

1 2	> seq(0,18,length.out=3) [1] 0 9 18

length.outの代わりにlength、lenとしてもよい。

> seq(0,10,length=3)
[1]  0  5 10
> seq(0,10,len=4)
[1]  0.000000  3.333333  6.666667 10.000000
> seq_len(0,10,5)

> seq(0,10,length=3)

[1] 0 5 10

> seq(0,10,len=4)

[1] 0.000000 3.333333 6.666667 10.000000

> seq_len(0,10,5)

seq_len()関数は、非負の整数で長さを与えて数列を生成。

> seq_len(5)
[1] 1 2 3 4 5

1 2	> seq_len(5) [1] 1 2 3 4 5

sequence()関数

seq()関数とは違う。

引数に整数値を持つベクトルc(n₁, n₂, …)を指定して、各要素を終了値としたベクトル1:n₁、1:n₂、…を結合したベクトルを返す。

> sequence(4)
[1] 1 2 3 4
> sequence(c(1,2,3,4))
 [1] 1 1 2 1 2 3 1 2 3 4
> sequence(c(4,3,2,1))
 [1] 1 2 3 4 1 2 3 1 2 1
> sequence(c(1,3,5))
[1] 1 1 2 3 1 2 3 4 5

> sequence(4)

[1] 1 2 3 4

> sequence(c(1,2,3,4))

[1] 1 1 2 1 2 3 1 2 3 4

> sequence(c(4,3,2,1))

[1] 1 2 3 4 1 2 3 1 2 1

> sequence(c(1,3,5))

[1] 1 1 2 3 1 2 3 4 5

乱数列の生成

Rでは、様々な確率分布に従う乱数列を生成することができる。

基本形はr<name>(n, <parameters>)の形を取り、<name>で確率分布を、<parameters>で確率分布に対応したパラメータを指定する。

一様分布

runif(n, min=a, max=b)

(a, b)の一様分布に従う乱数をn個生成する。

> runif(10, -1, 1)
 [1]  0.3411880 -0.6890782  0.5073527  0.6238675 -0.5368339
 [6]  0.3966292  0.4342262 -0.1191653  0.5260462  0.3214433
<\pre>

> runif(10, -1, 1)

[1] 0.3411880 -0.6890782 0.5073527 0.6238675 -0.5368339

[6] 0.3966292 0.4342262 -0.1191653 0.5260462 0.3214433

<\pre>

二項分布

rbinom(n, size=s, prob=p)

確率pの二項分布に従う乱数を生成する。たとえば不良品率pの部品からs個のサンプルを取り出し、不良品の数を数えるという試行をn回繰り返す、などの例に相当する。

> rbinom(10, 1, 0.2)
 [1] 0 1 0 0 0 0 0 0 1 0
> rbinom(10, 1, 0.2)
 [1] 1 0 0 0 0 0 0 0 0 1
> rbinom(10, 5, 0.2)
 [1] 2 1 0 0 0 3 2 0 3 1
> rbinom(10, 5, 0.2)
 [1] 1 2 2 2 0 0 1 2 4 0

> rbinom(10, 1, 0.2)

[1] 0 1 0 0 0 0 0 0 1 0

> rbinom(10, 1, 0.2)

[1] 1 0 0 0 0 0 0 0 0 1

> rbinom(10, 5, 0.2)

[1] 2 1 0 0 0 3 2 0 3 1

> rbinom(10, 5, 0.2)

[1] 1 2 2 2 0 0 1 2 4 0

指数分布

rexp(n, rate=l)

l= $\lambda$ の指数分布に従う乱数を生成する。

ポアソン分布

rpois(n, lambda=l)

l= $\lambda$ のポアソン分布に従う乱数を生成する。

正規分布

rnorm(n, mean=m, sd=s)

平均m、標準偏差sの正規分布に従う乱数を生成する。

R – 演算子

2016-09-08 / tau / コメントする

算術演算子

四則演算は他の言語と同じ。正数除算(%/%)と剰余(%%)は独特。

演算の優先順位は、

()
^
%/%, %%
*, /
+, –

> 11 + 3
[1] 14
> 11 - 2
[1] 9
> 11 * 3
[1] 33
> 11 / 3
[1] 3.666667
> 11 ^ 3
[1] 1331
> 11 %/% 3
[1] 3
> 11 %% 3
[1] 2

> 11 + 3

[1] 14

> 11 - 2

[1] 9

> 11 * 3

[1] 33

> 11 / 3

[1] 3.666667

> 11 ^ 3

[1] 1331

> 11 %/% 3

[1] 3

> 11 %% 3

[1] 2

比較演算子

> 1 == 1.0
[1] TRUE
> 1/3 == 0.3
[1] FALSE
> 1/3 == 0.333333333333333
[1] FALSE
> 1/3 == 0.3333333333333333
[1] TRUE

> 1 == 1.0

[1] TRUE

> 1/3 == 0.3

[1] FALSE

> 1/3 == 0.333333333333333

[1] FALSE

> 1/3 == 0.3333333333333333

[1] TRUE

他の言語と同じだが、数値の同一性は一定の誤差より小さい場合に「等しい」と判定される。

論理演算子

数値に対しては他の言語と同じ

> !T; !F
[1] FALSE
[1] TRUE
> T & T; T & F; F & F
[1] TRUE
[1] FALSE
[1] FALSE
> T | T; T | F; F | F
[1] TRUE
[1] TRUE
[1] FALSE
> xor(T,T); xor(T,F); xor(F,F)
[1] FALSE
[1] TRUE
[1] FALSE

> !T; !F

[1] FALSE

[1] TRUE

> T & T; T & F; F & F

[1] TRUE

[1] FALSE

> T | T; T | F; F | F

[1] TRUE

[1] FALSE

> xor(T,T); xor(T,F); xor(F,F)

[1] FALSE

[1] TRUE

[1] FALSE

ポアソン過程について

2016-09-02 / tau / コメントする

概要

次のような条件で、客の到着、トランザクションやトラブルの発生といったランダム事象の発生を考える。

各事象は、その前の事象の発生から影響を受けたり、その次の事象に影響を与えたりしない(無記憶性)
各事象が単位時間あたりに発生する数には偏りがない(均質性)

この場合、事象の発生は時間軸上で一様分布と考えられ、定常ポアソン過程／ポアソン到着と呼ばれる。

ポアソン過程の到着数

単位時間当たりの到着数が $\lambda$ で与えられたとき、観測時間 $t$ の間に $k$ 回の到着が観測される確率は、以下のPoisson分布で与えられる。

(1) $\begin{equation*} P_n (X = k; t; \lambda) = \frac{(\lambda t)^k}{k!} e^{- \lambda t} \end{equation*}$

この導出は”ポアソン過程の到着数～ポアソン分布“に示すが、その概略は、観測時間 $t$ を有限・等間隔の微小時間間隔 $\Delta t$ に分け、到着数が $k$ 回となる組み合わせの数を考慮しながら確率を計算し、 $\Delta t \to \infty$ の極限をとるという操作。

このとき、到着回数の平均と分散は以下のようになるが、その導出はPoisson分布の説明を参照。

(2) $\begin{equation*} E(k) = \lambda = {\rm const} \end{equation*}$

(3) $\begin{equation*} V(k) = \lambda = {\rm const} \end{equation*}$

ポアソン過程の到着時間間隔

単位時間当たりの到着数が $\lambda$ で与えられたとき、ある到着から次の到着までの時間間隔 $t$ は指数分布に従い、その確率密度 $f(t)$ と確率分布 $F(t \le T)$ は以下のように与えられる。

(4) $\begin{equation*} f(t) = \lambda e^{- \lambda t} \end{equation*}$

(5) $\begin{equation*} F(t \le T) = 1 - e^{- \lambda T} \end{equation*}$

この導出は”ポアソン過程の到着間隔～指数分布“に示すが、その概略は上のPoisson分布と同じく、観測時間 $t$ を微小時間間隔に分けて、その極限をとるという操作。

このとき、到着時間間隔の平均と分散は以下のようになるが、その導出は指数分布の説明を参照。

(6) $\begin{equation*} E(k) = \frac{1}{\lambda} \end{equation*}$

(7) $\begin{equation*} V(k) = \frac{1}{\lambda ^2} \end{equation*}$

このほか、それぞれの確率分布の形状や考察についても、Poisson分布の説明、指数分布の説明にまとめた。

ポアソン過程のシミュレーションと確認

一様分布による到着のシミュレーション

”ポアソン過程～一様分布によるシミュレーション“で、一様分布に従う乱数を発生させて、その到着時間間隔／到着回数が指数分布／Poisson分布に従っているか確認した。

ポアソン過程～一様分布によるシミュレーション

2016-09-02 / tau / コメントする

概要

一様分布をとる乱数によってランダムな到着をシミュレートし、ポアソン過程の時間間隔や到着回数を理論上の確率分布と比較してみた。

手順

表計算ソフトで0 ～1の間で平均して10個の発生するように乱数を発生させる
- 0～1で10個
- 0～10で100個
- 0～100で1000個
その値を昇順にソート(到着事象の到着時刻順にあたる)
到着時刻の分布を確認
到着時間間隔の分布と指数分布の比較
1. 各到着時刻間の時間間隔を計算
2. 時間間隔データの頻度(確率密度に相当)を棒グラフで表示し、指数分布の密度関数のグラフを重ねて比較
3. 時間間隔データの累積頻度(確率分布に相当)を棒グラフで表示し、指数分布関数のグラフを重ねて表示
到着回数の分布とPoisson分布の比較

指数分布のデータを得る手順については、 (参考)データ数10個の場合到着時間間隔の計算手順を参照。

Poisson分布のデータを得る手順については、(参考)データ数100個の場合の到着数の計算手順を参照。

以下、データ数10、100、1000の場合の到着時刻の分布、到着間隔の確率密度と確率分布について比較していく。

到着時刻の分布の比較

データ数10、100、1000に対する到着時刻の分布を以下に示す。

データ数が多くなるほどグラフは平滑になるが、1000個と比較的多くのデータでもばらつきが見られる。

データ数：10
random-arrival-verification-duration0010

データ数：100
random-arrival-verification-duration0100

データ数：1000
random-arrival-verification-duration1000

到着時間間隔と指数分布の比較

頻度(確率密度)

データ数10、100、1000に対する到着時間間隔の頻度分布とこれに相当する指数分布の確率密度関数を比較する。

到着時刻の分布に比べてデータ数100からシミュレーション値と理論値がよく整合している。

データ数：10
random-arrival-verification-exp-density0010

データ数100
random-arrival-verification-exp-density0100

データ数1000
random-arrival-verification-exp-density1000

累積頻度(確率分布)

同じデータの累積頻度とこれに相当する指数分布の確率分布関数を比較する。

確率分布に関しては、データ数10と少ないデータでも理論値と類似している。

データ数：10
random-arrival-verification-exp-prob0010

データ数：100
random-arrival-verification-exp-prob0100

データ数：1000
random-arrival-verification-exp-prob1000

到着数の分布とPoisson分布の比較

到着数の分布をみるためには、データ数10個の場合はデータセットにならないため、データ数100、1000で比較してみる。データ数100を10組に、データ数100は100組のデータセットを得ることができる。

到着数の分布と理論上のPoisson分布との整合性は、データセット数100でもそれほどよくない。

データ数：100(データセット数：10)
random-arrival-verification-poisson-prob0100

データ数：1000(データセット数100)
random-arrival-verification-poisson-prob1000

(参考)データ数10個の場合到着時間間隔の計算手順

以下の表は、0～1の範囲で10個の乱数を発生させ、それらを昇順にソートして発生間隔を計算したもの。

statistics-random-arrival-verification-1

まず、これらの発生時刻を、一定の階級幅で度数分布に分ける。データ数が10個の場合はばらつきが大きいが、データ数が多くなると、これが一様分布に近づくと期待される。

statistics-random-arrival-verification-2

次に、到着時間間隔の度数分布とその頻度を計算する。これを指数分布の密度関数と比較するため、時間間隔dtで除している。一番右の欄は、到着率10の指数分布の密度関数の値。

statistics-random-arrival-verification-3

最後に、指数分布の確率関数と比較をするため、到着間隔の累積度数と頻度を計算する。一番右の欄は、到着率10の確率分布関数の値。

statistics-random-arrival-verification-4

(参考)データ数100個の場合の到着数の計算手順

元データは、到着時刻を昇順にソートしたもの。

random-arrival-verification-poisson1

これらのデータの、0以上1未満、1以上2未満、・・・・、9以上10未満の10組ごとのデータ数を数える。

random-arrival-verification-poisson2

これらについて、到着回数ごとの頻度分布に整理する。下表一番右の欄は、到着率10に対するPoisson確率値。

random-arrival-verification-poisson3

ポアソン過程の到着数～ポアソン分布

2016-09-02 / tau / コメントする

単位時間当たりの到着率 $\lambda$ のポアソン過程において、 $t$ 時間の間に $k$ 回の到着が発生する確率を考える。

$t$ を $n$ 等分し、 $\Delta t = t / n$ とする。この間に $k$ 回の到着が発生し、 $(n - k)$ 回は到着が発生しないとすると、その確率は次の二項分布で表される。

(1) $\begin{equation*} P(k;t;n) = {}_n C_k (\lambda \Delta t)^k (1 - \lambda \Delta t)^{n-k} \end{equation*}$

この式を、以下のように展開しておく。

(2) $\begin{eqnarray*} P(k;t;n) &=& \frac{n!}{k! (n-k)!} (\lambda \Delta t)^k (1 - \lambda \Delta t)^{n - k} \\ &=& \frac{n(n - 1) \cdots (n - k + 1)}{k!} \left( \frac{\lambda t}{n} \right)^k \left(1 - \frac{\lambda t}{n} \right)^{n - k} \end{eqnarray*}$

ここで $n \to \infty (\Delta t \rightarrow 0)$ の極限を考える。まず前2項については、

(3) $\begin{eqnarray*} &&\lim_{n \to \infty} \frac{n(n - 1) \cdots (n - k + 1)}{k!} \left( \frac{\lambda t}{n} \right)^k \\ &=& \lim_{n \to \infty} \frac{1 \cdot \left(1 - \cfrac{1}{n} \right) \cdots \left(1 - \cfrac{k - 1}{n} \right)}{k!} (\lambda t)^k \\ &=& \frac{(\lambda t)^k}{k!} \end{eqnarray*}$

また3項目については $-n / \lambda t = r \to \infty$ とおいて、

(4) $\begin{equation*} \lim_{r \to \infty} \left( 1 + \frac{1}{r} \right)^{- r \lambda t - k} = \lim_{r \to \infty}\left[ \left(1 + \frac{1}{r} \right)^{- r \lambda t} \left( 1 + \frac{1}{r} \right)^{-k} \right] = e^{- \lambda t} \end{equation*}$

通常、Poisson分布の表現は上式において $\lambda t \to \lambda$ として表現されるが、この場合の $\lambda$ は観測時間内の平均到着数として考える。ここでは、単位時間あたりの到着率 $\lambda$ と観測時間 $t$ を明確にするため $\lambda t$ と表現した。

以上から、到着率 $\lambda$ のポアソン過程において、観測時間 $t$ の間に $k$ 回の到着が発生する確率は以下のPoisson分布で与えられる。

(5) $\begin{equation*} P(k; t; \lambda) = \frac{(\lambda t)^k}{k!} e^{- \lambda t} \end{equation*}$