R 備忘録 2013/05/10-2

ちなみに、「意味がわかる 統計解析」の流れに沿って、書いてあることが R で実行できるように R の機能を勉強しています。

確率分布

R には確率分布 DIST ごとに（DICT の累積分布関数を F(x)、確率密度関数を f(x) とおくと）

• dDIST(x)：f(x) の値を返す。
• pDIST(q)：F(x) の値を返す。
• qDIST(p)：F(x) の値が p となる x を求める。
• rDIST(n)：DIST に従って乱数を n 個生成し、長さ n のベクトルとして返す。

という４つの関数があります。
DIST にあたる部分には norm（正規分布）、unif（一様分布）、pois（ポアソン分布）などがあり、パラメータは引数で調節できます。
例えば

> dunif(0.0, -1, 1) # [-1:1] の一様分布
[1] 0.5
> pnorm(0.0)
[1] 0.5
> qnorm(0.5)
[1] 0
> runif(5)
[1] 0.58321736 0.05653621 0.51690395 0.63476821 0.69568012

となります。

乱数種（シード）の設定

set.seed を使います。

> set.seed(1)
> rnorm(3)
[1] -0.6264538 0.1836433 -0.8356286
> rnorm(3)
[1] 1.5952808 0.3295078 -0.8204684
> set.seed(1)
> rnorm(3)
[1] -0.6264538 0.1836433 -0.8356286

ちなみにベクトルを与えたら第一要素がシードになりました。

> set.seed(1:5)
> rnorm(3)
[1] -0.6264538 0.1836433 -0.8356286
> set.seed(5)
> rnorm(3)
[1] -0.8408555 1.3843593 -1.2554919
> set.seed(1)
> rnorm(3)
[1] -0.6264538 0.1836433 -0.8356286

基本統計量

mean, var sd などその名（省略名ですが）の通りの関数を使えば基本統計量を求めることができます。
なお var, sd は不偏分散、不偏標準偏差です。

> vec <- rnorm(1000)
> mean(vec)
[1] -0.02148036
> var(vec)
[1] 1.097021
> sd(vec)
[1] 1.047388
> median(vec)
[1] -0.03616252
> max(vec)
[1] 3.810277
> min(vec)
[1] -3.008049
> quantile(vec)
0% 25% 50% 75% 100%

• 3.00804860 -0.73133329 -0.03616252 0.69697293 3.81027668

やっかいなのは度数分布の基本統計量です。

> vec_freq <- table(cut(vec, breaks=seq(-5, 5, 0.1)))

今は個表データから度数分布を生成しているので本当は元の個表データから基本統計量を計算した方が単純なのですが、そういうのは無視して、上記のような度数分布があった場合、

> weight <- seq(-5, 4.9, 0.1) + 0.05
> weighted.mean(weight, vec_freq)
[1] -0.0232

というように階級値のベクトルを作ってあげれば、weighted.mean で度数分布の平均を求めることができます。
が、そこまでするなら自分で計算してあげた方が早い気もします。

> sum(weight * vec_freq) / sum(vec_freq)
[1] -0.0232
> weight %*% vec_freq / sum(vec_freq)
[,1]
[1,] -0.0232

なお、分散等を求めるための組み込み関数はありません。
重みを考慮しながら自分で計算するか、データ数がそれほど大きくない場合には扱いを簡単にするために仮想的に個票データを作ってしまうのが良いと思います。

> (vir_vec <- do.call(c, mapply(function(f, w) return (rep(w, times=f)), as.vector(vec_freq), weight)))
[1] -3.05 -2.95 -2.95 -2.85 -2.55 -2.55 -2.55 -2.55 -2.45 -2.45 -2.45 -2.45
[13] -2.45 -2.45 -2.35 -2.35 -2.35 -2.35 -2.25 -2.25 -2.25 -2.25 -2.25 -2.25
[25] -2.15 -2.15 -2.15 -2.05 -2.05 -2.05 -2.05 -2.05 -2.05 -2.05 -1.95 -1.95
[37] -1.95 -1.95 -1.95 -1.95 -1.95 -1.95 -1.85 -1.85 -1.85 -1.85 -1.85 -1.85
[49] -1.75 -1.75 -1.75 -1.75 -1.75 -1.75 -1.75 -1.65 -1.65 -1.65 -1.65 -1.65
# 長いので省略
> mean(vir_vec)
[1] -0.0232
> var(vir_vec)
[1] 1.097319

do.call はベクトルを関数の引数に展開するための関数で、例えば do.call(f, c(1, 2, 3)) は f(1, 2, 3) と同じ意味になります。

基本統計量の要約

summary 関数を使えば、基本統計量をいくつかまとめて出力することができます。

> summary(vir_vec)
Min. 1st Qu. Median Mean 3rd Qu. Max.

• 3.0500 -0.7500 -0.0500 -0.0232 0.6750 3.8500

summary も引数によって挙動がいろいろと変わるので総称関数っぽいですね。