確率・統計

(6) 標本分布

前章では、「中心極限定理」によって標本平均が正規分布に従うということ、また「大数の法則」によって標本数を増やすほど標本平均は真の平均に近づいてゆくことを紹介しました。例えば、サイコロを何回も投げた時に 1 の目が出る回数は、試行回数を増やすほど平均が 1 / 6 に近づき、その分布は正規分布に近づいていきます。「大数の法則」は経験上「当然」であるように見えますが、平均と分散を持った確率分布に従う事象ならば標本平均が正規分布に近似できる「中心極限定理」は驚くべき結果であって、正規分布が統計学の中で重要な位置を占めている理由の一つとなっています。「中心極限定理」は、標本平均の誤差が正規分布となることを意味していることにもなるので、正規分布は「誤差の分布」とも呼ばれます。

今回は、標本平均だけではなく、標本分散などの他の統計量の確率分布がどのようになるかを見ていきたいと思います。

1) χ²-分布(Chi-square Distribution)

まずは、確率変数 x が標準正規分布 p(x) = N( 0, 1 ) に従うとき、y = x² がどのような分布に従うのかを求めてみます。まず、求める分布を q(y) としたとき、y = x² の逆関数は二価関数なので、q(y)dy = 2p(x)dx となります。また、dy = 2xdx なので、

になります。

つぎに、独立した確率変数 y_i ( y_i = x_i² ; i = 1,2, ... N ) の和 Σ_i{1→N}( y_i ) の確率分布を考えると、独立な確率分布の和の分布は「畳み込み積分(Convolution)」で表されるので(「(4) 多変数の確率分布」の「3) 確率変数の変換」参照)、N 個の確率変数の和の分布を T_N( y ) とすれば、

が成り立つはずです。実際、N = 2 のときは

で、( u, v ) = ( y₁ + y₂, y₂ ) とすれば、y₁ = u - v, y₂ = v より、ヤコビアンは

となって、

が得られます。u に対する周辺分布が T₂( u ) なので、

となって、畳み込み積分で表すことができます。実際に上式を解いてみると、まず、y₁ = u - v > 0、y₂ = v > 0 より 0 < v < u なので積分の範囲は ( 0, u ) となって、v = usin²θ とすれば dv = 2usinθcosθdθ、v → 0 のとき θ → 0、v → u のとき θ → π / 2 なので、

になります。さらに T₃( y ) は

T₄( y ) は

t = ysin²θ とすれば dt = 2ysinθcosθdθ、t → 0 のとき θ → 0、t → y のとき θ → π / 2 なので、

cos2θ = 2cos²θ - 1 より cos²θ = ( cos2θ + 1 ) / 2 なので、

と求めていくことができます。N に対する T_N( y ) の式を下表にまとめておきます。

上記結果を見ると、e^-y/2 に変化はなく、y のべき乗は N に対して y^{( N - 2 ) / 2} と表せそうです。係数のところが複雑ですが、例えば N が 3, 5, 7 と変化した時に 1 / 3, 1 / 5 と小さくなったり、4, 6, 8 と変化した時には 1 / 4, 1 / 6 と小さくなったりと、変化に一定の法則がありそうです。また、sqrt( 2π ) があるのは N が奇数のみであるという特徴もあります。実際、これらを一般化した式は存在し、次のようになります。

これを、自由度 N の「χ²(カイ二乗)-分布(Chi-square Distribution)」といいます。上記説明から、y = Σ_N{1→N}( x_i² ) は 自由度 N の χ²-分布に従うことになります。

上式の中の Γ( N / 2 ) は「ガンマ関数」と呼ばれる次のような式になります。

ガンマ関数の主な性質として以下のようなものがあります。

• Γ( x + 1 ) = xΓ( x )
• Γ( n + 1 ) = n! ( n が正の整数のとき )
• Γ( 1 / 2 ) = √π

ガンマ関数は階乗 n! を正の実数(実際には実部が正の複素数)まで拡張した関数であり、n が正の整数であれば階乗と一致します。従って、Γ( N / 2 ) は、N が偶数のときは階乗 ( N / 2 - 1 )! 【 ( N / 2 )! ではないことに注意 】になり、奇数の場合は

と計算することができます。

χ²-分布の全範囲の積分値は

なので、確率分布として成り立っていることが分かります(ガンマ関数は常に正なので、T_N( x ) ≥ 0 も成り立っています)。また、累積分布関数 F(x) は

と表すことができます。ここで、γ( N / 2, x / 2 ) は「(第一種)不完全ガンマ関数 ( (Lower) Incomplete Gamma Function )」と呼ばれる以下の形をした関数になります。

ガンマ関数における積分範囲を ( 0, x ] に限定した形であり、逆に [ x, ∞ ) に限定したものは「(第二種)不完全ガンマ関数 ( (Upper) Incomplete Gamma Function )」と呼ばれます。

自由度 N の χ²-分布の積率母関数は

t = ( 1/2 - θ )x とすれば dt = ( 1/2 - θ )dx となって、

g(θ) を微分すると

なので、平均は

になります。もう一度微分して

となるので、E[x²] は

よって、分散 σ² は

になります。

以上まとめると、次のようになります。

