データ解析言語Rによる統計的プログラミング: 第 2 回 機能プログラミングとデータ操作

異常の検出と解析

Rは無料ソフトとしてリリースされている、機能豊富な統計環境です。この記事ではDavidとBradが3回シリーズの第2回として、このRの紹介を続けます。前回の記事でデータが整ったので、今回はこの言語の機能に入り込みます。

David Mertz, Ph.D (mertz@gnosis.cx), Author, Gnosis Software, Inc.

Photo of David MertzDavid Mertz氏は多くの分野で活躍しています。ソフトウェア開発や、それについて著述もしています。その他、学術政策理念について分野を問わず、関係する雑誌に記事も書いています。かなり以前には、超限集合論、ロジック、モデル理論などを研究していました。その後、労働組合組織者として活動していました。そして、David Mertz氏自身は人生の半ばにもまだ達していないと思っているので、これから何かほかの仕事をするかもしれません。



Brad Huntting (huntting@glarp.com)University of Colorado

Bradは過去14年間、それぞれ別の3社で、UNIXシステム管理とネットワーク・エンジニアリングに従事してきました。現在はボールダーにあるコロラド大学で応用数学の博士過程に取り組んでおり、コンピューター・サイエンス学部のUNIXサポートをすることで学費の支払いに充てています。



2004年 10月 07日

Rは非常に機能的なプログラミング言語であると同時に、データ・セットの統計分析用の汎用環境でもあります。環境としてRを使うことにより、カスタムのコマンドラインからデータを様々にグラフィック表現することができます。このシリーズ第1回の記事とCameron LairdによるRの記事(どちらも参考文献にあります)では、Rとその背景、またRが実行するプラットフォームや、R用として入手できるパッケージに関する情報などを説明していますので、併せてご覧ください。

PythonやRuby、Tcl/Tk用のwish、その他多くのLisp環境の対話型シェルと同じように、Rシェルは様々な操作を行う手段として非常に有効なものであり、コマンドの再呼び出しやコマンド編集を使って、コマンドの変種を試すことができます。他の言語での対話型シェルとは対照的に、ただしデータ中心的というRの性格に合わせて、Rシェルでは完全な環境(ワーキング・ディレクトリー毎に一つ)をオプションとして保存します。この中でコマンド・ヒストリーは便利な部分で、.Rhistoryファイルの中で検証したり修正したりすることができます。しかし、環境を保存する場合に最も重要なのは、バイナリー形式で .RDataに保存される動作データです。ちなみにrm()コマンドを注意して使うと、保存されるデータ環境が無制限に大きくなってしまうのを防ぐことができます(remove()rm()と同意語です)。

この記事はRに関するシリーズの第2回目です。第1回ではRのデータ・タイプの基本として、ベクターから始め、多次元配列(二次元配列は「マトリックス」としても知られています)、「スマート」テーブルのデータ・フレーム、異種集合のリストなどを紹介しました。前回の記事では、この記事の共著者であるBradが持っていた大きなデータ・セット(Bradの家の内外で3分毎に測定された、一年間に渡る温度変化)に対して、基本的な統計解析やグラフ化も行いました。実世界での数多くのデータ・セットと同じく、Bradの温度データには欠落や異常が含まれていることが分かっていますし、分かっていなくても予想ができます。そこで第2回目の今回は、データの中にあるおかしな部分を分析してみます。

この記事の概略として、次の2点がポイントです。まず、前回使ったものと同じデータを使って、R言語自体を前回よりも詳細に探ります。また、データの中にある一般的なパターンを調べ、そのパターンをRでどのように見つけ出すかを説明します。

プログラミング言語としてのR

(GPLとなっている)Rプログラミング言語には、S/S-PLUSとScheme (Lisp) という2つの親があります。Rの持つ機能プログラミングという側面はScheme側から来ているので、ここではSchemeとの関連を強調する必要があります。注意深い読者であれば、第1回の記事の例では不思議なことに、フロー制御が抜けていることに気がついたでしょう。この抜けは今回の記事でも続きます。実を言うと、Rにはきちんとifや、elsewhileそれにforなどのコマンドがあり、どれもPythonでの同じ綴りのコマンドとほとんど同じなのです(他の言語では同じコマンドを少し違った風に綴ります)。適切な冗長性を持たせるためには、repeatや、breaknextswitchなどを投げ込んでください。命令型言語(imperative languages)のプログラマーであれば、フロー制御が無くても様々なことができるのに驚くでしょうし、多くのタスクではむしろフロー制御が無い方が簡単なことに気が付くと、さらに驚くことでしょう。

