R – rbindについて

rbind()で困った点

データフレームにrbindでアイテムを追加する場合、紹介されている例ではみな、データをベクトルで与えている(注:文字列型が勝手にFactor型にされてしまうのを4行目で回避している)。

一見するとうまくいっているが、実は問題がある。ageの項目が数字のはずだと内容を参照したり計算しようとすると・・・

となってしまう。

これはRのベクトルの仕様で、1つのベクトルの要素の型がすべて同じでなければならず、文字列型と数値型が混在した場合は数値が文字列に変換されてしまうため。

これを回避するために幾つかの方法があるが、たとえばここではas.numeric()を用いてみる。

うむ、うまくいっている、と思って、これに新たなデータを追加すると・・・

またしても怒られた。これはデータ追加のたびにベクトルの段階で文字列への変換が行われてしまうためらしく、上の9~10行目で前に数値型にしたはずのデータまで文字列に変換されてしまっている。

ベクトルではなくリストを使う

そこで、rbindでデータ群を与えるのに、ベクトルではなくリストを使ってみた。

見た目はベクトルの時と同じく、うまくいっているように見える。そこで各要素をチェックしてみると・・・

ちゃんと数値として扱われている。それでは新たにデータを追加しても大丈夫か。

問題なし。

よって、文字列と数値が混在する場合は、rbindの引数にはリストを使うべき。

 

R – 因子型(Factor型)

準備

以下のように文字列型のベクトルを準備。

因子型への変換

factor()関数を使って因子型へ返還。このとき因子の順番はシステムで自動的に割り振られ、ソートするとその順番で並べ替えられる。

因子を指定した要素の抽出

要素の抽出は因子のラベルを指定して行う。
因子のオーダでの指定はできない。

因子の順序の指定

factor()関数のlevels指定で明示的に因子の順序を指定できる。

因子の大小関係の指定

上記のベクトルyやzの各因子は順序関係を持っているが、それらは順序については特定できるが大小判定は行えない。

大小判定可能な値とするにはordered()関数を通す必要がある。

Factorベクトルの新規生成

一つ目の要素を定義する場合

通常紹介されている方法は、文字列型などのベクトルがあらかじめ準備されていて、それをFactor型に変換するというもの。

一つ目の要素を定義して、以降付け足していきたい場合は以下のようにする。

またlevelsを指定しないと、初期値以外はNAとなって警告が出る。

またlevelsで指定した以外のラベルを指定するとNAとなり、最初の要素は無視され、その後追加しようとした要素はNAとして保存される。

要素数ゼロから定義する場合

最初からFactor型のベクトルを定義してゼロから要素を追加するには、以下のようにするとよい。

下記の例では、factor()関数の第1引数に仮のベクトルをnumeric(0)で与えているが、character(0)でも同じ結果となり、単なるプレースホルダとなっているらしい。

データフレーム生成の場合の因子化について

データフレーム生成時に、文字列型ベクトルが自動的に因子型に変換される。これを抑止する方法についてはこちらを参照。

 

R – 画面表示

オブジェクトの直接表示

オブジェクトを直接入力すると、その内容が表示される。

スクリプトファイルでも行頭にオブジェクトを書くと、その行の実行時点でオブジェクトの名用が表示される。

ただし条件判断やループの場合は、文でもブロックでも表示はされない。

例えば以下のスクリプトの場合、

表示されるのは1行目だけ。

print()~一般的な表示

オブジェクトの内容をそのまま表示する。

デフォルトでは文字列が””で囲われるが、quote=Fを指定すると””が表示されなくなる。

page()~別ウィンドウでの表示

page()関数ごとに別ウィンドウが立ち上がり、オブジェクトが表示される。大量のデータを表示するときに便利。

デフォルトではmethod="dput"が指定され、オブジェクトの定義表現が表示される。

オブジェクトの内容を表示するときは、明示的にmethod="print"を指定する。

1つ目のpage()関数で別ウィンドウが立ち上がり、以下のように表示される。

2つ目のpage()関数で2つ目の別ウィンドウが立ち上がり、以下のように表示される。

cat()~文字列の表示