χ²-分布のサンプル・プログラムを以下に示します。

/*
  ChiSquareDistribution : カイ二乗分布
*/
class ChiSquareDistribution : public ContDist
{
  unsigned int _n; // 自由度

  // 不完全ガンマ関数の計算
  static double iGamma( unsigned int n, double x );

public:

  /*
    コンストラクタ

    unsigned int n : 自由度
  */
  ChiSquareDistribution( unsigned int n )
    : _n( n ) {}

  // 確率変数 x における確率密度を返す
  double operator[]( double x ) const;

  // 区間 (-∞,a] における確率を返す
  double lower_p( double a ) const;

  double average() const { return( ( _n > 0 ) ? _n : NAN ); }      // 平均値
  double variance() const { return( ( _n > 0 ) ? 2 * _n : NAN ); } // 分散
};

/*
  ChiSquareDistribution::iGamma : 不完全ガンマ関数 Γ( n / 2, x ) の計算

  Γ( n / 2, x ) = ∫{0→x} t^(n/2) e^-t dt

  unsigned int n : 変数 n
  double x : 積分範囲

  戻り値 : 確率密度
*/
double ChiSquareDistribution::iGamma( unsigned int n, double x )
{
  if ( n == 1 ) {
    return( sqrt( M_PI ) * erf( sqrt( x ) ) );
  } else if ( n == 2 ) {
    return( 1.0 - exp( -x ) );
  } else if ( n > 2 ) {
    double a = (double)n / 2.0 - 1.0;
    return( a * iGamma( n - 2, x ) - pow( x, a ) * exp( -x ) );
  }

  return( 0 );
}

/*
  ChiSquareDistribution::operator[] : 確率変数 x における確率密度を返す

  double x : 確率変数

  戻り値 : 確率密度
*/
double ChiSquareDistribution::operator[]( double x ) const
{
  if ( _n == 0 ) return( NAN );
  if ( x == 0 && _n == 1 ) return( INFINITY );
  if ( x < 0 ) return( 0 );

  double n_2 = (double)_n / 2.0;

  return( pow( x, ( (double)_n - 2.0 ) / 2.0 ) * exp( -x / 2.0 ) / ( pow( 2, n_2 ) * tgamma( n_2 ) ) );
}

/*
  ChiSquareDistribution::lower_p : 区間 (-∞,a] における確率を返す

  double a : 区間の上限

  戻り値 : 確率
*/
double ChiSquareDistribution::lower_p( double a ) const
{
  if ( _n == 0 ) return( NAN );
  if ( a <= 0 ) return( 0 );

  return( iGamma( _n, a / 2.0 ) / tgamma( (double)_n / 2.0 ) );
}

χ²-分布の任意範囲の確率を計算するためには不完全ガンマ関数の計算が必要です。不完全ガンマ関数に対して部分積分法を用いると、以下の漸化式が成り立ちます。

χ²-分布で利用する不完全ガンマ関数は α = N / 2 の形を取るので、漸化式を使って α = 1 または 1 / 2 にすることができて、

と求めることができます。ここで、erf( x ) は正規分布のところでも紹介した「誤差関数(Error Function)」で、以下のような式になります。この関数で、標準正規分布の [ -x, x ] の範囲の確率を求めることができるので、ライブラリ関数として標準で用意されています。

不完全ガンマ関数を利用して、任意の値 α に対する [ 0, α ] の範囲の確率を求めることができるので、α < β に対して F(α), F(β) がそれぞれ [ 0, α ], [ 0, β ] の範囲の確率であるとすれば、F(β) - F(α) が [ α, β ] の範囲の確率を表すことになります。

下図は、自由度 N = 1, 2, 5, 10 に対するχ²-分布を示したものです。自由度が大きくなるほど、分布の高さは低く、横に広がった形に変化していきます。また、大きく偏った形から次第に左右対称な形に近づいていきます。

z = ( x - N ) / ( 2N )^1/2 としたとき、dz = dx / ( 2N )^1/2 より q(z) = p(x)dx/dz = ( 2N )^1/2p(x) なので、自由度 N の χ²-分布は

になります。スターリングの公式

は、ガンマ関数にも適用することができて(補足1)、

となるので、これを代入すると

ここで、

になります。また、f(x) = ln( x + 1 ) のマクローリン級数(Maclaurin Series) は

なので、x が充分に小さい場合は x + 1 ≅ exp( x - x² / 2 ) と近似することができます。これを利用すると、N が充分に大きければ

よって、

従って、z = ( x - N ) / ( 2N )^1/2 としたとき、χ²-分布は N が大きくなるに従って標準正規分布に近づいていきます。

平均 μ、分散 σ² を持つ同一確率分布上にある独立した確率変数 x = ( x₁, x₂, ... x_N ) に対して、

としたとき、N が大きくなるに従って、y は標準正規分布に近づいていきます。これを「中心極限定理(Central Limit Theorem)」というのでした。これは自由度 1 の χ²-分布に対しても成り立つので、x = Σ_i{1→N}( x_i ) とすると、μ = 1、σ = √2 より

となります。y が標準正規分布に従うことは先ほど示しましたが、中心極限定理からもこれは明らかであるということになります。

