一様乱数によるポアソン過程のシミュレーション

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

一様乱数に基づく分析

R言語を使って、到着率をλ、観測時間をtとしてその間の到着時刻を一様乱数で発生させ、それらの分布や到着時間間隔、到着回数の分布などを確認する。

到着時刻の分布

以下の手順で到着時刻の乱数データを得てヒストグラムを描き、結果の一様性を確認する。

一様分布による到着確率の考え方の概要やその他の考え方についてはRによるポアソン過程のシミュレーションを参照。

観測時間tの間の平均到着数はλt
ただしtの間で個数を指定して一様乱数を発生させると、全体の到着数を縛ることになる
そこでtと同じ長さの余裕分をtの前後にとり、区間[–t, 2t)の間で3λt個の乱数を発生させる(到着時間の確率変数)
そのデータから、区間[0, t)の範囲にあるデータを抽出し、ヒストグラムを描く
その際、データの階級のレンジ～1レンジ幅の平均到着数も変化させてみる

到着率λ = 1/10[回/秒]、観測時間5000[秒]、階級レンジ100秒としたときの、Rのコードは以下の通り。このとき各階級の平均的な到着回数は1/10 * 100 = 10回程度と予想される。

lambda <- 1 / 10
t.obs <- 5000
rank.width <- 100
num.data <- lambda * t.obs

org.data <- (runif(floor(num.data * 3)) * 3 - 1) * t.obs
data <- org.data[(org.data >= 0) & (org.data < t.obs)]

ranks <- seq(0, t.obs, rank.width)
hist(data, breaks=ranks)

lambda <- 1 / 10

t.obs <- 5000

rank.width <- 100

num.data <- lambda * t.obs

org.data <- (runif(floor(num.data * 3)) * 3 - 1) * t.obs

data <- org.data[(org.data >= 0) & (org.data < t.obs)]

ranks <- seq(0, t.obs, rank.width)

hist(data, breaks=ranks)

各変数の意味は、

`lambda`	到着率λ
`t.obs`	観測時間t
`rank.width`	ヒストグラムの階級幅(時刻幅)。この幅の中の平均到着数はλ*rank_width回。
`org.data`	前後に余裕幅を持たせた元データ
`data`	t.obsの範囲内だけを切り出したデータ

まず、rank.width = 100すなわち時刻幅あたりの平均到着数が10回の場合の分布をみてみる。

poisson-process-uniform-dist-n010

各時刻幅で到着数が変動していて、平均10回の到着が想定されるところに3～20回とばらついている。

次にrank.width = 500すなわち時刻幅あたりの平均到着数が50回の場合の分布をみてみると、今度は随分と各時刻幅での到着回数が揃ってきているのが確認できる。

poisson-process-uniform-dist-n050

到着時間間隔の分布

到着時間間隔は、時系列上としてソートされた到着データてn番目の要素からn – 1番目の要素を引いて得られる。この計算は、到着データの先頭を除いたベクトルから最後尾のデータを除いたベクトルを引いて行っている。