cat()関数は、括弧内の文字列をそのまま出力する。

catは最後に改行をつけないので、改行させたい時は”\n”を付ける。

 

str()~オブジェクト情報付き表示

要約されたオブジェクトの情報を付けて内容を表示する。

summary()~データの要約を表示する。

オブジェクトの内容に応じた要約情報を表示する。

summaryは特に統計データの要約に重要。たとえば▲データフレームの集計▲を参照。

 

 

R – データフレームの参照・変更

データフレーム全体の参照

データフレームの名前そのもので、データフレームの内容を参照できる。列名(name, age)や項目名(“1″~”3”)も含めて参照される。

データフレームの行数・列数の参照

nrow()関数、ncol()関数でデータフレームの行数、列数を参照できる。dim()関数は(行数, 列数)のベクトルを返す。

データフレームの内容の参照・変更

列名・項目名の参照・変更

列名はnames()関数かcolnames()関数で、項目名はrownames()関数で取得できる。

列名や項目名は、それぞれの参照関数にベクトルを代入することで変更できる。

データフレームの要素の参照・変更

行・列を直接指定して参照する場合、列の場合は列名を指定するか列番号で、行の場合は行番号を指定して参照する。行・列の番号を指定して1つのデータを取得することもできる。

行番号を指定する場合は後ろに、列番号を指定する場合は前に”,”をつける必要があり、これは後述の行の抽出の場合に重要になる。

なお、参照した要素の右に代入文を書くことで、その要素や業・列の内容を変更できる。

行・列の追加と削除

行・列の追加

行の追加はrbind()関数で、列の追加はcbind()関数で行う。

これらの関数は元のデータフレームを変更せず、新たなデータフレームを結果として返す。

rbindの注意点として、文字型を意図した項目はあらかじめFactor型から文字型に変更しておかないと、文字列を含んだデータを結合しようとするとエラーになる。

cbindの注意点としては、デフォルトでは引数の変数名が項目名にあてられる。直接c(…)と書いたりすると、それがそのまま項目名になってしまう。

【追記】 以下の例ではrbindの引数にベクトルを渡しているが、文字列と数値が混在している下記のような例では、これは危ない→rbindについてを参照

行・列の削除

行や列の番号を指定して削除する場合は、番号にマイナスをつける。範囲指定も可能。

行・列の追加と同じく、これらの操作も元のデータフレームを変更せず、新たなデータフレームを結果として返す。

データの抽出

列項目の条件を指定して、要素を抽出することができる。

指定の条件の後に”,”をつけるのを忘れないこと。これは行要素に対して条件指定していることを表している。

抽出操作はもとのデータフレームに影響を与えず、結果は新たなデータフレームとして返される。

 

R – データフレームのソート

order()関数

以下のデータを準備する。

このデータフレームをheightの項目の昇順のオーダーでソートしたいとする。

order()関数というのがあって、引数の項目でソートしたときに元データの項目番号がどういう順番で並ぶかを返してくれる。

たとえば1番目のデータは1番小さい157で、これは5番目のデータ。2番目のデータは2番目に小さい160で2番目のデータ・・・という風に、heightデータを昇順に並べた時の元データの項目番号の並びを教えてくれる。

rank()関数がデータを昇順に並べた時の、そのデータのソート後の順番を返すのとは違う点に注意。

この結果を元のデータフレームの行部分に使うと、その順番で並べ替えられたデータを返してくれる。

降順ソート

データを降順にソートしたい時は、decreasing=Tで指示する。

複数列のソートなど

以下の例は、数値以外でもソートできることと、複数のオーダーを組み合わせたソートの例を示している。

行番号の変更

ソート後に行ラベルを振り直したい場合は、rownames()の内容を変更する。

 

 

R – データフレームの集計

summary()関数で、データ全体の集計が行える。

by()関数によって、特定のデータの層別に集計することができる。

もちろん、個別の数値データについて各種代表値を計算することもできる。

 

 

R – データフレームの列の非因子化

ベクトルから生成する場合

データフレームを複数のベクトルから生成する場合、文字列のベクトルがFactor型(因子型)になってしまう。