ところで、正規分布 N( μ, σ² ) から N 個の標本 x₁, x₂, ... x_N を抽出して

としたとき、

とすれば、t_i は標準正規分布に従うので、y = Σ_i{1→N}( t_i² ) は自由度 N の χ²-分布に従うことになります。ところが、平均を母平均 μ から標本平均 m に置き換えて、y を

とすると、これは自由度 N - 1 の χ²-分布に従います。x_i - μ は N 個の独立な変数であったのに対し、

の関係式によって、最後の変数は独立ではなくなってしまいます。すなわち、その変数を x_N - m とすれば、

となるので、{ x_i - m } は「線形従属(Linear Dependence)」になります(補足2)。

また、正規分布 N( μ σ² ) からランダムに抽出した N 個の標本を使って標本平均 m を求めると、m はやはり正規分布 N( μ σ² / N ) に従います。よって、標本平均 m を 1 個抽出して

を求めると、y は自由度 1 の χ²-分布に従うことになります。

χ²-分布のその他の性質として「再生性(Reproductive Property)」があります。再生性とは、ある確率分布に従う確率変数の和が、また同じ確率分布に従うことで、例えば正規分布は再生性を持った確率分布です(「(5) 正規分布」の「4) 標本平均と標本分散」参照 )。

独立な確率変数 x_m、x_n がそれぞれ自由度 M、N の χ²-分布に従うとしたとき、x = x_m + x_n の確率分布 p(x) は

になります。ところが、χ²-分布の性質から

が成り立つので、

となって、p(x) = T_M+N( x ) が成り立ちます(補足3)。

以上の、χ²-分布に対する性質をまとめておきます。

2) F-分布(F-Distribution)

二つの独立した確率変数 t, u が、それぞれ自由度 M, N の χ²-分布に従うとき、

がどのような分布になるかを調べてみます。まずは t、u を変数とする同時分布を φ( t, u ) とすると、t、u は独立なので、

で表すことができます。φ( t, u ) の変数 t、u を

と変数変換すると、

となるので、ヤコビアンは

になります。従って、変数変換した分布を f( x, y ) とすれば、

f( x, y ) を y について ( -∞, ∞ ) の範囲で積分したものが x に対する周辺分布となるので、これを f_x( x ) としたとき、

と表すことができます。y = Mu より dy = Mdu なので、

これが、求めたい分布になります。f_x( x ) を f(x) に置き換えて、χ²-分布の式を代入すると次のようになります。

I = ∫{-∞→∞} u^(M+N)/2-1 exp( -u( Mx / N + 1 ) / 2 ) du として、t = u( Mx / N + 1 ) / 2 とすれば dt = du( Mx / N + 1 ) / 2 となって、

よって、

「ベータ関数(Beta Function)」を利用すると、

が成り立つので、ガンマ関数をベータ関数に置き換え、最後に N^(M+N)/2 を分母と分子に掛ければ

になります。この分布を自由度 ( M, N ) の「F-分布(F-distribution)」といいます。但し、x ≤ 0 では χ²-分布 T_N(x) = 0 なので、最初の定義から x ≤ 0 では f(x) = 0 になります。

以下、F-分布を次のように表します。

分布は x > 0 の範囲のみにあるので、F-分布の全事象の積分値は

で計算できます。ここで、y = x / { x + ( N / M ) } とすると、

で、x → 0 のとき y → 0、x → ∞ のとき y → 1 となるので、

∫{0→1} y^M/2-1( 1 - y )^N/2-1 dy = Β( M / 2, N / 2 ) なので積分値は 1 になります。また、ベータ関数は常に正の値を取るので、定義から G_M,N(x) ≥ 0 であり、確率分布として成り立っていることが分かります。もっとも、F-分布は、χ²-分布の積による同時分布から周辺分布を求めた結果なので、確率分布となっていることは明らかです。

上記結果を利用すると、累積分布関数 F(x) は

で表されます。ここで、Β_x( α, β ) は「不完全ベータ関数(Incomplete Beta Function)」を表し、

で定義されます。また、I_x( α, β ) は「正規化不完全ベータ関数(Regularized Incomplete Beta Function)」を表し、

になります。

平均値 μ は

ベータ関数の性質から

さらに、ガンマ関数の性質から

が成り立つので

と得られます。しかし、これは N > 2 の場合に成り立ちます。というのも、

が収束するためには N / 2 - 1 > 0 となる必要があるからです。N ≤ 2 の場合、平均は存在しません。簡単な例として M = N = 2 の場合を考えると

となって、積分値は発散してしまいます(補足4)。

E[x²] も同様に

と計算できるので、分散 σ² は

になります。但し、ここでも N > 4 の場合に成り立つという制約が付きます。

以上まとめると、次のようになります。

F-分布のサンプル・プログラムを以下に示します。

