データ解析言語Rによる統計的プログラミング: 第 1 回 豊富な統計機能で遊ぶ

生データを形にする

3回シリーズの第1回として、無料のソフトウェアとしてリリースされており、機能豊富な統計環境であるRをDavidとBradが紹介します。Rにはプログラミング言語が含まれており、また対話型のシェルや充実したグラフ機能があります。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年 9月 21日

R環境はいわゆるプログラミング言語として意図されたものではなく、むしろ対話型のツールと言うべきもので、データ・プロパティに対する広範なグラフィック表現の生成を含めて、データ・セットを解析調査するためのものです。生成されたグラフィックスも、セッション中に行ったステップのどちらも、後で使うために保存しておくことができます。これは、プロジェクト毎にそれまで様々に設定した動作環境を選択する上で特に便利です。デフォルトでRのコマンドはセッション・ヒストリーに保存されますが、特別に有効な命令シーケンスを.Rファイルに保存し、(source()を使って)セッション内でソースとして使うことができます。

Rの設計者達はその目標を「An Introduction to R」(Rの紹介:下記の参考文献に全文へのリンクがあります)の中で次のように述べています。

Rで行なう解析には別々の作業ディレクトリーを使うように推奨します。解析中にはxやyという名前を持ったオブジェクトがごく普通に生成されます。こうした名前は、単一の解析の場合には意味のある場合が多いのですが、同じディレクトリーで幾つかの解析が行われている時には、それが何かを判別するのは困難なものになりがちです。

各ディレクトリーには、セッション中の作業オブジェクトのバイナリー直列化を含む隠しファイルが生成されます。これによって、前回使ったアクティブ変数を全て持ったセッションが再開できるようになります。

(GPLライセンスの)Rプログラミング言語には2つの親があります。その一方である独自なS/S-PLUSプログラミング言語が、Rの構文の大部分の元になっており、またもう一方のSchemeプログラミング言語は、多くの(より詳細な)意味体系の元となっています。Sの起源は1984にさかのぼり、(S-PLUSを含めて)その後のバージョンでは多くの機能強化が行われています。Schemeはもちろん、(Lispと同様)ずっと大昔までさかのぼります。Rは1997年にS/S-frPLUSのフリー・ソフト・スーパーセットとして合体し、それ以来、活発なユーザー・コミュニティや開発者コミュニティができています。ただしRの素晴らしさを利用するにあたって、その歴史的価値を恐れる必要はありません。Cameron LairdがdeveloperWorksの記事データとの格闘にはRが便利 の中でRの背景や関連資料を説明しています。

RはLinuxやWindows、Mac OS XそれにMac OS Classicなど、様々なコンピューター・プラットフォーム用に実行ファイルが入手可能です。当然ながら、他のプラットフォームへのコンパイル用にソースを入手することも可能です(例えばこの記事の共著者であるBradは、FreeBSD用に難なくRをビルドしました)。ただし様々なプラットフォームで、Rには少し問題もあります。例えば、DavidのMac OS XマシンでQuartzを使ってプロット出力すると、反応のない表示ウィンドウができてしまいます。もっと悪いのはBradのFreeBSD/AMD AthlonマシンでRを終了すると、強制的にリブートをかけてしまうのです(これは恐らくSSEカーネル・オプションが正しくないことと関係しているのですが、それにしてもこの振る舞いには問題があります)。とは言えRは一般的には安定、高速であり、しかも統計や数学的機能の幅広さは驚くほどです。オプションのパッケージを追加すると、膨大な標準パッケージに対して、さらに幅広い機能が追加されることになります。

Rのデータ・モデル

Rでの基本的なデータ・オブジェクトはベクターです。ベクターの変種の幾つかによって、(多次元)配列やデータ・フレーム、(異種)リスト、マトリックスなどのような機能が追加されます。NumPy/NumArrayやMatlabなどと同様、ベクターやその兄弟分に関する操作は、メンバー・データに関して要素単位(elementwise)で行われます。Rが実際に使われている例をざっと見てみると、その構文の様子がつかめるでしょう(この最初の例では、シェル・プロンプトやレスポンスも含んでいます)。

