[アーティクル]
07/30/2015

2015 年 7 月

Volume 30 Number 7

R プログラミング言語 - C# プログラマ向けの R の手引き

R 言語 (以下 R) は、データの科学者やプログラマが使用する統計計算用の言語です。ソフトウェアシステムによって収集されるデータ量が増え、データを分析するニーズが高まっていることから、R は C# を扱う開発者の中で最も急速に広まっているテクノロジの 1 つです。この R に精通すれば、自身のスキルセットに貴重なスキルを加えることができるでしょう。

R は GNU プロジェクトの 1 つの無償ソフトウェアです。R は、1970 年代にベル研究所が作成した S 言語 ("Statistics" の頭文字) が起源です。R には優れたオンラインチュートリアルが数多くありますが、その大半は大学で統計を学ぶ学生を対象にしたものです。ここでは、C# プログラマができる限り早く R に慣れることを目指します。

今回の主旨を理解するには、サンプル R セッションを見るのが一番です (図 1 参照)。サンプルセッションには 2 つの独立したトピックがあります。前半にある少量のコマンド群は、一様分布のカイ 2 乗検定と呼ばれるものを示しています。2 つ目のコマンド群は線形回帰の例を示しています。線形回帰は、統計計算における Hello World 的なテクニック (私見) です。

図 1 サンプル R セッション

R の Web サイトは r-project.org (英語) にあります。このサイトには、R をダウンロードしてインストールできる複数のミラーサイトへのリンクがあります。インストールはシンプルな自己解凍型の実行ファイルで行われます。R は Windows XP 以降、および Windows 以外の著名なプラットフォームのほとんどで、公式にサポートされています。Windows 7 コンピューターにも Windows 8 コンピューターにも問題なく R をインストールできることを今回確認しています。既定では、インストールプロセスで 32 ビット版と 64 ビット版の両方が提供されます。