/*
  sinPowerInteg : sin関数のべき乗の積分値計算

  In(x) = ∫{0→arcsin(√x)} sinθ^n dθ

  unsigned int n : べき乗の指数
  double x : 積分範囲
*/
double sinPowerInteg( unsigned int n, double x )
{
  if ( x < 0 || x > 1 ) return( NAN );

  if ( n == 0 ) {
    return( asin( sqrt( x ) ) );
  } else if ( n == 1 ) {
    return( 1.0 - sqrt( 1 - x ) );
  } else {
    double d = -pow( x, ( (double)n - 1.0 ) / 2.0 ) * sqrt( 1 - x ) / (double)n;
    return( d + ( (double)n - 1.0 ) * sinPowerInteg( n - 2, x ) / (double)n );
  }
}

/*
  iBeta : 不完全ベータ関数の計算

  Bx( m / 2, n / 2 ) = ∫{0→x} t^(m/2-1) ( 1 - t )^(n/2-1) dt

  unsigned int m, n : 不完全ベータ関数のパラメータ(a=m/2,b=n/2)
  double x : 積分範囲
*/
double iBeta( unsigned int m, unsigned int n, double x )
{
  if ( x < 0 || x > 1 ) return( NAN );
  if ( m == 0 || n == 0 ) return( NAN );

  double a = (double)m / 2.0;
  double b = (double)n / 2.0;

  if ( n == 1 ) {
    return( sinPowerInteg( m - 1, x ) * 2 );
  } else if ( n == 2 ) {
    return( pow( x, a ) / a );
  } else {
    return( ( b - 1 ) * iBeta( m + 2, n - 2, x ) / a + pow( x, a ) * pow( 1 - x, b - 1 ) / a );
  }
}

/*
  FDistribution : F-分布
*/
class FDistribution : public ContDist
{
  unsigned int _m;  // 自由度 M
  unsigned int _n;  // 自由度 N

public:

  /*
    コンストラクタ

    unsigned int m, n : 自由度
  */
  FDistribution( unsigned int m, unsigned int n )
    : _m( m ), _n( n ) {}

  // 確率変数 x における確率密度を返す
  double operator[]( double x ) const;

  // 区間 (-∞,a] における確率を返す
  double lower_p( double a ) const;

  double average() const { return( ( _n > 2 ) ? _n / ( _n - 2 ) : NAN ); } // 平均値
  double variance() const; // 分散
};

/*
  FDistribution::operator[] : 確率変数 x における確率密度を返す

  double x : 確率変数

  戻り値 : 確率密度
*/
double FDistribution::operator[]( double x ) const
{
  if ( _m == 0 || _n == 0 ) return( NAN );
  if ( x == 0 && _m == 1 ) return( INFINITY );
  if ( x == 0 && _m == 2 )
    return( 2.0 * tgamma( (double)_n / 2.0 + 1 ) / ( (double)_n * tgamma( (double)_n / 2.0 ) ) );
  if ( x <= 0 ) return( 0 );

  double a = (double)_m / 2.0;
  double b = (double)_n / 2.0;
  double s = (double)( _m + _n ) / 2.0;

  double beta = tgamma( a ) * tgamma( b ) / tgamma( s );

  return( pow( _m, a ) * pow( _n, b ) * pow( x, a - 1 ) / ( beta * pow( (double)_m * x + (double)_n, s ) ) );
}

/*
  FDistribution::lower_p : 区間 (-∞,a] における確率を返す

  double a : 区間の上限

  戻り値 : 確率
*/
double FDistribution::lower_p( double a ) const
{
  if ( _m == 0 || _n == 0 ) return( NAN );
  if ( a <= 0 ) return( 0 );

  double beta = tgamma( (double)_m / 2.0 ) * tgamma( (double)_n / 2.0 ) / tgamma( (double)( _m + _n ) / 2.0 );

  if ( _m >= _n )
    return( iBeta( _m, _n, a / ( a + (double)_n / (double)_m ) ) / beta );
  else
    return( 1.0 - iBeta( _n, _m, 1.0 - a / ( a + (double)_n / (double)_m ) ) / beta );
}

/*
  FDistribution::variance : 分散を返す

  戻り値 : 分散
*/
double FDistribution::variance() const
{
  if ( _n <= 4 ) return( NAN );

  double m = _m;
  double n = _n;

  return( 2 * pow( n, 2 ) * ( m + n - 2 ) / ( m * pow( n - 2, 2 ) * ( n - 4 ) ) );
}

F-分布の任意範囲の確率を計算するためには不完全ベータ関数の計算が必要です。不完全ベータ関数に対して部分積分法を用いると、以下の漸化式が成り立ちます。

漸化式によって β = 1 または 1 / 2 にすることができて、

但し、Β_x( α, 1 / 2 ) では t = sin²θ で変数変換して導いています。さらに α = M / 2 ならば

となるので、積分値の部分を I_M-1、y = sin^-1√x とすると、

より、

のような漸化式が得られます。最終的には M - 1 は 0 か 1 になって、

より、不完全ベータ関数を解くことができます。

α < β のとき、α と β を交換して Β_x( β, α ) として計算した方が早く処理ができます。残念ながら、ベータ関数とは異なり、不完全ベータ関数の場合は Β_x( α, β ) ≠ Β_x( β, α ) なので、単純に交換して処理するだけでは正しい結果は得られません。Β_x( α, β ) と Β_x( β, α ) の関係式は次のように得られます。

において、u = 1 - t とすると、du = -dt、t → 0 のとき u → 1、t → x のとき u → 1 - x になるので、

