SAS でアニメーション作成 ~中心極限定理のデモ~

なぜ中心極限定理か

SAS は 9.4 以降から使用可能になった ods animation オプションでアニメーション作成もできます (といっても静止画をコマ送りにするだけですが)。最近、 『可視化で理解する中心極限定理 #rstatsj - qiita』という記事で R 言語で中心極限定理のアニメーションが作成されてるのを見かけたので、ods animation を利用して SAS でも似たようなことをしてみました。
まず、中心極限定理の簡単な説明をします。中心極限定理を知っているという方は次のセクションまで読み飛ばしてかまいません。
中心極限定理というのは、「一定の条件をみたすならどんな確率分布のデータでも、その平均値を以下のような形で変換すれば、観測数 (レコード件数) が多いほど標準正規分布に近づくようになる」

というものです。X̄_N というのが実際のデータの算術平均 (このデータには N 件のレコードがあるとします。N 個の値の平均です) で、μ は理論上の平均、 σ は標準偏差です。

とか

のような書き方をしている教科書もありますが、どれも同じ意味です。中心極限定理というものが存在するので、元の分布が分からなくても、(レコード数が多ければ) 正規分布に近づくので、有意検定とか、カイ二乗検定などいろいろな仮説検定が可能になるので、推理統計学 (数理統計学) では非常に重要な理論です。
しかし、分布が近づく、というのは非常にイメージしにくいでしょう。「大数の法則」ならば、レコード数が増えれば算術平均が真の値に収束する、ということでイメージしやすいのですが、中心極限定理は点ではなく分布、つまり「バラつき」がある形に収束する、ということで、いくらか抽象的です。そこで、アニメーションを用いて中心極限定理の視覚的なイメージを与えてみましょう。

SAS コード

まず、SAS でアニメーションを作成する方法は、SAS公式のリファレンスで確認できます (2015年5月12日現在英語のみ)。 日本語で解説したものとしては SAS忘備録さんの 『SASでアニメーションを作る方法(SAS9.4以降)』 などが参考になります。これを元に、少しづつレコードを増やした場合のヒストグラムを次々表示することで、ヒストグラムが標準正規分布の形に近づいていく様子を見せたいと思います。
しかし、あるデータから計算できる平均値は1つの数字だけです。つまり、1つの点でしかないのでヒストグラムで分布を表現できません。そこで、SF ちっくな表現ですが、あたかも平行世界があるように、「あるデータ」を何度も複製して、それぞれの平均をとることで平均値の分布を取ります。つまり、N 件レコードのあるデータを M 回複製する、ということです。もちろん、それぞれのデータの中身が全く同じであっては意味がありませんので、以下のようなイメージで乱数で作成します。

データ1	データ2	...
1.0	1.1	1.2	1.5
2.0	2.1	1.9	2.5
3.0	3.1	3.5	3.0
1	0	1	1

平均値をいくつも計算する必要があるだけなので、実際には N×M 件のレコードがあれば別々のデータセットにする必要はありません。そのため、この後に掲載するコードでは1つのデータセットだけを作成しています。
さて、今回は、どんな分布でも標準正規分布になることを確認したいので、正規分布、カイ二乗分布、ベルヌーイ分布の3種類で試してみます。これら3種類の分布はどれも、理論上の平均値と分散がわかっています (標準偏差の2乗が分散なので、逆に分散の平方根をとれば標準偏差がわかります)。

正規分布は、パラメータとして平均と分散の値を自由に指定することができます。ただし、SAS の関数では、分散の代わりに標準偏差を指定する必要があります。正規分布は左右対称な山なりの分布で、山の頂点がちょうど平均値になります。つまり、平均値で指定した値が発生しやすく、平均値から離れた値ほど発生しにくくなります。平均がゼロ、分散が1の正規分布を特に標準正規分布といいます。中心極限定理は標準正規分布に収束するという話なので、今回は平均が10、分散が25の正規分布にします。

カイ二乗分布は、正の値しかないような分布です。そのため、分布は山なりですが、正規分布のように対称ではありません。k というパラメータで形が決まります。今回は 2 にしました。カイ二乗分布の平均は k、分散は 2k であることがわかっています。
ベルヌーイ分布というのは、 0 か 1 の2通りしか発生しない分布です。たとえば、コイン投げの結果も表か裏の2通りしかないので、ベルヌーイ分布と考えられます。ベルヌーイ分布は、1 になる確率、コイン投げでいうなら表が出る確率 p の値で分布の形が決まります。0.5 にしてもいいのですが、面白くないので 0.01 というとても不公平な値に設定しましょう。これで裏 (ゼロ) が出る確率が 99.9 パーセントになるので、非常にゼロに偏った分布になります。ベルヌーイ分布の平均は p、分散は p(1-p) になることがわかっています。
今回の場合、3つの分布はだいたいこんな感じになります。

(正規分布)

(カイ二乗分布)

(ベルヌーイ分布)

以上を踏まえて、以下のコードでデータを複製します。とりあえず、500件のレコードのデータを500個複製してみました。それぞれ最初のマクロ変数 size と rep で数を調整できますが、大きくし過ぎると計算に時間がかかります。

/* 保存先フォルダ。末尾に\ や / の区切り記号をつけてください*/
%let path =hogehoge;