意図するところを宣言する

機能プログラミングの伝統そのままに、Rでは単純な宣言ステートメントを使って、ほとんど何でも行うことができます。Rの持つ2つの機能のおかげで、大部分の場合には命令型フロー制御が不要になります。第一に、集合オブジェクトに対する操作の大部分は要素単位(elementwise)で動作することを既に見てきました。ベクターの全要素に対して何かを行うのは簡単なので、データのベクターを手動でループする必要はありません。

リスト1. Rでの要素単位操作
> a = 1:10            # Range of numbers 1 through 10
> b = a+5             # Add 5 to each element
> b                   # Show b vector
 [1]  6  7  8  9 10 11 12 13 14 15

あるいは「インデックス配列」を使うことで、あるインデックスを持つ要素に対してのみ、選択的に操作することもできます。

リスト2. インデックス配列を使って要素を選択する
> c = b               # Make a (lazy) copy of b
> pos = c(1,2,3,5,7)  # Change the prime indices
> c[pos] = c[pos]*10  # Reassign to indexed positions
> c                   # Show c vector
 [1]  60  70  80   9 100  11 120  13  14  15

あるいは恐らく一番良いのは、HaskellやPythonでのリスト内包とよく似た構文を使い、必要なプロパティを持つ要素に対してのみ操作を行うようにすることでしょう。

リスト3. 述部を使って要素を選択する
> d = c
> d[d %% 2 == 0] = -1 # Reassign -1 to all even elements
> d
 [1] -1 -1 -1  9 -1 11 -1 13 -1 15

非常に鋭い読者であれば、これらの例では= を使っており、前回の記事では<- を使っていたことに気が付いているでしょう。等号は割り当て矢印と同じことをしますが、最上位レベルでのスコープでのみ使用され、ネストした表現の中では使われません。よく分からなければ、矢印の方が安全です。

温度データを探る

第1回の記事では、4つの測定点それぞれにおける最初の約17万点の読み取り値を、ファンクションread.table()を使ってデータ・フレームに読み込みました。ただし、個々の温度シリーズに対してアクセスしやすいように、測定場所の名前(outside,livingroom,basement,lab: それぞれ戸外、リビングルーム、地下室、作業室)をつけたベクターに一連のデータをコピーしています。各温度シーケンスでは幾つかのデータ点が欠けていることも忘れないでください。この原因としては、データを記録しているコンピューターがダウンしていたとか、測定器が故障していた、4つの温度計全てが全く同じ日にオンラインになるわけではないこと、などが考えられます。つまり私達が使ったデータは、皆さんが扱うであろう実世界のデータとよく似た性質のものだと言えるのです。

温度データのセットを調べた前回では、データをヒストグラムにプロットしてみると異常があることに気がつきました。戸外温度の摂氏24度のごく近くに、大きな突出があるようです。この原因として私達は最初、時々温度読み取り値の伝送順序が変わってしまい、室内温度の一つが戸外温度のログに記録されてしまう(そして、これに対応して戸外温度が室内温度ログの一つとして記録されてしまう)のだと考えました。Rを試す一つの手段として、この説明が正しいかどうかを見てみることにしましょう。

異常を検出する

異常なデータを探る第一歩は、とにかくそれを見つけることです。具体的に言うと、単独の点の異常は、その前後のデータ点に対する突然の温度変化として分かります。さらに具体的に言うと、異常点らしき点の前後の点が、その異常点らしき点よりも大幅に高いか低いであろうと予測できます。例えば、(3分間隔での)温度の単純シーケンス「20, 16, 13」は異常に早い温度低下を表していますが、真ん中の読み取り値が単一のエラーだとは示していません。もちろん、だからといって私達は単なる孤立データ以上に関わるような、他のタイプのエラーや異常の存在を事前排除するわけではありません。