よって、α < β の場合は Β_1-x( β, α ) を計算した上でその結果を Β( α, β ) から減算することで同じ結果が得られます。求めたいのは正規化された不完全ベータ関数 I_x( α, β ) なので、

と計算すればよいことになります。

F-分布のグラフは次のような形になります。

N または M を 1 に固定した上で、もう一方の自由度を変化させた時の F-分布のグラフを下図に示します。

F-分布 ( N = 1 )	F-分布 ( M = 1 )

分散の等しい二つの正規母集団 N( μ₁, σ² ), N( μ₂, σ² ) からそれぞれ大きさ M, N の標本 ( x₁₁, x₁₂, ... x_1M ) と ( x₂₁, x₂₂, ... x_2N ) をランダムに抽出して、標本平均 m₁, m₂ と標本分散 s₁², s₂² を得たとします。このとき、

はそれぞれ、自由度 M - 1, N - 1 の χ²-分布に従うのでした。ここで、N( μ₁, σ² ), N( μ₂, σ² )の分散の不偏推定量をそれぞれ u₁², u₂² としたとき、

になるので、これらの比 y = u₁² / u₂² を計算すると、

となって、y は自由度 ( M - 1, N - 1 ) の F-分布に従うことになります。

次に、上記と同じ条件で抽出した標本の平均と母集団の平均との差 m₁ - μ₁ と m₂ - μ₂ を求めると、これらはそれぞれ正規分布 N( 0, σ² / M ), N( 0, σ² / N ) に従うことになります。正規分布に従う確率変数の和の分布は、それらの分布の平均・分散の和を平均・分散とする正規分布に従うので(「(5) 正規分布」の「4) 標本平均と標本分散」参照)、( m₁ - μ₁ ) - ( m₂ - μ₂ ) は N( 0, ( 1 / M + 1 / N )σ² ) に従い、χ²-分布の性質から

は自由度 1 の χ²-分布に従うことになります。また、Ms₁² / σ², Ns₂² / σ² はそれぞれ自由度 M - 1, N - 1 の χ²-分布に従い、この二つの確率変数の和

は χ²-分布の再生性から自由度 M + N - 2 の χ²-分布に従います。従って、

は自由度 ( 1, M + N - 2 ) の F-分布に従うことになります。

また、正規母集団 N( μ, σ² ) から N 個の標本を抽出して得られた標本平均 m を使って、

と定義すると、χ₁² は自由度 1 の χ²-分布に従うのでした。同じ N 個の標本に対して標本分散 s² を求めると、

は自由度 N - 1 の χ²-分布に従うので、

は自由度 ( 1, N - 1 ) の F-分布に従うことになります。この値を計算してみると、

{ Ns² / ( N - 1 ) } は不偏分散を表すので u² とすれば、

となって、N( m - μ )² / u² が自由度 ( 1, N - 1 ) の F-分布に従うことになります。

最後に、確率変数 x が自由度 ( M, N ) の F-分布に従うとき、1 / x がどのような分布に従うかを調べてみます。y = 1 / x とした時、dy = -dx / x² なので、y が従う分布を p(y) とすれば p(y)dy = G_M,N(x)dx より

よって、1 / x は自由度 ( N, M ) の F-分布に従うことになります。dx / dy の符号を勝手に外してもいいのかということになりますが、例えば y = x と y = -x のどちらで変換しても分布が変わらないのに、後者では分布が負の方向に反転したような形になってしまいます。x の分布 p(x) の微小領域 dx における確率が y の分布 q(y) の微小領域 dy の確率と等しくなるわけなので、ここでは dx と dy のどちらも正でなければならず、絶対値で考える必要があるわけです。

以上の、F-分布に対する性質をまとめておきます。

二項分布 B_N,p(r) に従う確率変数 r において、r が 0 から k までの値を取る確率は Σ_r{0→k} B_N,p(r) で求められます。また、逆に r が k + 1 から N までの値を取る確率は Σ_r{k+1→N} B_N,p(r) を計算すれば得られます。意外なことに、これらの値は F-分布を利用して求めることができてしまいます。

まず、

を部分積分によって展開していくと、

よって、

各項に対して N! / r!( N - r - 1 )! を掛けると

となるので、

が成り立ちます。但し、q = 1 - p としています。これで、二項分布の 0 から r までの和が積分によって計算できることが示されたわけです。上式の左辺の係数 N! / r!( N - r - 1 )! をガンマ関数で表すと

となるので、r + 1 = m / 2、N - r = n / 2 とおくと、N + 1 = ( m + n ) / 2 となって、

次に、積分の部分を x = my / ( my + n ) で変数変換すると、y = nx / m( 1 - x ), dx = { mn / ( my + n )² }dy で x → p のとき y → np / mq、x → 1 のとき y → +∞ になるので、

となって、この式を代入すれば

となります。また、1 - Σ_i{0→r}( B_N,p( i ) ) = Σ_i{r+1→N}( B_N,p( i ) ) なので、

が成り立ちます。