リスト1. ベクターと要素単位操作
> a <- c(3.1, 4.2, 2.7, 4.1)  # Assign with "arrow" rather than "="
> c(3.3, 3.4, 3.8) -> b       # Can also assign pointing right
> assign("c", c(a, 4.0, b))      # Or explicitly to a variable name
> c                              # Concatenation "flattens" arguments
[1] 3.1 4.2 2.7 4.1 4.0 3.3 3.4 3.8
> 1/c                            # Operate on each element of vector
[1] 0.3225806 0.2380952 0.3703704 0.2439024 0.2500000 0.3030303 0.2941176
[8] 0.2631579
> a * b                          # Cycle shorter vector "b" (but warn)
[1] 10.23 14.28 10.26 13.53
Warning message:
longer object length
        is not a multiple of shorter object length in: a * b
> a+1                            # "1" is treated as vector length 1
[1] 4.1 5.2 3.7 5.1

R構文でのインデックス化、スライシング、名前付き引き数やオプション引き数、その他の要素については下記に挙げます。Rのシェル・プロンプトは、特にGNU readlinesをインストールしてある場合は、様々なことを探るには素晴らしいインターフェースです。作業をしながら詳しく学ぶには、help(function)コマンドを使うことを頭に置いておいてください(?functionを使うこともできます)。PythonシェルのユーザーにとってはRシェルはすぐに分かるでしょうし、どちらのユーティリティも便利に使えるでしょう。


温度データ・セット

Bradは自分の家の内外に置いてある4つの温度計の温度をほとんど一年に渡って集めており、スライド式窓のような読み取り値を自動的にコンパイルし、GnuPlotを使ってWebからアクセスできるグラフを作っています(Gnuplotに関する詳しい情報については参考文献にリンクがあります)。こうした、ハッカーのようなデータ集合は汎用的な科学用途にはあまり意味がないかも知れませんが、科学データに類似している幾つかの優れた特徴もあるのです。

データは3分ごとに収集されるので、一年間では大きなデータ点の集まりになります(4つの測定点で750,000)。温度計や伝送チャネルの故障、データを記録しているコンピューターの故障などのために、幾つかのデータは欠けています。また、稀ではありますがタイミング・エラーのために、単線伝送チャネルでは同時に読み取った値の順序を変えて伝送する場合があることが知られています。つまり、かなり良いのだが多少欠けている点や不完全な点がある、という意味で、Bradの温度データは実世界の科学データに非常によく似ているのです。

データを読む

温度データは測定点をファイル名とする4つの別々なデータ・ファイルに集められ、それぞれのデータは次のようなフォーマットを持っています。

リスト2. 元々の温度データ・ファイルのフォーマット
2003 07 25 16 04 27.500000
2003 07 25 16 07 27.300000
2003 07 25 16 10 27.300000

このデータを読むための最初のパスは次のようなものになります。

リスト3. 温度データを読むための最初のパス
> lab <- read.fwf('therm/lab', width=c(17,9))   # Fixed width format
> basement <- read.fwf('therm/basement', width=c(17,9))
> livingroom <- read.fwf('therm/livingroom', width=c(17,9))
> outside <- read.fwf('therm/outside', width=c(17,9))
> l_range <- range(lab[,2]) 	  # Vector of data frame: entire second column
> b_range <- range(basement[,2])    # range() gives min/max
> v_range <- range(livingroom[,2])
> o_range <- range(outside[,2])
> global <- range(b_range, l_range, v_range, o_range)
> global          			  # Temperature range across all sites
[1] -19.8  32.2

データをきれいにする

純粋な初期データのフォーマットには少し問題があります。その一つとして、欠けているデータが明示的に示されておらず、単純に、行とタイムスタンプが無いことでしか判別できません。さらに、日付が(ISO8601/W3Cではなく)空白の入った非標準フォーマットで保存されています。また、小さなことですが、4つのファイルでタイムスタンプを繰り返しているのは記憶容量の点で非効率です。R自体でデータをきれいにすることは当然できるのですが、私達は「R Data Import/Export」(参考文献)という文書にある、Rの設計者達の助言に従うことにしました。一般的にテキスト処理を一番うまく行うには、テキスト処理に特化した言語を使うことです。私達の場合はRを使ってそのまま読み取れる統一データファイルを生成するために、Pythonスクリプトを書きました。例えば、新しいデータファイルglarp.tempsの初めの数行は次のようなものになります。

リスト4. 統一温度データ・フォーマット
timestamp      basement     lab     livingroom  outside
2003-07-25T16:04 24.000000      NA     29.800000  27.500000
2003-07-25T16:07 24.000000      NA     29.800000  27.300000
2003-07-25T16:10 24.000000      NA     29.800000  27.300000