おかしな読み取り値を特定しようと私達が最初に考えた方法は、複雑なものでした。つまりシーケンスに対するフーリエ変換における高周波成分を見ればよいだろう、とか、温度ベクターに対して導関数を適用して解析する(多分、2次導関数で変化点が得られる)、さらには温度ベクターに対して畳み込み演算を行ってみる、などです。これらはどれもRに組み込まれています。しかしこうした高度な案をあれこれ考えた挙げ句、私達は足を地に着けることにしました。私達が探している単純なパターンというのは、前後の値と比べて、大きくかけ離れた値です。つまり、引き算の魔術を使えばよいのです。

ここでは隣接データ点との差を検出するために、ある行における列がそれぞれ、前のデータ点、現在のデータ点、後のデータ点に対応するようなデータ・フレームを作ります。これによってデータ・フレーム全体に対してフィルターをかけ、関心対象の行を選択できるようになります。

リスト4. 戸外気温での単一点の異常を検出する
> len <- length(outside)    # Shorthand name for scalar length of vector
> i <- 1:(len-2)            # Shorthand vector of data frame rows nums
> # Create a data frame for window of measurements per row
> odf <- data.frame(lst=outside[1:(len-2)],
+                   now=outside[2:(len-1)],
+                   nxt=outside[3:len] )
> # Create vector of local temperature fluctuations, add to data frame
> odf$flux <- (odf[i,"now"]*2) - (odf[i,"lst"] + odf[i,"nxt"])
> odf2 <- odf[!is.na(odf$flux),]  # Exclude "Not Available" readings
> oddities <- odf2[abs(odf2$flux) > 6,] # It's odd if flux is over 6

さて、フィルターをかけた後に何かあるのかを見てみましょう。

リスト5. 戸外温度での「異常」に注目する
> oddities
        lst  now  nxt flux
2866   21.3 15.0 14.9 -6.2
79501  -1.5 -6.2 -4.1 -6.8
117050 21.2 24.6 21.6  6.4
117059 20.6 23.4 20.1  6.1
127669 24.1 21.2 24.7 -6.4
127670 21.2 24.7 21.5  6.7

前後の点と大きく変化している読み取り値が幾つかあることが分かります。ただしその大部分は、順序を変更された読み取り値ではないようです。例えば時間ステップ79501で、温度6.2の前後は明らかにもっと高い温度です。しかし、これらのどれも、室内温度ではなさそうであり、特別寒い風が吹いたと考えた方が妥当なようです。

ただやはり、順序入れ換えが、時間ステップ117059または12670のあたりで起きたように思えます。真ん中の温度は正に室内温度の24に近く、その前後の読み取り値は氷点下ではありませんが、前後の値よりは明らかに低くなっています。どうやら春の順序入れ換えが起きているようです。

便利な操作をファンクションにラップする

ここで知りたいことは、順序入れ換えらしき幾つかの読み取り値が、他の読み取り値にも現れているかどうかです。室内測定点の測定データの中に「欠落している」戸外温度が紛れ込んでいるのが分かれば、私達の仮定が強く支持されることになります。しかし、測定点名outsidelabで置き換えるためだけに全ステップを打ち直すことはしたくありません。本来は全ステップのシーケンスを、再利用可能なファンクションの形でパラメーター化すべきだったのです。

リスト6. flux(fluctuation: 変動)検出をファンクションとしてパラメーター化する
oddities <- function(temps, flux) {
  len <- length(temps)
  i <- 1:(len-2)
  df <- data.frame(lst=temps[1:(len-2)],
                   now=temps[2:(len-1)],
                   nxt=temps[3:len])
  df$flux <- (df[i,"now"]*2) - (df[i,"lst"] + df[i,"nxt"])
  df2 <- df[!is.na(df$flux),]
  oddities <- df2[abs(df2$flux) > flux,]
  return(oddities)
}

一旦ファンクションができれば、他のデータ・セットに対してフィルターをかけたり、変動のしきい値にフィルターをかけたりすることがずっと容易になります。oddities(lab,6)を実行しても、戸外温度での順序入れ換えらしきものに一致する時間ステップは何も得られません(livingroombasementでも同様です)。ところがlab(作業室)での温度からは、ちょっと驚くような結果が出てくるのです。

リスト7. Labでの大きな温度変化
> oddities(lab, 30)
       lst  now  nxt  flux