F-分布の "F" は、イギリスの統計学者「ロナルド・フィッシャー(Sir Ronald Aylmer Fisher)」のイニシャルを表しているらしいのですが、本当のところは不明です。また、F-分布は別名「スネデカーの F-分布(Snedecor's F-distribution)」や「フィッシャー - スネデカー分布(Fisher-Snedecor Distribution)」などと呼ばれることもあるようです。"スネデカー"とは統計学者の「ジョージ・スネデカー(George W. Snedecor」のことで、おそらく分布を最初に発表したのがスネデカーではないかと思いますが、本当のところは調べきれませんでした。

3) t-分布(t-Distribution)

今度は、x が標準正規分布 N( 0, 1 ) に従い、y が自由度 N の χ²-分布 T_N(y) に従うとします。x と y が独立ならば、同時分布 p( x, y ) は

y ≤ 0 のときは p( x, y ) = 0 になります。ここで、t = x / ( y / N )^1/2、u = y / N と変数変換すれば、x = t√u、y = Nu となって、ヤコビアン det(J) は

となるので、q( t, u ) = p( x, y )|det(J)| より

但し、y > 0 より u > 0 のときのみ成り立ちます( u ≤ 0 なら y ≤ 0 なので、q( t, u ) = 0 です)。t に対する周辺分布を f_N(t) とすると、

v = ( t² + N )u / 2 と変数変換すると、

なので、

Γ( 1 / 2 ) = √π なので、

となって、

が得られます。これを自由度 N の「t-分布(t-distribution)」といいます。

x は標準正規分布に従うので、x² は自由度 1 の χ²-分布に従うことになります。よって、

とすれば、s は自由度 ( 1, N ) の F-分布に従います。ところが、t = x / ( y / N )^1/2 だったので t² = s という関係式が成り立ちます。従って、t が自由度 N の t-分布に従うならば、s = t² は自由度 ( 1, N ) の F-分布に従うことになります。

分布は t = 0 を中心に左右対象となるため、t-分布の全事象の積分値は

で求められます。u = ( 1 + t² / N )^-1 とすれば

となって、t → 0 のとき u → 1、t → ∞ のとき u → 0 なので

となります。累積分布関数 F(t) は、x < 0 ならば t = -{ N( 1 - u ) / u }^1/2 となることに注意して

ここで、( 1 + x² / N )^-1 = α とすれば

但し、Β_α( N / 2, 1 / 2 ) は不完全ベータ関数、I_α( N / 2, 1 / 2 ) はその正規化されたものを表しています。x ≥ 0 ならば

となりますが、具体的には [ x, ∞ ) を計算する代わりに ( -∞, -x ] を計算して全体からその値を減算しているのと同じことになります。

分布の形状から平均 μ は 0 になります。実際に計算すると、

より、全事象の積分と同様に u = ( 1 + t² / N )^-1 とすれば

なので、平均はゼロになります。しかし、これは N = 1 の時には成り立たず、( -∞, 0 ] の範囲は ∞、[ 0, ∞ ) の範囲は -∞ となって、平均は得られません。N = 1 の場合の t-分布は次のような式になります。

これは「標準コーシー分布」の式であり、平均の存在しない確率分布の代表です。

分散 σ² は次のように計算できます。

ベータ関数をガンマ関数に分解して変形すると、

となるので、

しかし、ここでも N = 1 の場合は当然成り立たず、N = 2 のときは、

より v² = 1 - u とすれば du = -2vdv で、u → 0 のとき v → 1、u → 1 のとき v → 0 なので、

つまり、N = 2 のときは分散は発散してしまいます。

以上まとめると、次のようになります。

t-分布のサンプル・プログラムを以下に示します。

/*
  TDistribution : t-分布
*/
class TDistribution : public ContDist
{
  unsigned int _n;  // 自由度 N

public:

  /*
    コンストラクタ

    unsigned int n : 自由度
  */
  TDistribution( unsigned int n )
    : _n( n ) {}

  // 確率変数 x における確率密度を返す
  double operator[]( double x ) const;

  // 区間 (-∞,a] における確率を返す
  double lower_p( double a ) const;

  double average() const { return( ( _n > 1 ) ? 0 : NAN ); } // 平均値
  double variance() const // 分散
  { return( ( _n > 2 ) ? (double)_n / (double)( _n - 2 ) :
            ( ( _n == 2 ) ? INFINITY : NAN ) ); }
};

/*
  TDistribution::operator[] : 確率変数 x における確率密度を返す

  double x : 確率変数

  戻り値 : 確率密度
*/
double TDistribution::operator[]( double x ) const
{
  if ( _n == 0 ) return( NAN );

  double beta = tgamma( (double)( _n + 1 ) / 2.0 ) / ( tgamma( (double)_n / 2.0 ) * sqrt( M_PI ) );

  return( beta / ( sqrt( _n ) * pow( 1 + pow( x, 2 ) / (double)_n, (double)( _n + 1 ) / 2.0 ) ) );
}

/*
  TDistribution::lower_p : 区間 (-∞,a] における確率を返す

  double a : 区間の上限

  戻り値 : 確率
*/
double TDistribution::lower_p( double a ) const
{
  if ( _n == 0 ) return( NAN );

  double beta = 2 * ( tgamma( (double)_n / 2.0 ) * sqrt( M_PI ) ) / tgamma( (double)( _n + 1 ) / 2.0 );
  double p = iBeta( _n, 1, (double)_n / ( (double)_n + pow( a, 2 ) ) ) / beta;
  if ( a > 0 ) p = 1.0 - p;

  return( p );
}