改善されたデータセットで作業してみることにしましょう。

リスト5. 統一温度データ・フォーマットで作業する
> glarp <- read.table('glarp.temps', header=TRUE, as.is=TRUE)
> timestamps <- strptime(glarp[,1], format="%Y-%m-%dT%H:%M")
> names(glarp)              # What column names were detected?
[1] "timestamp"  "basement"   "lab"        "livingroom" "outside"
> class(glarp[,'basement'])  # Kind of data is in basement column?
[1] "numeric"
> basement <- glarp[,2]         # index by position
> lab <- glarp[,'lab']          # index by name
> outside <- glarp$outside      # equiv to prior indexing
> livingroom <- glarp$living    # name with unique initial name
> summary(glarp)  		# Handy built-in to describe most R objects
  timestamp            basement            lab            livingroom
 Length:171349      Min.   :   6.40   Min.   :  -6.40   Min.   :   7.20
 Class :character   1st Qu.:  17.00   1st Qu.:  16.60   1st Qu.:  18.10
 Mode  :character   Median :  19.10   Median :  17.90   Median :  20.30
                    Mean   :  18.88   Mean   :  18.12   Mean   :  20.17
                    3rd Qu.:  20.50   3rd Qu.:  19.50   3rd Qu.:  22.00
                    Max.   :  27.50   Max.   :  25.50   Max.   :  31.30
                    NA's   :1854.00   NA's   :2406.00   NA's   :1855.00
    outside
 Min.   : -19.800
 1st Qu.:   2.100
 Median :   9.800
 Mean   :   9.585
 3rd Qu.:  17.000
 Max.   :  32.200
 NA's   :1858.000

基本的な統計解析

range()ファンクションは既に見ました。min()max()は、データの範囲での最小値と最大値を見つけます。summary()も明らかにこの情報を表示するのですが、他の計算で直接使える形式では表示しません。このデータの持つ、非常に基本的な統計プロパティをさらに見つける所から始めましょう。

リスト6. 温度データに対する基本的な統計計算
> mean(basement)        	# Mean fails if we include unavailable data
[1] NA
> mean(basement, na.rm=TRUE)
[1] 18.87542
> sd(basement, na.rm=TRUE)  	# Standard deviation must also exclude NA
[1] 2.472855
> cor(basement, livingroom, use="all.obs")   # All observations: no go
Error in cor(basement, livingroom, use = "all.obs") :
        missing observations in cov/cor
> cor(basement, livingroom, use="complete.obs")
[1] 0.9513366
> cor(outside, livingroom, use="complete.obs")
[1] 0.6446673

直感的に考えられる通り、2つの室内温度データの間には、戸外の温度どちらよりも相関があります。それにしても、チェックは簡単です。

温度の分布

平均と標準偏差は見ましたので、皆さんは直感的に温度が正規分布していると期待するかも知れません。これを見てみましょう。

リスト7. 短い一行でヒストグラムを生成する
> hist(outside)

様々なRのコマンドによって、データ・セットをプロットしたもの、グラフ化したもの、図表化したものなどを含んだ2番目のウィンドウがポップアップします。これがどのように行われるかの詳細はプラットフォームやそれぞれの設定によって異なります。こうしたグラフィックスは、後で使えるように外部ファイルにリダイレクトすることもできます。上記のコマンドは、次のようなhist()グラフを生成します。

図1.
Default histogram of outside temperatures

最初としては悪くありません。幾つかのパラメーターを使うと、丸めのしきい値を狭めることができます

リスト8. ヒストグラムの丸め密度を変更する
hist(outside, breaks=1000, border="blue")
図2.
Dense histogram of outside temperatures

密度の濃いヒストグラムにはデータの荒れがあり、大まかに7度から12度の間の範囲で特定な測定値の頻度が非常に高く、またその他に、予想外に低い頻度もあることに気がつきます。こうした大きな不連続は、恐らく温度計の特性によるサンプルの偏りを示しているのだと思います。その一方、ちょうどサーモスタットで調節された室内温度に近い24度付近で、狭いながら頻度が突出しています。これは温度計の伝送チャネルに関して先に述べた、測定値の伝送順序が変えられたことによる可能性が高いようです。いずれにせよ、こうしたグラフィック表示を見ると、さらに詳しく調べたり解析したりしたくなるものが見えてきます。