47063 19.9 -2.6 19.5 -44.6
47268 17.7 -2.6 18.2 -41.1
84847 17.1 -0.1 17.0 -34.3
86282 14.9 -1.0 14.8 -31.7
93307 14.2 -6.4 14.1 -41.1

これは私達の予想していた変化とは違います。せいぜい説明できるとしたら、Bradが冬の真っ只中に作業室の窓を開けたのでしょう。もしそうだとしたら、3分以内に元の温度に回復してしまうという、Bradの家の暖房装置にDavidは驚嘆するばかりです。

本当の説明

おかしな時間ステップのデータで、隠れている全データを見れば何かが明らかになるでしょう。

リスト8. ステップ47063付近の全温度セット
> glarp[47062:47066,]
             timestamp basement  lab livingroom outside
47062 2003-10-31T17:07     21.5 20.3       21.8    -2.8
47063 2003-10-31T17:10     21.3 19.9       21.2    -2.7
47064 2003-10-31T17:13     20.9 -2.6       20.9    -2.6
47065 2003-10-31T17:16     20.8 19.5       20.8    -2.6
47066 2003-10-31T17:19     20.5 19.4       20.7    -2.8

まず、oddities()ファンクションが戻す時間ステップには、一ステップずれている、というエラーがあることに気がつきます。そう言えば、各データ・フレームの行に対するウィンドウを取得するためにオフセットを使ったのでした。しかし全データ自体は「ドアが開いていた」という説を裏付けるようです。コロラドでハロウィーンの日(10月31日)の夕方5時であれば、かなり寒いと考えられ、トリック・オア・トリートの子供達(trick-or-treaters)が実際Bradの家の玄関口に現れたのかも知れません(そして3分間に渡ってキャンディーを受け取っていたのかも知れません)。ですから、この秘密は解決されたのかも知れません。

しかしそれにしても、そもそも私達が調べ始める原因となった秘密はどうなのでしょう。なぜ戸外温度の24度付近に突出した点があるのでしょう。どうやらヒストグラムをもう少し詳しく見るべきなのかも知れません。叙述的判断基準(predicative criteria)を使ってベクターをインデックス化し、関心のある温度範囲のみにヒストグラムを狭めてみましょう。

リスト9. 温度ヒストグラムを狭める
hist(outside[outside < 26 & outside > 23],
     breaks=90, col="green" border="blue")
図1. 狭い温度範囲における戸外温度のヒストグラム
Outside histogram for narrow temperature range

突出した温度付近を拡大して見ると、突出のすぐ前後に、明らかな落ち込みがあるのが分かります。24.7度のすぐ近くでデータを平滑化すれば、ほとんど想定の結果にたどり着けそうに見えます。ここで私達は、温度測定装置内部での、おかしなデータ丸め動作を思い出すのです。ただし、恐らく、であって、まだ確かではありません。また、もう少し平滑化したとしても、やはり少し飛び出しはあるのです。


中間的統計解析

Rの強みの一つは、線形回帰モデルも非線形回帰モデルも計算できることです。簡単な例を見るため、2つのベクターを作りましょう。xはデータ・セットの開始から測った、一日の中でのサンプル時間であり、yはそれに対応する戸外温度です。

リスト10. 回帰ベクターを作る
> y <- glarp$outside
> x <- 1:length(y)/(24*60/3)

次で、このデータに直線を当てはめることができます。

リスト11. 連続データに線形回帰を適用する
> l1 = lm(y ~ x)
> summary(l1)
Call:
lm(formula = y ~ x)
Residuals:
     Min       1Q   Median       3Q      Max
-29.4402  -7.4330   0.2871   7.4971  23.1355
Coefficients:
        Estimate Std. Error t value Pr(>|t|)
(Intercept) 10.2511014  0.0447748  228.95   <2e-16 ***
x           -0.0037324  0.0002172  -17.19   <2e-16 ***
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 9.236 on 169489 degrees of freedom
Multiple R-Squared: 0.00174,    Adjusted R-squared: 0.001734
F-statistic: 295.4 on 1 and 169489 DF,  p-value: < 2.2e-16