不完全ベータ関数を求めるところは、F-分布で利用したヘルパ関数の iBeta をそのまま利用しています。

t-分布は下図に示すような分布になります。

見ると正規分布に非常によく似ていることがわかります。実際、Β( 1 / 2, N / 2 ) においてスターリングの公式を利用すると N が十分大きければ

よって、

従って、N が大きい場合は標準正規分布に近似されることになります。

正規母集団 N( μ, σ² ) から大きさ N の標本 x = ( x₁, x₂, ... x_N ) を抽出して

とします。このとき、t は自由度 N - 1 の t-分布に従います。これを証明してみましょう。

m は標本平均で、正規分布 N( μ, σ² / N ) に従うので、( m - μ ) / ( σ / √N ) は標準正規分布に従います。これを y とすれば、

となるので、{ z / ( N - 1 ) }^1/2 = u / σ としたときに z が自由度 N - 1 の χ²-分布に従うことが示されれば、t が自由度 N - 1 の t-分布に従うことが証明されたことになります。この式を z について解くと

となって、χ²-分布の性質から z は自由度 N - 1 の χ²-分布に従います。よって、t は自由度 N - 1 の t-分布に従うことが示されました。

以上の、t-分布に対する性質をまとめておきます。

t-分布は、「ウィリアム・ゴセット(William Sealy Gosset)」によって「平均値の誤差の確率分布(The probable error of a mean)」という論文の中で 1908 年に発表されたのが最初です。当時、ゴセットの務めていたビール醸造会社のギネス社(あのギネス・ブックで有名な会社です)が社員による論文の発表を禁止していたため、ゴセットは「スチューデント(Student)」というペンネームを使って論文を発表していました。そのため、t-分布は「スチューデントのt-分布(Student's t-distribution)」という名でも知られています。

今回は、χ²-分布、F-分布、t-分布の三つの確率分布を中心に紹介しました。どの分布の確率密度関数もかなり複雑でしたが、統計学における推定や検定の問題の基礎となる分布としてどれも重要なものです。統計解析用のソフトウェアは、有名な R をはじめ様々なものがあり、Excel などにも関数として用意されているので、仕事などで使われている方もいらっしゃると思います。様々な統計解析に対してどのようにデータ処理をすればよいか、ある程度はその手法が決まっていてルーチンワークで作業することが多いと思いますが、その基礎となる考え方を理解すると、今までとは違った観点からデータを扱えるようになることも期待できます。そのためにも、今回紹介した分布の意味するところを理解することは重要だと考えています。

補足1) ガンマ関数に対するスターリングの公式

ガンマ関数 Γ( x + 1 ) を次のように表します。

右辺の積分に対して t = x( 1 + u ) で変数変換を行うと、dt = x・du で、t → 0 のとき u → -1、t → ∞ のとき u → ∞ なので、

積分の範囲を [ -1, -ε ), [ -ε, ε ), [ ε, ∞ ) の三つに分割します。但し、ε はゼロに非常に近い任意の正数とします。すると、x が充分に大きければ積分値は [ -ε, ε ) に集中するので、

と近似することができます。u は非常に小さいので、f(u) = ln( 1 + u ) の マクローリン級数(Maclaurin Series) から

と三次以降の項を無視して近似することができます。よって、

積分範囲を ( -∞, ∞ ) に戻していますが、[ -ε, ε ) 以外の値はほとんどゼロで無視できることから、全範囲で積分しても近似できることを意味しています。この積分はガウス積分 ∫{-∞→∞} exp( -u² ) du = √π と変数変換を利用すれば ( 2π / x )^1/2 になるので、

これはスターリングの公式そのものです。

補足2) y = Ns² / σ² が自由度 N - 1 の χ²-分布に従うことの証明

確率変数 x = ( x₁, x₂, ... x_N ) が全て互いに独立で、標準正規分布 N( 0 , 1 ) に従うと仮定します。このとき、x に対する同時分布 p(x) は

になります。任意の直交行列 Q によって、

と変換した時、z = ( z₁, z₂, ... z_N ) は x を回転または鏡映した形になり、そのノルムは変わりません。つまり、

が成り立つことになります。変数を x から z に変換した時、z = Qx を実際に計算すると

より z_i = u_i1x₁ + u_i2x₂ + ... + u_iNx_N なので、∂z_i / ∂x_j = u_ij となって、ヤコビ行列は Q と等しくなります。したがって、ヤコビアンの絶対値は |det( Q )| = 1 になり(「(5) 正規分布」の「補足 4) 行列の積の行列式」参照)、z に対する同時分布が q(z) ならば p(x) = q(z) となって、||z|| = ||x|| より

つまり、z = ( z₁, z₂, ... z_N ) は互いに独立で、標準正規分布 N( 0, 1 ) に従うことを示しています。

直交行列 Q の N 行目の行ベクトルを u_N = ( 1 / √N, 1 / √N, ... 1 / √N ) としても、