%let m    = 10;    /* 正規分布の平均 */
%let s    = 5;    /* 正規分布の標準偏差 */
%let k    = 2;    /* カイ2乗分布の自由度 */
%let p    = .01;   /* ベルヌーイ試行の成功確率 */
%let size = 500; /* 最大サンプルサイズ*/
%let rep  = 500; /* 試行回数*/
%let f    = 20;   /* アニメーションのコマ数*/

%let d=%eval(&size./&f.) /* 1コマごとに増やすサンプルサイズ*/;

/* サイズ size のサンプルを rep 回作成 */
data work.random(keep= j normal chisq bernoulli);
    do i=1 to &rep.;
        do j=1 to &size.;
            /* 3種類の乱数を生成 */
            r_normal    = rand('Normal', &m., &s.);
            r_chisq     = rand('ChiSqaure', &k.);
            r_bernoulli = rand('bernoulli', &p.);

            /* 平均を求めるために累積和を計算 */
            c_normal    + r_normal;
            c_chisq     + r_chisq;
            c_bernoulli + r_bernoulli; 
            /* 算術平均 */
            m_normal     = c_normal    / j;
            m_chisq      = c_chisq     / j;
            m_bernoulli  = c_bernoulli / j;
            /* 標準化 */
            normal    = sqrt(j) * (m_normal - &m.) / &s.;
            chisq     = sqrt(j) * (m_chisq - &k.) / sqrt(2*&k.);
            bernoulli = sqrt(j) * (m_bernoulli - &p.) / sqrt(&p.*(1-&p.));
            output;
        end;
        c_normal    = 0;
        c_chisq     = 0;
        c_bernoulli = 0; 
    end;
run;

ループ処理で作成しているので、 j がレコード件数に対応しています。500回複製しているので、レコード j 件のデータの平均値は 500通り存在することになりますので、平均値のヒストグラムが作成できます。
これで必要なデータが出来ましたので、次に以下のようなマクロでヒストグラムのアニメーションを作成します。

%macro clt(var);
    options  animation=start animduration=0.2   printerpath=GIF;
    ods printer file="&path.clt_demo_&var..gif"  ;
    /* 表示範囲をそろえるために最大絶対値を取得*/
    proc summary data=work.random;
        var &var.;
        output out=work.maxinterval min=min max=max;
    run;
    data _null_;
        set work.maxinterval;
        interval = ceil(max(abs(min), max));
        call symput('interval', interval);
        stop;
    run;    
    /* 最初のコマを長く表示*/
    %do k=1 %to 5;
        %hist(1);
    %end;
    /* 各コマのヒストグラムを順に表示*/
    %do i =1 %to &f.;
        %let l = %eval(&i.*&d.);
        %hist(&l.);
    %end;
    /* 最後のコマを長く表示*/
    %do k=1 %to 5;
        %hist(&size.);
    %end;

    ods printer close ;
    options  animation=stop ;
%mend;

/* 各コマのヒストグラムと分布曲線を描くマクロ*/
%macro hist(x);
    title "N=&x.";
    proc univariate data=work.random(where=(j=&x.)) noprint;
        var &var.;
        histogram / midpoints = -&interval. to &interval. by .5
                    vaxis     = 0 to 50 by 10
                    kernel
                    normal (noprint w=1 l=1 mu=0 sigma=1);
        inset mean std;
    run;
%mend;

univariate プロシジャで、さきほど作成したデータセットから、j=x に対応するレコードを抜き出して、 x 件のレコードを持つデータの平均値のヒストグラムを作成しています。また、ヒストグラムに重ねて理論上の分布曲線を描きます。1つは normal というオプションで描かれる標準正規分布、もう1つは kernel オプションで描かれる、ヒストグラムを滑らかにした曲線です。これはカーネル法と呼ばれる方法で作成できる曲線で、そのデータの擬似的な確率分布とみなせます。よって、中心極限定理がうまくはたらいているなら、このカーネル法による曲線が標準正規分布の曲線に重なっていくはずです。また、左上に平均と標準偏差を表示しました。中心極限定理がはたらいているなら、標準正規分布の平均と標準偏差であるゼロと 1 にそれぞれ近づくはずです。
アニメーションの1コマの表示時間は animduration= で指定できます (単位は秒です) が、一律同じなので、特定のコマだけ長く表示、ということは出来ません。そこで、最初と最後のコマは同じヒストグラムを何度も表示することで、長めに表示しているように見せています。
このマクロは、3つの分布でのアニメーションを別々に作成しますので、以下のように3回呼び出して使います。

%clt(normal); /* 正規分布の場合*/
%clt(chisq); /* カイ二乗分布の場合*/
%clt(bernoulli); /* ベルヌーイ分布の場合*/

すべて完了するまで数分かかるかもしれませんが、これで3種類のアニメーションが完成します。まずは正規分布の場合です。青い曲線が標準正規分布で、赤い曲線がカーネル法による曲線です。

元の分布が正規分布なので、最初からあまり代わりません。最後が少し歪んでいますが、まあ乱数なのでそういうこともあります...。次にカイ二乗分布の場合。

こちらは、最初は元のカイ二乗分布のように左に歪んでいますが、すぐに対称な正規分布の形に近づいていきます。
そして、一番おもしろい動きをするのがベルヌーイ分布です。

ということで ベルヌーイ分布であっても中心極限定理が成立することが視覚的にわかりました。
2015/5/15 片桐