「~」構文は公式オブジェクトを表します。これは実質的にRに対して、sum((y[i]-(A*x[i]+B))^2)を最小にする係数、AとBを見つけるように要求しています。最適なのはAが-0.0037324(傾斜がほとんど平坦に近い)で、Bが10.2511014の場合です。残留標準誤差は9.236であり、これは元々、yでの標準偏差とほとんど同じことに注意してください。これから、単純な時間の線形関数は戸外温度に対するモデルとして非常に不適切である、と言うことができます。

より良いモデルとしては、1日や1年間の周期を持つ、正弦、余弦関数を使うものです。このためには公式を次のように変更します。

リスト12. 三角関数曲線に合わせてみる
> l2 = lm(y ~ +I(sin(2*pi*x/365)) +I(cos(2*pi*x/365))
+             +I(sin(2*pi*x)) +I(cos(2*pi*x)) )

この公式の構文はちょっと間違いやすいものです。I()コール内部での数学記号は、その記号が通常表す意味を持っています。ですから例えば最初のI()では、xのpi倍の2倍を365で割ったもの、という意味になります。I()コールの前の「+」記号は加算ではなく、に続く因子はモデルに含まれる必要がある、ということを意味します。この結果はいつものように、summary()コマンドで表示することができます。

リスト13. 三角関数回帰での結果
> summary(l2)
Call:
lm(formula = y ~ +I(sin(2 * pi * x/365)) +I(cos(2 * pi * x/365))
                 +I(sin(2 * pi * x)) +I(cos(2 * pi * x)))
Residuals:
     Min       1Q   Median       3Q      Max
-21.7522  -3.4440   0.1651   3.7004  17.0517
Coefficients:
               Estimate Std. Error t value Pr(>|t|)
(Intercept)             9.76817    0.01306  747.66   <2e-16 ***
I(sin(2 * pi * x/365)) -1.17171    0.01827  -64.13   <2e-16 ***
I(cos(2 * pi * x/365)) 10.04347    0.01869  537.46   <2e-16 ***
I(sin(2 * pi * x))     -0.58321    0.01846  -31.59   <2e-16 ***
I(cos(2 * pi * x))      3.64653    0.01848  197.30   <2e-16 ***
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Residual standard error: 5.377 on 169486 degrees of freedom
Multiple R-Squared: 0.6617,     Adjusted R-squared: 0.6617
F-statistic: 8.286e+04 on 4 and 169486 DF,  p-value: < 2.2e-16

残留誤差は相変わらず大きい(5.377)のですが、私達はコロラドの天気は予測しにくいことで有名なのだと自分達を慰めることにします。

Rは時系列解析用に他にもツールがあります。例えばリビングルームの温度に対して自動相関関数をプロットすることもできます。

リスト14. リビングルームの温度に対する自動相関
> acf(ts(glarp$livingroom, frequency=(24*60/3)),
+     na.action=na.pass, lag.max=9*(24*60/3))
図2. リビングルームに対する自動相関グラフ
Autocorrelation graph for livingroom

埋め込みのts()ファンクションは、ベクターglarp$livingroomから時系列オブジェクトを作ります。サンプリングは一日当たりのサンプル数で規定されます。当然ですが、サンプリングを日数の整数倍にすると、温度には強い相関が現れます。また、7日間でわずかなピークがあることにも注意してください。これはBradのサーモスタットの設定が、(Bradが普通、日中家を留守にする)週末では、平日と異なっていることによるものです。これによって、7日間のウィンドウで、わずかに大きな相関が起きています。


まとめ

私達はここまでで、データ・セットの中にあるパターンや異常の解析にRを使ってきました。このデータ・セットには、実世界の科学データほとんど全てに共通するようなデータ欠落や潜在的な問題が含まれています。今回の記事ではまた、パターンやデータを素早く検証する上で、また解析シナリオで、Rの持つ機能プログラミング・スタイルがいかに有用かを調べてきました。第3回の記事では、大きなデータ・セットにおけるパターンを検出するための方法として、より高度な統計技法を使ってみます(それでもRの持つ豊富な機能のごく一部なのですが)。