2、3他のグラフも見てみると、室内の温度分布が分かります。

リスト9. リビング・ルーム温度のヒストグラム
> hist(livingroom, breaks=40, col="blue", border="red")
> hist(livingroom, breaks=400, border="red")
図3.
40 step histogram of living room temperatures
図4.
400 step histogram of living room temperatures

リビング・ルームの温度分布は、戸外よりも妥当なようです。分解能を上げると不連続な点がいくつか現れますが、これは温度計の持つ、規模の小さな偏りによるもののようです。しかし一般的なパターンとしては、(Bradの)タイマー制御によるサーモスタットから想定される3峰分布(trimodal distribution)に従っています(大きな山が21度付近にあり、小さな山が16度と24度付近にあります)。


データをさらに視覚化する

各測定点は温度の値のリニアー・ベクターですが、データには、日と年という、主に2つの周期があることが直感的に想定されます(つまり夜や冬は寒い)。

最初の問題は1次元のデータを2次元のデータ点マトリックスに変換することです。次にこの2次元データ・セットを視覚化することにします。

リスト10. ベクターを再整形し、温度をプロットする
> oarray <- outside[1:170880] # Need to truncate a few last-day readings
> dim(oarray) <- c(480,356)   # Re-dimension the vector to a 2-D array
> plot(oarray[1,], col="blue", type="l", main="4 p.m. outside temp",
+      xlab="Day of year (starting July 25, 2003)",
+      ylab="Temperature (Celsius)")
図5.
4 p.m. outside temperature

ベクターを一旦「時間 X 日」マトリックス(2次元配列)に変換すれば、それぞれの日に対する温度から年間パターンのグラフを作るのが自然なのですが、他の方法でもできます。つまり、ベクターから480番目毎の点を抽出するのです。Rによるこの方法は、ずっとすっきりしています。

3次元データ

丸1年間の温度測定値の表現はどうでしょう? 一つの方法として、次のように色でコード化した熱分布グラフを作ることができます。

リスト11. 熱分布グラフを作る
> x <- 1:480                  # Create X axis indices
> y <- 1:356                  # Create Y axis indices
> z <- oarray[x,y]            # Define z-axis (really same as oarray)
> mycolors <- c(heat.colors(33),topo.colors(21))[54:1]
> image(x,y,z, col=mycolors,
+       main="Outside temperatures near Glarp", # Brad's name for his house
+       xlab="Minutes past 4 p.m",
+       ylab="Days past July 25, 2003" )
> dev2bitmap(file="outside-topo.pdf", type="pdfwrite")
図6.
Thermal map of daily x yearly temperatures

今ここで行ったことに関して幾つか注釈を付けておきましょう。座標や点の定義は、数字のリストを作るためのPython風のスライス操作が一旦分かってしまえば、極めて明確です。zを作る上でxyでインデックスをつけると、そのインデックスの幅と高さを持った配列ができます。この場合では、zoarrayと同じで、どうでも良いのですが、値や、値に到達するために使われるオフセットをシステム的に変更するのはごく簡単です。カラー・マップmycolorを作るには、少し試行錯誤をしています。赤や黄色を使うと「暖かい」温度(つまり摂氏0度以上)には良いようですが、低い温度には適切ではないようです。そこでカラー・ベクターに少し青/緑を連結させます。私達は標準のcolormapファンクションで生成される順序と逆の順序でも色が欲しかったのですが、インデックスによると、これも簡単でした。

熱分布マップはこれまでのグラフに比べて少しシャープに描かれていることに気が付くかも知れません。image()コマンドを使うと、適切ながらこれほど印象的ではないグラフがスクリーン上に描かれます。ここでも分かる通り、「current image」を外部ファイルにエクスポートした方が良い結果が得られる場合がよくあります。

私達の個人的な好みとしては先に示した平坦な熱分布マップが好きなのですが、擬似遠近法で3次元データにした方がうまく情報が伝わると思う人も多いかも知れません。Rでは、ほとんど手をかけることなく、驚くような3次元データの立体遠近法表現が生成できるのです。

リスト12. 立体的な表面グラフを作る
> persp(x,y,z, theta=10, phi=60, ltheta=40, lphi=30, shade=.1, border=NA,
+       col=mycolors[round(z+20)], d=.5,
+       main="Outside temperatures near Glarp",
+       xlab="Minutes past 4 p.m",
+       ylab="Days past July 25, 2003",
+       zlab="Temperature", )
図7.
Topographic map of daily x yearly temperatures