なので、Q は直交行列として成り立っています。このとき、

になります。但し、m は標本平均で m = Σ_i{1→N}( x_i / N ) とします。y = Σ_i{1→N}( ( x_i - m )² ) とすれば、

になり、z_i は N( 0, 1 ) に従うので、y は自由度 N - 1 の χ²-分布に従うことになります。

次に、x が正規分布 N( μ, σ² ) に従う時を考えると、

とすれば z_i は N( 0 , 1 ) に従うことになります。m_z = Σ_i{1→N}( z_i / N ) とすれば

但し、m_x = Σ_i{1→N}( x_i / N ) とします。このとき、

なので、y = Σ_i{1→N}( ( x_i - m_x )² ) / σ² とすれば

になって、Σ_i{1→N}( ( z_i - m_z )² ) は自由度 N - 1 の χ²-分布に従うので、y = Σ_i{1→N}( ( x_i - m_x )² ) / σ² も同様に自由度 N - 1 の χ²-分布に従うことになります。

補足3) 畳み込み積分の性質

畳み込み積分(または合成積 ; Convolution)は、画像のサンプル補間の中で「5) 補間関数と畳み込み積分」として紹介しています。入力信号に対する出力が時間に応じて線形的に変化するとき、出力信号の重ね合わせが畳み込み積分の形になり、次のような式で表されます。

f(x) * g(x) のことを ( f * g )(x) と表したり、変数 x を省略して単に f * g と表す場合もあります。

畳み込み積分は、通常の積と同様に以下の性質を持っています。

交換律の場合、畳み込み積分において u = x - t と変数変換すれば t → ∞ のとき u → -∞、t → -∞ のとき u → ∞、du = -dt になるので

と証明することができます。結合律は、

として v = t + u と変数変換すれば、

最後の分配律は、積分の加法性から簡単に示すことができます。以上の性質と χ²-分布の性質

から、T_M * T_N は

となって、これを繰り返すことで T_M * T_N = T_M+N を示すことができます。

補足4) ベータ関数が収束することの証明

ベータ関数 Β( α, β ) = ∫{0→1} t^α-1( 1 - t )^β-1 dt の被積分関数 t^α-1( 1 - t )^β-1 は、α ≥ 1 かつ β ≥ 1 ならば [ 0, 1 ] の範囲で有限なので、積分値も収束します。ところが、0 < α < 1 のときは t → 0 のとき、0 < β < 1 のときは t → 1 のときに、被積分関数は発散します。にもかかわらず、積分値は収束することを示すことができます。

まずは、ベータ関数の積分範囲を二つに分割します。

( 0, 1/2 ] の範囲において、t^α-1( 1 - t )^β-1 dt ≤ Mt^α-1 を満たす定数 M が存在します。また、[ 1/2, 1 ) の範囲において、t^α-1( 1 - t )^β-1 dt ≤ N( 1 - t )^β-1 を満たす定数 N が存在するので、結局 0 < α < 1 で ∫{0→1/2} t^α-1 dt が、0 < β < 1 で ∫{1/2→1} ( 1 - t )^β-1 dt が収束することが証明できればいいことになります。ところが、∫{1/2→1} ( 1 - t )^β-1 dt において 1 - t = u とすれば t → 1/2 のとき u → 1/2、t → 1 のとき u → 0、du = -dt となって、

となるので、結局 ∫{0→1/2} t^α-1 dt が収束することだけを示せばよいことになります。まず、0 < ε < 1/2 として ∫{ε→1/2} t^α-1 dt を計算すると、

となります。よって、∫{0→1/2} t^α-1 dt は収束することが示され、ベータ関数もまた収束することになります。ところが、α < 0 の場合、ε^α は発散してしまうので、ベータ関数も収束しなくなってしまいます。

ガンマ関数が収束することは以前証明しています。ベータ関数とガンマ関数の間には

の関係が成り立っているので、任意の α に対して Γ(α) > 0 であることから、Β( α, β ) が収束することは自明であるともいえます。

1. 「確率・統計入門」 小針あき宏著 (岩波書店)
2. 「統計数学入門」 本間鶴千代著 (森北出版)
3. 「基礎課程 解析入門」 野本久夫/岸正倫 共著 (サイエンス社)
4. Wikipedia

◆◇◆更新履歴◆◇◆

χ²-分布のサンプルプログラムにある iGamma 関数に誤りが見つかったため修正しました(2010-10-24)

F-分布の性質に関する説明で誤りがあったため修正しました(2010-11-07)

(誤)[ ( M + N - 2 ){ ( m₁ - μ₁ ) + ( m₂ - μ₂ ) }² ] / [ ( 1 / M + 1 / N )( Ms₁² + Ns₂² ) ] は自由度 ( 1, M + N - 2 ) の F-分布に従う

(正)[ ( M + N - 2 ){ ( m₁ - m₂ ) - ( μ₁ - μ₂ ) }² ] / [ ( 1 / M + 1 / N )( Ms₁² + Ns₂² ) ] は自由度 ( 1, M + N - 2 ) の F-分布に従う

サンプル・プログラムを少し見直しました(2010-12-19)

前に戻る

タイトルに戻る