設定条件は到着率λ = 1/10[回/秒。観測時間は分布形には影響を与えず、データ数にのみ関係する。今回はt = 1000秒で、データ数はλt – 1 = 99個。

また、理論上の指数分布の密度関数の曲線をヒストグラムに重ね合わせている。

lambda <- 1 / 10
t.obs <- 1000
num.data <- lambda * t.obs

org.data <- (runif(floor(num.data * 3)) * 3 - 1) * t.obs
data <- sort(org.data[(org.data >= 0) & (org.data < t.obs)])

interval <- data[-1] - data[-length(data)]

hist(interval, prob=TRUE, freq=FALSE)
curve(dexp(x,rate=lambda), add=TRUE)

lambda <- 1 / 10

t.obs <- 1000

num.data <- lambda * t.obs

org.data <- (runif(floor(num.data * 3)) * 3 - 1) * t.obs

data <- sort(org.data[(org.data >= 0) & (org.data < t.obs)])

interval <- data[-1] - data[-length(data)]

hist(interval, prob=TRUE, freq=FALSE)

curve(dexp(x,rate=lambda), add=TRUE)

この表示結果は以下のようになる。データ数は99だが、比較的よく指数分布の密度関数と整合している。

poisson-process-uniform-dist-interval100

到着数の分布

以下の手順で観測時間tの間の到着数のデータを用意する。

必要回数分、以下を繰り返す
- 観測時間をtとして、区間[–t, 2t)で乱数を発生
- 0以上、t未満のデータを除いたデータを得てソート
- その個数を数える

Rのコードは以下の通り。到着率はこれまでと同じ1/10、観測時間は100秒で、この観測時間内の期待到着数は10回だが、これがどの程度ばらつくかを確認するため、このデータを1000回発生させて、ヒストグラムを描かせている。

lambda <- 1 / 10
t.obs <- 100
num.data <- lambda * t.obs
n.loop <- 1000

arrivals <- rep(0, num.data)

for (i in 1:n.loop) {
  org.data <- (runif(floor(num.data * 3)) * 3 - 1) * t.obs
  data <- org.data[(org.data >= 0) & (org.data < t.obs)]
  arrivals[i] <- length(data)
}

hist(arrivals, prob=TRUE, xlim=c(0,20), breaks=seq(0,20,2.5))

lambda <- 1 / 10

t.obs <- 100

num.data <- lambda * t.obs

n.loop <- 1000

arrivals <- rep(0, num.data)

for (i in 1:n.loop) {

org.data <- (runif(floor(num.data * 3)) * 3 - 1) * t.obs

data <- org.data[(org.data >= 0) & (org.data < t.obs)]

arrivals[i] <- length(data)

}

hist(arrivals, prob=TRUE, xlim=c(0,20), breaks=seq(0,20,2.5))

到着率、観測時間が同じパラメータのPoisson分布もあわせて描いたのが以下のグラフ。

poisson-process-uniform-dist-poisson-n010

到着率を変えずに観測時間tを100→1000秒にした場合の分布は以下のようになる。

poisson-process-uniform-dist-poisson-n100

ここで縦軸に注意。頻度／確率のスケールが初めのものより小さくなっている。つまり、このグラフは最初のグラフに比べて、右の方へずれて、高さはぐっと低くなっている(この傾向はPoisson分布の傾向)。

以下に、パラメータλtのtをj変化させたときのヒストグラムを重ね合わせてみる。

generate <- function(dummy, lambda, t) {
  n <- lambda * t
  org.data <- (runif(floor(n * 3)) * 3 - 1) * t
  data <- org.data[(org.data >= 0) & (org.data < t)]
  return (length(data))
}

lambda <- 1 / 10
t.obs <- 100

arrivals1 <- sapply(rep(0, 1000), generate, lambda, 50)
arrivals2 <- sapply(rep(0, 1000), generate, lambda, 100)
arrivals3 <- sapply(rep(0, 1000), generate, lambda, 200)

xlm <- c(0, 35)
ylm <- c(0, 0.25)
brk <- seq(0, 35, 2.5)
hist(arrivals1, prob=TRUE, col="#ff000040", xlim=xlm, ylim=ylm, breaks=brk)
hist(arrivals2, prob=TRUE, col="#00ff0040", xlim=xlm, ylim=ylm, breaks=brk, add=TRUE)
hist(arrivals3, prob=TRUE, col="#0000ff40", xlim=xlm, ylim=ylm, breaks=brk, add=TRUE)

generate <- function(dummy, lambda, t) {

n <- lambda * t

org.data <- (runif(floor(n * 3)) * 3 - 1) * t

data <- org.data[(org.data >= 0) & (org.data < t)]

return (length(data))

}

lambda <- 1 / 10

t.obs <- 100

arrivals1 <- sapply(rep(0, 1000), generate, lambda, 50)

arrivals2 <- sapply(rep(0, 1000), generate, lambda, 100)

arrivals3 <- sapply(rep(0, 1000), generate, lambda, 200)

xlm <- c(0, 35)

ylm <- c(0, 0.25)

brk <- seq(0, 35, 2.5)

hist(arrivals1, prob=TRUE, col="#ff000040", xlim=xlm, ylim=ylm, breaks=brk)

hist(arrivals2, prob=TRUE, col="#00ff0040", xlim=xlm, ylim=ylm, breaks=brk, add=TRUE)

hist(arrivals3, prob=TRUE, col="#0000ff40", xlim=xlm, ylim=ylm, breaks=brk, add=TRUE)

poisson-process-uniform-dist-poisson-3params

Rによるポアソン過程のシミュレーション

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

概要

一様乱数を基にしたポアソン過程のシミュレーションは前にまとめたが、その時には表計算ソフトを使ったので操作性も悪く、限られたケースのみ計算したレベルだった。

その後、統計分析用のソフト”R“を少し勉強したので、より掘り下げた確認をしてみる。

前提条件と手法

平均すると10秒に1回、客が訪れる事象を考える。単位時間を1秒とすると、到着率は $\lambda = 1/10$ となる。この条件下で、以下のことを確認する。

一様乱数によったときの確率分布
- 時系列上での到着時刻の分布～一様分布
- 到着時間間隔の分布～指数分布
- 到着数の分布～Poisson分布
時刻制御(time driven)によったときの各種確率分布
- 同上
イベント制御(event driven)ときの各種確率分布
- 時系列上での到着時刻の分布～一様乱数
- 到着数の分布～Poisson分布

考え方

一様乱数による方法

下図のように、観測時間tの間で一度に一様乱数を発生させ、それを昇順にソートしたものを到着時刻の時系列データとする。

poisson-process-fig-uniform-dist

単に到着率λと観測時間tを掛けた平均到着数で乱数を発生させた場合、観測時間内の到着数が一定となってしまうため、観測時間と同じ時間幅で前後に余裕をとり、その中で一様乱数を発生させたうえで両端を切り落としている。

時刻制御とイベント制御

下図に、時刻制御(time driven)とイベント制御(event driven)の考え方の違いを比較して示している。

poisson-process-fig-time-and-event-driven

時刻制御については、観測時間を小さな時間間隔に区切って、区間ごとに到着確率に沿って到着の有無を決めながら、時刻を前に進めていく。

イベント制御の場合は、ポアソン過程で到着時間間隔が指数分布に従うことを利用して、指数乱数を発生させながら順次に到着時刻を算出していく。

R – sapplyについて

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

ベクトル演算とsapplyの違い

ベクトルの各要素に関数を作用させたいとき、普通は次のようにする。

> f <- function(x) x*x
> f(seq(1, 5))
[1]  1  4  9 16 25

> f <- function(x) x*x

> f(seq(1, 5))

[1] 1 4 9 16 25

これを、sapplyを使って次のようにも書ける。

> f <- function(x) x*x
> sapply(seq(1, 5), f)
[1]  1  4  9 16 25

> f <- function(x) x*x

> sapply(seq(1, 5), f)

[1] 1 4 9 16 25

一見すると何ら違いはないように見えるが、次のような場合に違いが出てくる。

> f <- function(x) x * runif(1)
> f(rep(1, 5))
[1] 0.7131735 0.7131735 0.7131735 0.7131735 0.7131735
> sapply(rep(1, 5), f)
[1] 0.56394702 0.71369838 0.09517609 0.72681560 0.47276722

> f <- function(x) x * runif(1)

> f(rep(1, 5))

[1] 0.7131735 0.7131735 0.7131735 0.7131735 0.7131735

> sapply(rep(1, 5), f)

[1] 0.56394702 0.71369838 0.09517609 0.72681560 0.47276722

ベクトル演算の場合は、全ての要素に対して同じ演算結果が適用されるる。このためrunif(1)が1度だけ計算され、各要素にその結果が適用されるので、全要素が同じ値になる。

sapplyの場合は、各要素に対して毎回演算が行われる、このため各要素ごとにrunif(1)が計算されて適用されるので、各要素の値が異なる。

各要素に新規の計算結果を入れたいとき

各要素の初期値と無関係に要素に計算結果を適用したいと意図したとき、関数にダミー変数を書いておいて処理する。

例えば以下の例では、5つの要素に新たな乱数を発生させてセットしようとしている。

> f <- function(dummy) runif(1)
> sapply(rep(0, 5), f)
[1] 0.2031135 0.4134661 0.9572596 0.6462498 0.2886842
> f(rep(0, 5))
[1] 0.4594526
> rep(f(0), 5)
[1] 0.7015942 0.7015942 0.7015942 0.7015942 0.7015942

> f <- function(dummy) runif(1)

> sapply(rep(0, 5), f)

[1] 0.2031135 0.4134661 0.9572596 0.6462498 0.2886842

> f(rep(0, 5))

[1] 0.4594526

> rep(f(0), 5)

[1] 0.7015942 0.7015942 0.7015942 0.7015942 0.7015942

sapplyの場合は期待通りの結果。この場合は関数に引数が要求されるので、dummyという引数を与えている。

関数の引数にベクトルを入れると、その引数に対して1回だけ乱数を発生させるので、結果は1個の数値。

rep()関数の引数にf()を適用すると、1回だけ発生した乱数がすべての要素に適用されている。

複数の引数を持つ関数の場合

次の例は、各要素それぞれに0～2の範囲の乱数をセットしている。

> f <- function(dummy, a) runif(1) * a
> sapply(rep(0, 5), f, 2)
[1] 0.8868690 1.0057694 1.3814147 1.4157728 0.2052112

> f <- function(dummy, a) runif(1) * a

> sapply(rep(0, 5), f, 2)

[1] 0.8868690 1.0057694 1.3814147 1.4157728 0.2052112

sapplyに渡す関数が2つ以上の引数を持つときは、関数に続けて残りの引数を渡せばよい。

R – ヒストグラム

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

基本形

ヒストグラムは、ベクトルでデータと階級を与えて表示させる。

data <- c(42,68,57,48,75,81,65,93,51,66,79,69,63,80,53)
hist(data)

1 2	data <- c(42,68,57,48,75,81,65,93,51,66,79,69,63,80,53) hist(data)

r-hist-ex-smp

オプション

階級値の指定と度数／頻度の選択

breaks属性に指定したベクトルで階級値を指定
ylim属性に指定したベクトルで縦軸の範囲を指定
prob=TRUEを指定することで縦軸を度数から頻度に変更

data <- c(42,68,57,48,75,81,65,93,51,66,79,69,63,80,53)
ranks <- seq(40,100,15)
density <- c(0, 0.03)
hist(data, breaks=ranks, ylim=density, prob=TRUE)

data <- c(42,68,57,48,75,81,65,93,51,66,79,69,63,80,53)

ranks <- seq(40,100,15)

density <- c(0, 0.03)

hist(data, breaks=ranks, ylim=density, prob=TRUE)

r-hist-ex-smp-axis

ラベルや色の指定

main属性とxlab属性に文字列を指定して、表題と横軸のラベルを設定
col属性でヒストグラムの色指定が可能

data <- c(42,68,57,48,75,81,65,93,51,66,79,69,63,80,53)
hist(data, main="Data", xlab="Range", col="#008888")

1 2	data <- c(42,68,57,48,75,81,65,93,51,66,79,69,63,80,53) hist(data, main="Data", xlab="Range", col="#008888")

r-hist-ex-smp-lblcol

複数グラフ

横に並べる

data1 <- c(42,68,57,48,75,81,65,93,51)
data2 <- c(66,79,69,63,80,53)

par(mfrow=c(1,2))

ranks <- seq(40,100,15)
hist(data1, breaks=ranks)
hist(data2, breaks=ranks)

data1 <- c(42,68,57,48,75,81,65,93,51)

data2 <- c(66,79,69,63,80,53)

par(mfrow=c(1,2))

ranks <- seq(40,100,15)

hist(data1, breaks=ranks)

hist(data2, breaks=ranks)

r-hist-ex-col2

縦に並べる

data1 <- c(42,68,57,48,75,81,65,93,51)
data2 <- c(66,79,69,63,80,53)

par(mfrow=c(2,1))

ranks <- seq(40,100,15)
hist(data1, breaks=ranks)
hist(data2, breaks=ranks)

data1 <- c(42,68,57,48,75,81,65,93,51)

data2 <- c(66,79,69,63,80,53)

par(mfrow=c(2,1))

ranks <- seq(40,100,15)

hist(data1, breaks=ranks)

hist(data2, breaks=ranks)

r-hist-ex-row2

関数曲線との重ね合わせ

ヒストグラムの定義域の数値をそのまま使い、関数を渡して描画させることができる。

data <- runif(100)
hist(data, freq=FALSE)
curve(x*2, add=TRUE, col="blue", lwd=2)
curve(-log(x), add=TRUE, col="red", lwd=3)

data <- runif(100)

hist(data, freq=FALSE)

curve(x*2, add=TRUE, col="blue", lwd=2)

curve(-log(x), add=TRUE, col="red", lwd=3)

要点は以下の通り

hist()関数でfreq=FALSEを指定
curve()関数の第一引数で描画したい関数を指定(独立変数をxとする)
curve()関数でadd=TRUEを指定

r-hist-ex-superimpose

これを利用して、ヒストグラムに確率密度関数を重ねて描画することができる。

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の正規分布に従う乱数を生成する。

一様乱数に基づく分析

到着時刻の分布

到着時間間隔の分布

到着数の分布

関連リンク

概要

前提条件と手法

考え方

一様乱数による方法

時刻制御とイベント制御

関連リンク

ベクトル演算とsapplyの違い

各要素に新規の計算結果を入れたいとき

複数の引数を持つ関数の場合

基本形

オプション

階級値の指定と度数／頻度の選択

ラベルや色の指定

複数グラフ

横に並べる

縦に並べる

関数曲線との重ね合わせ

概要

基本形

タイプ指定

plot()のタイプ

折れ線グラフ

関数を直接与える場合

sort()関数

rank()関数

数値演算

ベクトル同士の数値演算

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

比較演算

数学関数

集合演算

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

単一要素の取得

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

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

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

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

要素・ベクトルの置換

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

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

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

1ずつ増減

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

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

sequence()関数

乱数列の生成

一様分布

二項分布

指数分布

ポアソン分布

正規分布