まとめ

この記事で私達が行った幾つかの統計的解析や生成した基本的なグラフは、Rの持つ統計機能の、ほんのごく一部にすぎません。ほとんどあらゆる科学分野に対して、あるいは(知られているもの、知られていないものも含めて)ほとんど全ての統計手法に対して、標準機能が用意されているか、またはその機能に関連した数学的手法をサポートする拡張パッケージがあるのです。

この記事で、Rを使っての作業がどんなものかを理解できたと思います。Rにはここで紹介した他にも数多くの機能を持っています。Rに関するシリーズ第2回目の次回は、中級から上級の手法として、線形、非線形回帰に入って行くことにします。

参考文献

  • データ解析言語Rによる統計的プログラミング: 第 2 回 機能プログラミングとデータ操作」(developerWorks, 2004年10月)は第 1 回への追加として、この言語の機能に入り込んでいます。
  • データ解析言語Rによる統計的プログラミング: 第 3 回 再利用可能なオブジェクト指向プログラミング」(developerWorks, 2004年10月)は R のオブジェクト指向と R のその他の一般的なプログラミング概念を考察します。
  • R Project for Statistical Computingのホームページを見てください。このRのWebサイトにはチュートリアルから完璧なAPI記述まで、充実したドキュメンテーションが揃っています。Rが初めての読者にとっては、特に次の2つの資料が役に立つでしょう。
  • Davidによる魅力的なPythonコラムの読者の何人かはPythonのユーザーであることから、PythonからRへのバインディングを特に気に入っています。実はバインディングにはRPyと、それよりも古いRSPythonの2種類があるのです。RSPythonも良いのですが、Davidの印象としてはRpyバインディングの方が良いようです。これらのバインディングいずれかを使い、Pythonオブジェクトをファンクションへの引き数として使うことで、Rのファンクションの全てをPythonコードから透過的に呼ぶことができます。
  • Davidが魅力的なPythonコラムの記事として、Rと似たような使い心地や機能を数多く持ったNumerical Pythonについて書いています(ただし、Rの方がずっと機能が充実しています)。
  • ご存じないかも知れませんが、ほとんどのシステムではRのコマンドラインでhelp.start()を入力すると、Rのドキュメンテーション用に生成されたHTMLページをブラウザーで開くことができます。
  • summary of the history of S and S-Plus(SとS-Plusの歴史の要約)はオンラインで見ることができます。
  • Cameron Lairdによる、サーバー・クリニック: データとの格闘にはRが便利(developerWorks, 2003年)はR(とS)の概要を解説しており、参考資料も豊富に挙げられています。
  • Gnuplotもデータをプロットしたり解析したりするための役に立ちます。Nishanth Sastryが記事Visualize your data with gnuplot(developerWorks, 2004年)の中でGnuplotを紹介し、欠けているデータ点の問題を解説しています。
  • データに関する基準は存在するのですが、あまり知られていません。例えばBradは初めからNational Space Science Data Center(米国宇宙科学データセンター)の Common Data Format (CDF) を使った方が賢明だったのかも知れません。ところがそれは良くある間違いです。皆さんも同じような状況になったら、また生データをRやNumPy、Gnuplotなどのようなツールで使えるような形式に変換する必要があるならば、PerlやPython、それにGNUなどのテキスト・ユーティリティが役に立ちます。この最後のGNUに関して、DavidがdeveloperWorksのチュートリアルUsing the GNU text utilities(developerWorks, 2004年)で紹介しています。
  • DavidとBradがこの記事で使った温度データを、変換したり解析したりする練習をすることができます。(最初のログ・フォーマットを見やすい表形式に変換するPythonスクリプトもここにあります)
  • developerWorksのLinuxゾーンには、Linux開発者のための資料が他にも豊富に取り揃えられています。
  • Linux上で実行する、選りすぐりのdeveloperWorks Subscription製品の無料の試用版をダウンロードして下さい。developerWorksのSpeed-start your Linux appセクションからWebSphere Studio Site DeveloperやWebSphere SDK for Web services、WebSphere Application Server、DB2 Universal Database Personal Developers Edition、Tivoli Access ManagerそれにLotus Domino 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=230109
ArticleTitle=データ解析言語Rによる統計的プログラミング: 第 1 回 豊富な統計機能で遊ぶ
publish-date=09212004