これを文字型に直すのに、以下の3つの方法がある。

as.charactor()関数

as.charactor()関数で指定した列がcharacter化されるので、それを元の列に代入。

levels()関数

仕組みはよくわかっていない。

transform()関数

データフレーム全体に適用。

ファイルから読み込む場合

ファイルから読み込んだ後に上記の操作を行ってもよいが、読み込み時にFactor型への変換を抑止するため、引数にstringsAsFactors = Fを指定する方法もある。

個の引数指定は、read.table、read.delimread.csvの何れでも指定できる。

 

 

R – データフレームの生成

ベクトル列から生成する方法

同じ長さのベクトルを各列として、データフレームを生成可能。

上記の結果は以下の通り。

ファイルからの読み込み

通常、作業ディレクトリかその下にあるディレクトリにあるファイルを扱う。扱いやすいのはタブ区切りテキスト、CSV。

read.tableによる方法

TAB区切り/ヘッダなし

read.tableでそのまま読み込む。項目名はRが勝手につけてくれる。

デフォルトでTABが区切り文字になってるが、明示する場合は引数にsep="\t"を追加。

CSV/ヘッダなし

read.tableで引数にsep=","を付けて読み込む。

TAB区切り/ヘッダあり

read.tableheader=Tを付けて読み込む。sep="\t"を指定してもよい。

この形式の場合、wrapperクラスのread.delimでファイルのみ指定して読み込める。

CSV/ヘッダあり

read.tablesep=","header=Tを指定して読み込む。

この形式の場合、wrapperクラスのread.csvでファイルのみ指定して読み込める。

注意点:Factor型への自動変換

ベクトルからの生成にしても、ファイルからの読み込みにしても、文字列のデータを持ったベクトルをデータフレームに取り込むと、自動的にFactor型(因子型)にされてしまう。

これを解決するための方法は列の非因子化を参照。

 

イベント制御によるポアソン過程のシミュレーション

イベント制御による分析

時刻制御(time driven)では微小な時間間隔の間において到着イベントを確率的に発生させた。そこでは結果としての到着時間間隔の指数分布や到着数のPoisson分布は前提にはおかれていない。

イベント制御(event driven)の考え方では、ある到着があった後、次の到着までの時間間隔を指数分布に従う乱数を発生させて計算する。計算にはR言語を使う。

イベント制御の考え方の概要やその他の考え方についてはRによるポアソン過程のシミュレーションを参照。

到着イベントの時系列

ここでは、到着率はこれまでと同じλ = 1/10として、時刻t = 0からt ≤ 1000秒の間、到着の時間間隔が確率密度e^{- \lambda t}の指数分布に従うような乱数を発生させていく。

指数分布に従う乱数を生成するには逆関数法を使う。

指数分布の確率分布関数は以下の通り。

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

一様乱数runifを与えて指数分布に従う乱数を得るには、一様乱数を分布関数の逆関数に適用すればよい。

(2)    \begin{equation*} F^{-1}(t) = - \frac{1}{\lambda}\ ln(1 - {\rm runif(1)}) \end{equation*}

時刻0を起点として、1番目の到着時刻をこの指数乱数で決定し、その次の到着時刻も指数乱数で・・・と繰り返して、時刻が上限を超えるまで時系列として記録する。

一様乱数による場合と同じく、到着率1/10、観測時間5000秒として、階級幅100秒、500秒の場合の時系列を生成したのが下図。階級あたりのデータ数と分布の平滑さの関係は、一様乱数による場合やtime drivenの場合と同じ。

poisson-process-event-driven-arrival100

poisson-process-event-driven-arrival500

到着数の分布

到着数の確率分布は、一様乱数による場合やtime drivenと同じように集計する。

結果は、Poisson分布の理論式ともよくあっている。

poisson-process-event-driven-poisson

関連リンク

 

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

一様乱数に基づく分析

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 到着率λ
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個。

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

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

poisson-process-uniform-dist-interval100

到着数の分布

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

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

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

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

poisson-process-uniform-dist-poisson-n010

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

poisson-process-uniform-dist-poisson-n100

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

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

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

関連リンク