今回は、中級以上のプログラミングスキル (C# と R との類似点/相違点についての説明が理解できる程度) があることを前提にしていますが、R について何も知らなくても問題ありません。今回は C# のデモプログラムをすべて示していますが、本稿付属のコードダウンロードでも入手できます。

R を使用したカイ 2 乗検定

図 1 を見て最初に気付くことは、R の使い方が C# とは大きく異なることです。R のスクリプトを作成することはできますが、ほとんどの場合、コマンドシェルにおいて対話形式で使われます。最初の R の例は、普通の 6 面サイコロの目の出方に偏りがあるかどうかを確かめる分析です。サイコロに偏りがなければ、何度転がしても、6 つの目が出る結果の数は、おおよそ同数になると期待されます。

R のプロンプトはシェル内にトークン (>) で示されます。図 1 で入力された最初のステートメントは文字 (#) で始まっていますが、この文字はコメントを表す R のトークンです。

サンプルセッションの実際の最初のコマンドは次のとおりです。

> observed <- c(20, 28, 12, 32, 22, 36)

これにより、(連結用の) 関数 c を使用する、observed というベクトルが作成されます。ベクトルは同じデータ型のオブジェクトを保持します。これとほぼ等しい C# ステートメントは次のようになります。

var observed = new int[] { 20, 28, 12, 32, 22, 36 };

先頭の値 20 は 1 の目が出た回数、次の 28 は 2 の目が出た回数というように続きます。6 項目の合計は 20 + 28 + ..+ 36 = 150 になります。偏りのないサイコロであれば、期待されるそれぞれの結果は約 150/6 = 25 回です。しかし、観測された 3 の目の数 (12) の低さが問題になりそうです。

ベクトルを作成後、関数 chisq.test を使用してカイ 2 乗検定を行います。

> chisq.test(observed)

R では、読みやすい変数名や関数名を作成するために、ほとんどの場合、アンダースコア (_) ではなくドット (.) を使用します。関数 chisq.test を呼び出した結果はかなり長いテキストになります。

Chi-squared test for given probabilities

data:  observed
X-squared = 15.28, df = 5, p-value = 0.009231

C# の観点から考えると、大半の R 関数はデータ構造を返します。このデータ構造に含まれるデータの大半は無視できますが、Console.WriteLine に相当する出力を公開するステートメントも多く含まれています。R の出力の意味をどのように解釈するかは利用者次第です。今回の後半では、C# コードそのもの (ライブラリではなく) を使用して、同等のカイ 2 乗検定を作成する方法を示します。

今回の例の "X-squared" 値 15.28 は、計算されたカイ 2 乗統計です (ギリシャ文字のカイは大文字の X に似ていることからこのように表現されています)。値 0.0 は、サイコロに偏りがないと想定したときと観測値がまったく等しくなったことを表します。カイ 2 乗値が大きいほど、サイコロに偏りがない場合に、出目の観測値に偶然性がない可能性が高いことを示します。df 値 5 は "自由度" で、観測値の数より 1 小さくなります。df は今回の分析であまり重要ではありません。

p-value の 0.009231 は、サイコロに偏りがない場合の観測回数の偶然性の確率です。p-value が 1% 未満と非常に小さいため、観測値が偶然出現している可能性が非常に低いことから、問題のサイコロにはほぼ確実に偏りがあるという統計的根拠が存在すると結論づけられます。

R を使用した線形回帰分析

図 1 中の 2 つめのステートメント群は、線形回帰の例を示しています。線形回帰とは、数値変数 (統計では従属変数) と、数値にもカテゴリにもなり得る 1 つ以上の説明変数 (独立変数) との関係を表すために使用される統計手法です。独立した説明/予測変数が 1 つしかない場合、この手法を単線形回帰と呼びます。デモの例のように独立変数が 2 つ以上ある場合、この手法を多重線形回帰と呼びます。

線形回帰分析を行う前に、ディレクトリ C:\IntroToR に、DummyData.txt という名前で、8 つの項目を持つカンマ区切りのテキストファイルを作成します。ファイルの内容は次のとおりです。

Color,Length,Width,Rate
blue, 5.4, 1.8, 0.9
blue, 4.8, 1.5, 0.7
blue, 4.9, 1.6, 0.8
pink, 5.0, 1.9, 0.4
pink, 5.2, 1.5, 0.3
pink, 4.7, 1.9, 0.4
teal, 3.7, 2.2, 1.4
teal, 4.2, 1.9, 1.2

このファイルでは、花の色、花びらの長さと幅、成長率に関する花のデータを表しています。Color、Length、Width の各値から Rate 値 (最終列) を予測するというのが考え方です。コメントステートメントの後、次のような線形回帰分析に関する最初の R コマンドが 3 つあります。

> setwd("C:\\IntroToR")
> data <- read.table("DummyData.txt",
  header=TRUE, sep=",")
> print(data)

先頭のコマンドで作業ディレクトリを設定し、ソースデータファイルへのフルパスを指定しなくても済むようにしています。C# で一般的な (\\) トークンを使用する代わりに、Windows 以外のプラットフォームで一般的な (/) を使用してもかまいません。

次のコマンドは、data というテーブルオブジェクトのデータをメモリに読み込みます。R では名前付きパラメーターを使用します。header パラメーターは、データの 1 行目がヘッダー情報か (TRUE または短縮形の T) そうでないか (FALSE または F) を表します。R は大文字と小文字を区別します。通常、値の代入には (<-) 演算子または (=) 演算子のどちらかを使用できます。どちらを選択するかは好みの問題です。オブジェクトの代入には (<-) を、パラメーター値の代入には (=) を使うのが一般的です。

sep (セパレーター) パラメーターは、各行の値の分割方法を示します。たとえば、(\t) は値がタブで区切られていることを表し、("") は空白で区切られていることを示します。

関数 print はメモリ内のデータテーブルを表示します。関数 print にはオプションのパラメーターが数多くあります。図 1 の出力では、1 から始まるデータ項目のインデックスを表示しています。R では配列、行列、オブジェクトのインデックスが 1 から始まります。これに対して、C# 言語のインデックスは 0 から始まります。

線形回帰分析は以下の 2 つの R コマンドによって実行されます。

> model <- lm(data$Rate ~ (data$Color + data$Length + data$Width))
> summary(model)

最初のコマンドは、「予測のための従属変数をテーブルオブジェクト (data$Rate) の Rate 列とし、独立予測変数を Color、Length、Width として、関数 lm (線形モデル) の分析結果を model というオブジェクトに格納する」と解釈できます。2 つ目のコマンドは、「model というオブジェクトに格納された分析結果をそのまま表示する」という意味です。

関数 lm は膨大な量の情報を生成します。入力値を Color = pink、Length = 5.0、Width = 1.9 としたときの Rate 値を予測するとします (この入力値はデータ項目 [4] と同じ値のため、実際の Rate 値は 0.4 です)。予測を立てるには、列 Estimate の値を使用します。

Coefficients:
               Estimate Std. Error t value Pr(>|t|)
(Intercept)    -0.14758    0.48286  -0.306  0.77986
data$Colorpink -0.49083    0.04507 -10.891  0.00166 **
data$Colorteal  0.35672    0.09990   3.571  0.03754 *
data$Length     0.04159    0.07876   0.528  0.63406
data$Width      0.45200    0.11973   3.775  0.03255 *

X が従属変数値を表し、Y が予測後の Rate 値を表すと、以下のようになります。

X = (blue = NA, pink = 1, teal = 0, Length = 5.0, Width = 1.9)
Y = -0.14758 + (-0.49083)(1) + (0.35672)(0) + (0.04159)(5.0) + (0.45200)(1.9)
  = -0.14758 + (-0.49083) + (0) + (0.20795) + (0.85880)
  = 0.42834

予測した Rate (0.43) は実際のレート 0.40 にかなり近い値です。

これを言葉で表現すると、「モデルを使用して予測を立てるには、値 Estimate にそれぞれ対応する X 値を乗算した結果の線形和を求める」となります。値 Intercept は、どの変数にも関連付けられない定数です。カテゴリ形式の説明変数がある場合、値のうち 1 つ (この例では blue) が除外されます。

出力表示の最下部の情報は、独立変数の Color、Length、Width によって、従属変数の Rate がどれだけ精密に説明されるかを示します。

Residual standard error: 0.05179 on 3 degrees of freedom
Multiple R-squared:  0.9927,    Adjusted R-squared:  0.9829
F-statistic: 101.6 on 4 and 3 DF,  p-value: 0.00156

値 Multiple R-squared (0.9927) は、独立変数を直線状に組み合わせることで説明される従属変数の変化をパーセンテージにしたものです。少し別の言い方をすると、R-squared は 0 ～ 1 の値になり、この値が高いほどモデルの予測性が高いことになります。上記の結果では、R-squared 値が非常に高いため、Color、Length、Width によって Rate をかなり正確に予測できることが示されています。F-statistic、Adjusted R-squared 値、p-value は、モデルの適合性に関するその他の尺度です。

この例のポイントの 1 つは、R を使用してプログラミングする場合、この言語機能の裏側にある統計を理解することが群を抜いて難しいことです。R を段階的に学習する人のほとんどが、いくつか特定の問題を解決するニーズに応じて、1 つずつテクニックの知識を増やしています。C# を学習する人にも同じようなことが当てはまります。Hashtable や Queue クラスなど、Collections.Generic 名前空間にあるさまざまなコレクションオブジェクトについて学習する場合、たいていの開発者は 1 つずつデータ構造を学んでいきます。すべてのクラスについての知識を得てからプロジェクトで使用することを考える開発者はあまりいません。

他のカイ 2 乗検定

図 1 の種類のカイ 2 乗検定は、観測されたすべてのデータが等しい回数になるかどうか、つまりデータが一様分布しているかどうかをテストするため、多くの場合一様分布検定と呼ばれます。他にも、独立性検定と呼ばれるカイ 2 乗検定など、数種類のカイ 2 乗検定があります。

100 人から成るグループがあり、性別 (男性/女性) と所属政党 (民主党、共和党、その他の政党) との間に関係があるかどうかを調べるとします。データは分割表で表されるとします (図 2 参照)。表では、男性の方が女性よりも共和党員である確率が高いことが示されています。

図 2 分割表

	民主	共和	その他
男性	15	25	10	50
女性	30	15	5	50
	45	40	15	100

性別と所属の 2 つの要素が統計的に無関係かどうかを R で検定するには、まず次のコマンドで、データを数値行列に配置します。

> cm <- matrix( c(15,30,25,15,10,5), nrow=2, ncol=3 )

R では行列のデータが、C# のように行単位 (左から右へ、その後上から下へ) ではなく、列単位 (上から下へ、その後左から右へ) で格納されます。

次のコマンドで R のカイ 2 乗検定を行います。

> chisq.test(cm)

検定の結果、p-value が 0.01022 になるため、この 2 つの要素が無関係ではないことが 5% の有意水準で示されます。つまり、性別と所属政党には統計上関係があります。

最初のサイコロに関するカイ 2 乗検定の例の入力パラメーターはベクトルですが、今回の性別と所属に関する例では入力が行列です。関数 chisq.test には合計 7 つのパラメーターがあり、うちの 1 つ (ベクトルまたは行列) は必須で、残りの 6 つはオプションの名前付きパラメーターです。

Microsoft .NET Framework 名前空間の多くの C# メソッドと同様、R のほとんどの関数には多くのオーバーロードがあります。C# のオーバーロードでは、通常、同じ名前で複数のメソッドを実装し、異なるパラメーターを受け取るようにします。ジェネリックを使用するのもオーバーロードの 1 つの形式です。R のオーバーロードは 1 つの関数を使用して実装され、多くのオプションの名前付きパラメーターを使用します。

R によるグラフ

C# プログラマであれば、いくつかのプログラム出力データをグラフ化するとき、通常はプログラムを実行して出力したデータをコピーし、Ctrl + V でメモ帳に貼り付けてから不自然な制御文字を取り除き、Excel にそのデータをコピー、ペーストして Excel でグラフを作成します。この手順はやや煩雑ですが、ほとんどの状況で問題なく機能します。

R システムの強みの 1 つは、グラフを生成するネイティブな機能があることです。図 3 の例を見てください。このグラフは、Excel ではアドインがなければ作成できない種類のグラフです。

図 3 R を使用した 3D グラフ

図 1 で紹介したシェルプログラムに加え、R には準 GUI インターフェイスの RGui.exe もあり、グラフを作成する場合に使用できます。

図 3 のグラフは関数 z = f(x,y) = x * e^(-(x^2 + y^2)) を表します。最初の 3 つの R コマンドによってグラフを生成します。

> rm(list=ls())
> x <- seq(-2, 2, length=25)
> y <- seq(-2, 2, length=25)

関数 rm はメモリ内の現在のワークスペースからオブジェクトを削除します。使用しているコマンドは、すべてのオブジェクトを削除する R の魔法の呪文です。2 つ目と 3 つ目のコマンドは、-2 ～ +2 を両端を含めて 25 に均等分割したベクトルを作成します。代わりに関数 c を使うこともできます。

次の 2 つのコマンドは以下のとおりです。

> f <- function(x,y) { x * exp(-(x^2
  + y^2)) }
> z <- outer(x,y,f)

最初のコマンドは、R で function キーワードを使用して関数を定義する手法を示しています。outer という R の組み込み関数は、ベクトル x と y 、および関数定義 f を使用して値の行列を作成します。この結果は 25 行 25 列の行列となり、各セルの値は x と y に対応する関数 f の値になります。

次の 2 つのコマンドは以下のとおりです。

> nrz <- nrow(z)
> ncz <- ncol(z)

関数 nrow と ncol は、行列引数の列数と行数を返します。ここではどちらの値も 25 になります。

次の R コマンドは関数 colorRampPalette を使い、グラフに色付けするために、ユーザー設定の色のグラデーションパレットを作成します。

> jet.colors <- colorRampPalette(c("midnightblue", "blue",
+ "cyan", "green", "yellow", "orange", "red", "darkred"))

R で長い行を入力しているときに Enter キーを押すと、システムがカーソルを 1 つ下の行に移動し、コマンドがまだ完了していないことを示すプロンプトとして + 文字を表示します。次のコマンドは以下のとおりです。

> nbcol <- 64
> color <- jet.colors(nbcol)

上記 2 つのコマンドの結果は、非常に濃い青から、緑、黄色、非常に濃い赤までの範囲で 64 種類の色値をもつ、color というベクトルになります。次のコマンドは以下のとおりです。

> zfacet <- z[-1,-1] + z[-1,-ncz] + z[-nrz,-1] + z[-nrz,-ncz]
> facetcol <- cut(zfacet,nbcol)

図 3 のグラフを詳しく見ると、表面が 25 × 25 = 625 個の小さな四角形 (小平面) から成り立っているのがわかります。上記の 2 つのコマンドにより、各セルの値が対応する表面の小平面の色付けに使用される 64 色のどれか 1 つとなる、25 行 25 列の行列が作成されます。

最後に関数 persp (透視図) を使用して、3D グラフを表示します。

> persp(x,y,z, col=color[facetcol], phi=20, theta=-35,
+ ticktype="detailed", d=5, r=1, shade=0.1, expand=0.7)

関数 persp にはオプションの名前付きパラメーターがたくさんあります。パラメーター col は、(使用する) 色を指します。パラメーター phi と theta はグラフに対する視野角 (左右および上下) です。パラメーター ticktype は、x、y、z 軸の値の表示形式を制御します。パラメーター r と d は、視点からグラフまでの距離と、視覚の 3D 効果を制御します。パラメーター shade は、仮想光源からの陰影のシミュレーションを制御します。パラメーター expand はグラフの幅と高さの比率を制御します。関数 persp にはこれ以外にもさまざまなパラメーターがありますが、大半の状況は上記のパラメーターで十分対応できます。

上記の例から、R には強力で柔軟なネイティブのグラフ描画機能があっても、比較的機能レベルが低く、かなりの作業を要することがわかります。

C# でのカイ 2 乗検定

R 言語と C# 言語のプログラミングの類似点と相違点を理解するためには、C# でカイ 2 乗検定を実装してみるのが効果的です。また、こうした C# コードは、個人のコードライブラリの強化に役立ちます。図 4 の C# デモプログラムを見てください。

図 4 C# を使用したカイ 2 乗検定

このデモプログラムは、図 1 で示したサイコロに関する R 言語のカイ 2 乗検定に似せたものです。C# デモの出力値は、R セッションの値とまったく同じです。

デモプログラムを作成するには、Visual Studio を起動して、ChiSquare という新しい C# コンソールアプリケーションプロジェクトを作成します。テンプレートコードがエディターに読み込まれたら、ソリューションエクスプローラーウィンドウで Program.cs ファイルを右クリックし、名前を「ChiSquareProgram.cs」に変更します。Program クラスの名前が Visual Studio によって自動的に「ChiSquareProgram」に変更されます。

図 5 の完全な C# デモプログラのリストを参考にしてください。ソースコードは予想以上に長くなります。C# そのものを使って統計プログラムを実装するのはそれほど難しくありませんが、コードが長くなる傾向があります。

図 5 C# のカイ 2 乗デモプログラム

using System;
namespace ChiSquare
{
  class ChiSquareProgram
  {
    static void Main(string[] args)
    {
      try
      {
        Console.WriteLine("\nBegin Chi-square test using C# demo\n");
        Console.WriteLine(
          "Goal is to see if one die from a set of dice is biased or not\n");
        int[] observed = new int[] { 20, 28, 12, 32, 22, 36 };
        Console.WriteLine("\nStarting chi-square test");
        double p = ChiSquareTest(observed);
        Console.WriteLine("\nChi-square test complete");
        double crit = 0.05;
        if (p < crit)
        {
          Console.WriteLine("\nBecause p-value is below critical value of " +
            crit.ToString("F2"));
          Console.WriteLine("the null hypothsis is rejected and we conclude");
          Console.WriteLine("the data is unlikely to have happened by chance.");
        }
        else
        {
          Console.WriteLine("\nBecause p-value is not below critical value of " +
            crit.ToString("F2"));
          Console.WriteLine(
            "the null hypothsis is accepted (not rejected) and we conclude");
          Console.WriteLine("the observed data could have happened by chance.");
        }
        Console.WriteLine("\nEnd\n");
        Console.ReadLine();
      }
      catch (Exception ex)
      {
        Console.WriteLine(ex.Message);
        Console.ReadLine();
      }
    } // Main
    static void ShowVector(int[] v)
    {
      for (int i = 0; i < v.Length; ++i)
        Console.Write(v[i] + " ");
      Console.WriteLine("");
    }
    static double ChiSquareTest(int[] observed)
    {
      Console.WriteLine("Observed frequencies are: ");
      ShowVector(observed);
      double x = ChiSquareStatistic(observed);
      Console.WriteLine("\nX-squared = " + x.ToString("F2"));
      int df = observed.Length - 1;
      Console.WriteLine("\ndf = " + df);
      double p = ChiSquareProb(x, df);
      Console.WriteLine("\np-value = " + p.ToString("F6"));
      return p;
    }
    static double ChiSquareStatistic(int[] observed)
    {
      double sumObs = 0.0;
      for (int i = 0; i < observed.Length; ++i)
        sumObs += observed[i];
      double expected = (int)(sumObs / observed.Length);
      double result = 0.0;
      for (int i = 0; i < observed.Length; ++i)
      {
        result += ((observed[i] - expected) *
         (observed[i] - expected)) / expected;
      }
      return result;
    }
    public static double ChiSquareProb(double x, int df)
    {
      // x = a computed chi-square value. df = degrees of freedom.
      // output = prob. the x value occurred by chance.
      // So, for example, if result < 0.05 there is only a 5% chance
      // that the x value occurred by chance and, therefore,
      // we conclude that the actual data which produced x is
      // NOT the same as the expected data.
      // This function can be used to create a ChiSquareTest procedure.
      // ACM Algorithm 299 and update ACM TOMS June 1985.
      // Uses custom Exp() function below.
      if (x <= 0.0 || df < 1)
        throw new Exception("parameter x must be positive " +
        "and parameter df must be 1 or greater in ChiSquaredProb()");
      double a = 0.0; // 299 variable names
      double y = 0.0;
      double s = 0.0;
      double z = 0.0;
      double e = 0.0;
      double c;
      bool even; // Is df even?
      a = 0.5 * x;
      if (df % 2 == 0) even = true; else even = false;
      if (df > 1) y = Exp(-a); // ACM update remark (4)
      if (even == true) s = y; else s = 2.0 * Gauss(-Math.Sqrt(x));
      if (df > 2)
      {
        x = 0.5 * (df - 1.0);
        if (even == true) z = 1.0; else z = 0.5;
        if (a > 40.0) // ACM remark (5)
        {
          if (even == true) e = 0.0;
          else e = 0.5723649429247000870717135;
          c = Math.Log(a); // log base e
          while (z <= x)
          {
            e = Math.Log(z) + e;
            s = s + Exp(c * z - a - e); // ACM update remark (6)
            z = z + 1.0;
          }
          return s;
        } // a > 40.0
        else
        {
          if (even == true) e = 1.0;
          else e = 0.5641895835477562869480795 / Math.Sqrt(a);
          c = 0.0;
          while (z <= x)
          {
            e = e * (a / z); // ACM update remark (7)
            c = c + e;
            z = z + 1.0;
          }
          return c * y + s;
        }
      } // df > 2
      else
      {
        return s;
      }
    } // ChiSquare()
    private static double Exp(double x) // ACM update remark (3)
    {
      if (x < -40.0) // ACM update remark (8)
        return 0.0;
      else
        return Math.Exp(x);
    }
    public static double Gauss(double z)
    {
      // input = z-value (-inf to +inf)
      // output = p under Normal curve from -inf to z
      // e.g., if z = 0.0, function returns 0.5000
      // ACM Algorithm #209
      double y; // 209 scratch variable
      double p; // result. called 'z' in 209
      double w; // 209 scratch variable
      if (z == 0.0)
        p = 0.0;
      else
      {
        y = Math.Abs(z) / 2;
        if (y >= 3.0)
        {
          p = 1.0;
        }
        else if (y < 1.0)
        {
          w = y * y;
          p = ((((((((0.000124818987 * w
            - 0.001075204047) * w + 0.005198775019) * w
            - 0.019198292004) * w + 0.059054035642) * w
            - 0.151968751364) * w + 0.319152932694) * w
            - 0.531923007300) * w + 0.797884560593) * y * 2.0;
        }
        else
        {
          y = y - 2.0;
          p = (((((((((((((-0.000045255659 * y
            + 0.000152529290) * y - 0.000019538132) * y
            - 0.000676904986) * y + 0.001390604284) * y
            - 0.000794620820) * y - 0.002034254874) * y
            + 0.006549791214) * y - 0.010557625006) * y
            + 0.011630447319) * y - 0.009279453341) * y
            + 0.005353579108) * y - 0.002141268741) * y
            + 0.000535310849) * y + 0.999936657524;
        }
      }
      if (z > 0.0)
        return (p + 1.0) / 2;
      else
        return (1.0 - p) / 2;
    } // Gauss()
  } // Program class
} // ns

Main メソッドは大部分が WriteLine ステートメントで構成されています。基本は以下のコード呼び出しで成り立っています。

int[] observed = new int[] { 20, 28, 12, 32, 22, 36 };
double p = ChiSquareTest(observed);
double crit = 0.05;
if (p < crit) {
  // Messages
} else {
  // Messages
}

C# のカイ 2 乗検定は観測値の配列を受け取り、カイ 2 乗の統計値と自由度を計算および表示し、さらに p-value を計算して返します。

ChiSquareTest メソッドは次の 3 つのヘルパーメソッドを呼び出します。

ShowVector(observed);
double x = ChiSquareStatistic(observed);
int df = observed.Length - 1;
double p = ChiSquareProb(x, df);

ShowVector メソッドは、入力パラメーター値をエコーする R 関数 chisq.test のアプローチに似た方法で、入力ベクトルを表示します。ChiSquareStatistic メソッドは、計算したカイ 2 乗 (R の出力では "X-squared") を返し、ChiSquareProb メソッドは ChiSquareStatistic から返ってきた結果を使用して確率 (R の出力では "p-value") を計算します。

ChiSquareStatistic メソッドは、一様分布用の単純な検定です。カイ 2 乗統計量を求める統計式は次のようになります。

chi-squared = Sum( (observed - expected)^2 / expected )

そのため、カイ 2 乗を計算する前に、観測値に関連付けられる期待値を計算する必要があります。前述のように、観測された配列のカウント値を合計し (デモでは 150)、配列中の値の個数 (デモでは 6) でその合計を除算することで、すべてのセルに期待値 (25) を提供します。

カイ 2 乗検定を作成する際に最も難しいのは、p-value の計算です。図 5 の ChiSquareProb メソッドの定義に目を通すとすぐに、その理解には深い専門知識が必要なことがわかります。さらに、このメソッドは Gauss というヘルパーメソッドを呼び出しますが、こちらのメソッドも同じく複雑です。

実はほとんどの関数がここ数十年で解決されたため、複雑な数値関数でも非常に簡単にコーディングできます。中でも、Association for Computing Machinery (ACM) が収集したアルゴリズムのリポジトリと、縮めて "A&S" と呼ばれることで有名な、Abramowitz と Stegun による『Handbook of Mathematical Functions』(Dover Publications、1965) の 2 つのリソースが最もよく使われています。どちらの資料も、Web 上の複数のサイトから無償で入手できます。

たとえば、ChiSquareProb メソッドは ACM Algorithm #299、Gauss ヘルパーメソッドは ACM Algorithm #209 をそれぞれ実装しています。 ACM #209 は、A&S の equation #7.1.26 に簡単なラッピングステートメントを加えたものです。

ACM Algorithms と A&S の関数の多くは、既存の C# ライブラリに実装されていますが、そのようなライブラリの大半は非常に大きくなります。可能であれば、外部に頼らないようにしながら必要なメソッドを使い、コードを最初から作成することをお勧めします。

コメントをいくつか

本稿で紹介したのは R 言語のほんの一部分ですが、使いこなせるようになるためには十分な内容だと考えています。C# プログラマは、次の 2 つの状況で R に触れる機会が多くなります。まず、独自で統計分析を行う場合に R を使用するような状況です。今回述べたように、背景にある統計を理解すれば、R は非常に簡単な言語です。次に、R を主な言語として使用する人と関わる必要があるために、R の環境と用語を理解しなければならない状況です。

開発ツールにはわずかながら R 言語が統合されているものがあり、それを活用できます。有名なのは、RStudio と呼ばれるプロジェクトです。一般的にはネイティブな R コンソールが好んで使われます。複数のマイクロソフト製品が R に対応しています。たとえば、Microsoft Azure Machine Learning サービスは、R 言語のスクリプトを使用することができます。またいずれ、Visual Studio に R へのサポートが追加されるとの推測もあります。

Dr. James McCaffrey は、ワシントン州レドモンドにある Microsoft Research に勤務しています。これまでに、Internet Explorer、Bing などの複数のマイクロソフト製品にも携わってきました。McCaffrey 博士の連絡先は、jammc@microsoft.com (英語のみ) です。

この記事のレビューに協力してくれたマイクロソフト技術スタッフの Dan Liebling および Robin Reynolds-Haertle に心より感謝いたします。
Robin Reynolds-Haertle は Visual Studio によるクロスプラットフォーム開発に関するドキュメントを執筆しています。彼女は、現在 .NET Core、C#、Swift、および R に関心を寄せています。