参考文献

  • データ解析言語Rによる統計的プログラミング: 第 1 回 豊富な統計機能で遊ぶ」(developerWorks, 2004年9月)は R を紹介し、また基本的な統計分析との関連を解説しています。
  • データ解析言語Rによる統計的プログラミング: 第 3 回 再利用可能なオブジェクト指向プログラミング」(developerWorks, 2004年10月)は R のオブジェクト指向と R のその他の一般的なプログラミング概念を考察します。
  • Cameron Lairdによる、サーバー・クリニック: データとの格闘にはRが便利(developerWorks, 2003年)はR(とS)の概要を解説しており、参考資料も豊富に挙げられています。
  • R Project for Statistical Computingのホームページを見てください。このRのWebサイトにはチュートリアルから完璧なAPI記述まで、充実したドキュメンテーションが揃っています。Rが初めての読者にとっては、特に次の2つの資料が役に立つでしょう。
  • Davidによる魅力的なPython コラムの読者の何人かはPythonのユーザーであることから、PythonからRへのバインディングを特に気に入っています。実はバインディングにはRPyと、それよりも古いRSPythonの2種類があるのです。RSPythonも良いのですが、Davidの印象としてはRpyバインディングの方が良いようです。これらのバインディングいずれかを使い、Pythonオブジェクトをファンクションへの引き数として使うことで、Rのファンクションの全てをPythonコードから透過的に呼ぶことができます。
  • ご存じないかも知れませんが、ほとんどのシステムではRのコマンドラインでhelp.start()を入力すると、Rのドキュメンテーション用に生成されたHTMLページをブラウザーで開くことができます。
  • Summary of the history of S and S-Plus(SとS-Plusの歴史の要約)はオンラインで見ることができます。
  • DavidがCharming Python installment on 魅力的なPython コラムの記事として、Rと似たような使い心地や機能を数多く持ったNumerical Pythonについて書いています(ただし、Rの方がずっと機能が充実しています)。
  • DavidとBradがこの記事で使った温度データを、変換したり解析したりする練習をすることができます。(最初のログ・フォーマットを見やすい表形式に変換するPythonスクリプトもここにあります)
  • developerWorksのLinuxゾーンには、Linux開発者のための資料が他にも豊富に取り揃えられています。
  • Linux上で実行する、IBMミドルウェア製品の無料試用版をダウンロードしてください。developerWorksのSpeed-start your Linux appセクションから、WebSphere® Studio Application Developer, WebSphere Application Server, DB2® Universal Database, Tivoli® Access Manager, そしてTivoli Directory Serverが入手できます。またハウ・ツー記事や技術サポート情報もご利用ください。
  • developerWorks blogsに参加して、developerWorksコミュニティに加わってください。
  • Developer BookstoreのLinuxセクションではLinux関連の書籍が値引きして購入できますのでご利用ください。

コメント

developerWorks: サイン・イン

必須フィールドは(*)で示されます。


IBM ID が必要ですか?
IBM IDをお忘れですか?


パスワードをお忘れですか?
パスワードの変更

「送信する」をクリックすることにより、お客様は developerWorks のご使用条件に同意したことになります。 ご使用条件を読む

 


お客様が developerWorks に初めてサインインすると、お客様のプロフィールが作成されます。会社名を非表示とする選択を行わない限り、プロフィール内の情報(名前、国/地域や会社名)は公開され、投稿するコンテンツと一緒に表示されますが、いつでもこれらの情報を更新できます。

送信されたすべての情報は安全です。

ディスプレイ・ネームを選択してください



developerWorks に初めてサインインするとプロフィールが作成されますので、その際にディスプレイ・ネームを選択する必要があります。ディスプレイ・ネームは、お客様が developerWorks に投稿するコンテンツと一緒に表示されます。

ディスプレイ・ネームは、3文字から31文字の範囲で指定し、かつ developerWorks コミュニティーでユニークである必要があります。また、プライバシー上の理由でお客様の電子メール・アドレスは使用しないでください。

必須フィールドは(*)で示されます。

3文字から31文字の範囲で指定し

「送信する」をクリックすることにより、お客様は developerWorks のご使用条件に同意したことになります。 ご使用条件を読む

 


送信されたすべての情報は安全です。


static.content.url=http://www.ibm.com/developerworks/js/artrating/
SITE_ID=60
Zone=Linux
ArticleID=230110
ArticleTitle=データ解析言語Rによる統計的プログラミング: 第 2 回 機能プログラミングとデータ操作
publish-date=10072004