確率・統計 (19) ロジスティック回帰 (Logistic Regression Model)

前章で、「一般化線形モデル (Generalized Linear Model)」の概要について紹介をしました。「線形重回帰モデル (Linear Multiple Regression Model)」は一般化線形モデルの一つであり、逆に考えれば、線形重回帰モデルを一般化したものが一般化線形モデルであるともいえます。今回は、一般化線形モデルの代表である「ロジスティック回帰 (Logistic Regression)」について紹介したいと思います。

1) ロジスティックモデル (Logistic Model)

ある試行が成功・失敗するというような事象を考えます。例えば、サイコロを二つ振ったとき、目の数が一致したら成功とするような事象や、コインを投げて表が出たら成功とする事象などです。この時の確率変数は二つのみであり、例えば

と表すことができます。このような変数は「二値確率変数 (Binary Random Variable)」と呼ばれます。試行が成功する確率 P(x=1) が π であるとしたとき、当然 0 ≤ π ≤ 1 であり、試行が失敗する確率 P(x=0) は 1 - π になります。この確率はまとめて

このような試行を N 回繰り返したとき、j 回目 ( 1 ≤ j ≤ N ) の試行に対する確率変数 x_j は 0, 1 のいずれかの値を取ります。j 回目の試行が成功する ( x_j = 1 になる ) 確率を π_j としたとき、同時確率は

と表されます。これは、h(x) = 1, η_j(π) = log( π_j / ( 1 - π_j ) ), T_j(x) = x_j, A(π) = -Σ_j{1→N}( log( 1 - π_j ) ) としたとき N-母数指数型分布族に属します ( 但し、x, π は N 個の x_j, π_j からなるベクトルとします )。もし、全ての j に対して π_j が等しければ、その値を π としたとき上式は

と表すことができます。y = Σ_j{1→N}( x_j ) としたとき、この値は試行が成功した回数を表し、同時確率は

になります。y = Σ_j{1→N}( x_j ) となる x の値は、N 個の要素から y 個を選ぶときの組み合わせ分だけあるので、y を確率変数とした時の周辺確率は、上記同時確率と組み合わせの数 _NC_y の積で求められ

という結果になります。これは、二項分布 B_N,π(y) と同じ確率密度関数です。さらに一般的な場合として、K 個のグループがあって、各グループにおいて n_i 回の試行を行い、成功した回数が y_i 回だとしたとき、その同時確率は

二項分布 B_N,π(y) の期待値は E[y] = Nπ であり、N 回の試行で y 回成功した時の割合は p = y / N なので、その期待値は E[p] = E[y] / N = π です。そこで、確率 π を連結関数 g(π) によって

の形のモデル式に当てはめることを検討します。もし、g(π) = π (恒等関数) とした場合、π は x による線形式で表現できることになりますが、x の定義域が有限でなければ π そのものも有限となりません。特に、π は確率を表すので 0 ≤ π ≤ 1 である必要があり、x は有限の値を取らなければなりません。例えば、

というモデル式に対して、0 ≤ π ≤ 1 を満たすためには 0 ≤ α₀ + α₁x ≤ 1 でなければならないので、

という制限が必要になります。そこで、0 ≤ π ≤ 1 を満たすために最適な連結関数として累積分布関数がよく用いられます。すなわち、任意の t に対して f(t) ≥ 0 を満たす関数 (これは確率密度関数が満たすべき条件の一つです) を使い、

と表します。但し、P(∞) = ∫{-∞→∞} f(t) dt = 1 となる必要があります (確率密度関数はこれも満たす必要があります)。例えば、区間 [ a, b ] 上の一様分布 P_a,b(t) の累積分布関数は

と表されるので、α₀ = -a / ( b - a )、α₁ = 1 / ( b - a ) とすれば π = α₀ + α₁x の形になり、連結関数を恒等関数とした時のモデル式と一致します。しかし前述の通り x の値が制限されるため、このモデル式はほとんど使われないようです。代わりに、正規分布 N( μ, σ² ) を確率密度関数として用いたモデルを考えます。このときの累積分布関数は、

となり、s = ( t - μ ) / σ とすると ds = dt / σ で、t → -∞ のとき s → -∞、t = x のとき s = ( x - μ ) / σ なので、

と表すことができます。但し、Φ(x) は標準正規分布 N( 0, 1 ) の累積分布関数を表します。P(x) = π とすれば、

なので、α₀ = -μ / σ、α₁ = 1 / σ とすれば、これは連結関数 g(π) を標準正規分布の逆累積分布関数 Φ^-1(π) とした一般化線形モデルとなります。さらに一般化して

で表される一般化線形モデルを「プロビット・モデル (Probit Model)」といいます。このとき、π = P(x) は

プロビット・モデルはモデル式に積分を含むので、いわゆる「初等関数 (Elementary Function)」とは異なる関数です。このモデル式によく似た形状を持った初等関数として、「ロジスティック・モデル (Logistic Model)」または「ロジット・モデル (Logit Model)」があります。ロジスティック・モデルの場合、正規分布の代わりに以下のような関数を使います。

s = e^t とすれば、f(s) = s / ( 1 + s )²、ds / dt = e^t = s で、t → -∞ のとき s → 0、t = x^Tα のとき s = exp( x^Tα ) なので、

と求めることができます。P(x) = 1 - 1 / [ 1 + exp( x^Tα ) ] より P(x) は単調増加であり、x^Tα → +∞ のとき P(x) → 1 となることから、f(t) は確率密度関数として成り立っていることになります。P(x) = π とすれば、

となるので、連結関数を g(π) = log( π / ( 1 - π ) ) とすれば一般化線形モデルとなります。この連結関数は「ロジット関数 (Logit Function)」と一般的に呼ばれています。

その他に、「極値分布 (Extreme Distribution)」と呼ばれる確率密度関数の一つである「標準ガンベル分布 (Standard Gumbel Distribution)」を使った以下のようなモデル式もあります。

f(t) = e^t・exp( -e^t ) より s = e^t とおくと、ds / dt = s、f(s) = se^-s、t → -∞ のとき s → 0、t = x^Tα のとき s = exp( x^Tα ) なので、

となります。P(x) は単調増加であり、x^Tα → +∞ のとき P(x) → 1 となることから、極値分布は確率密度関数として成り立ち、P(x) = π とすれば、

となるので、連結関数を g(π) = log( -log( 1 - π ) ) とすれば一般化線形モデルとなります。この連結関数は「Complementary Log-log 関数」といいます。

これまで紹介した累積分布関数は、グラフに表すと全て S 字形の曲線になるため、ギリシャ文字の ς に似た曲線という意味で「シグモイド曲線 (Sigmoid Curve)」と呼ばれます。「シグモイド関数 (Sigmoid Function)」という言葉もありますが、こちらは関数が一意に定義されていて、

となります。これは、ロジスティック・モデルにおける累積分布関数の特殊形 (切片がゼロのモデル) と等しくなります。

まずは、今まで紹介した連結関数用のクラスを以下に示します。

/*
  ProbitModelFunc : プロビット・モデル連結関数
*/
class ProbitModelFunc : public LinkFunction_IF
{
  NormalDistribution norm_;

public:

  ProbitModelFunc()
    : norm_( 0, 1 ) {}

  // 連結関数 g(x)
  virtual double operator()( double x ) const
  {
    if ( x < 0 || x > 1 ) return( NAN );
    if ( x < 0.5 )
      return( -binSearch( norm_, 0.5 - x ) );
    else
      return( binSearch( norm_, x - 0.5 ) );
  }

  // 導関数 g'(x)
  virtual double df( double x ) const
  { return( exp( pow( (*this)( x ), 2 ) * 0.5 ) * sqrt( 2 * M_PI ) ); }

  // 逆関数 g^-1(y)
  virtual double invf( double y ) const
  { return( norm_.lower_p( y ) ); }

  // 属性を表す文字列
  virtual string ident() const
  { return( "Probit Model Function" ); }
};

/*
  LogitFunc : ロジット連結関数
*/
struct LogitFunc : public LinkFunction_IF
{
  // 連結関数 g(x)
  virtual double operator()( double x ) const
  {
    if ( x < 0 || x > 1 ) return( NAN );
    if ( x == 0 ) return( -INFINITY );
    if ( x == 1 ) return( INFINITY );
    return( log( x / ( 1 - x ) ) );
  }

  // 導関数 g'(x)
  virtual double df( double x ) const
  {
    return(
           ( x < 0 || x > 1 ) ? NAN :
           ( ( x == 0 || x == 1 ) ? INFINITY : 1 / ( x * ( 1 - x ) ) )
           );
  }

  // 逆関数 g^-1(y)
  virtual double invf( double y ) const
  { return( 1 - 1 / ( 1 + exp( y ) ) ); }

  // 属性を表す文字列
  virtual string ident() const
  { return( "Logit Function" ); }
};

/*
  LoglogFunc : Complementary Log-log 連結関数
*/
struct LoglogFunc : public LinkFunction_IF
{
  // 連結関数 g(x)
  virtual double operator()( double x ) const
  {
    if ( x < 0 || x > 1 ) return( NAN );
    if ( x == 0 ) return( -INFINITY );
    if ( x == 1 ) return( INFINITY );
    return( log( -log( 1 - x ) ) );
  }

  // 導関数 g'(x)
  virtual double df( double x ) const
  {
    return(
           ( x < 0 || x > 1 ) ? NAN :
           ( x == 0 || x == 1 ) ? INFINITY : 1 / ( ( x - 1 ) * log( 1 - x ) )
           );
  }

  // 逆関数 g^-1(y)
  virtual double invf( double y ) const
  { return( 1 - exp( -exp( y ) ) ); }

  // 属性を表す文字列
  virtual string ident() const
  { return( "Complementary Log-log Function" ); }
};

プロビット・モデルにおいて、連結関数の値のほかに、導関数と逆関数の値を返す処理が必要になります。プロビット・モデルは

に対する逆関数 y = Φ^-1(x) が連結関数なので、標準正規分布 N( 0, 1 ) の区間 -∞ ≤ t ≤ y における確率が x になるような y を「二分探索」や「ニュートン-ラフソン法」などを使って求める必要があります。区間 -∞ ≤ t ≤ y における確率は、連続確率密度関数を表すクラス ContDist 用のメンバ関数 lower_p で求めることができます (*1-1)。また、導関数に対しては、

となるので、やはり確率が x になるような y を求めてからそれを上式に代入して結果を得ます。最後の逆関数は簡単で、lower_p を使って区間 -∞ ≤ t ≤ y における確率を求めるだけです。

ロジスティック・モデルの場合、y = log( x / ( 1 - x ) ) なので連結関数の値を求めるのは簡単です。導関数は、

最後の、Complementary Log-log 関数 y = log( -log( 1 - x ) ) の導関数は、

これらの関数群を、前回作成したスコア法用のサンプル・プログラムに適用すればいいわけですが、一つ問題となるのが π = 0, 1 になる時の処理で、ロジット関数と Complemenmtary Log-log 関数を利用した場合は計算途中で非数 (nan) が発生するため正しい結果が得られません。そのため、0 や 1 に近い値に置き換えて計算する必要があります。プロビット・モデルの場合は計算方式の理由で非数が発生しづらいようです (あくまでサンプル・プログラムの実装での話です)。

二項分布を一般化した場合、対数尤度関数は以下のような式で表されるのでした。

飽和モデルは、パラメータ π を使ったモデル式であり、π_i による偏微分が

と求められることから、∂l / ∂π_i = 0 のとき π_i = y_i / n_i なので、対数尤度の最大値を l( π|y )_max = l( p|y ) としたとき、

です。独立変数 x_i = ( x_i1, x_i2, ... x_ip )^T を使い、g(π_i) = x_i^Tα が成り立つような p 個のパラメータ α = ( α₁, α₂, ... α_p )^T が存在すれば、π_i = g^-1(x_i^Tα) より ∂π_i / ∂α_j = x_ij / g'(π_i) なので、α_j による偏微分 u_j = ∂l / ∂α_j は

となります。任意の指数型分布族 exp( η(θ)T(y) - A(θ) + B(y) ) に対し、u_j = Σ_i{1→K}( ( y_i - μ_i )x_ij / V[y_i]g'(μ_i) ) で表されるのでした。ここでは μ_i = n_iπ_i、V[y_i] = n_iπ_i( 1 - π_i ) であり、g'(μ_i) = ∂g / ∂μ_i = ( ∂g / ∂π_i )( ∂π_i / ∂μ_i ) = g'(π_i) / n_i となるので、両式は一致します。y_i = n_ip_i と表せば ( つまり、実測値を使って確率を計算した結果を p_i とすれば )、上式は

となるので、もし n_i が全て等しければ、その値を N としたとき上式は

となり、連立方程式 u_j = 0 において N は全て除去できるので、ちょうど N = 1 の二項分布 (すなわちベルヌーイ分布) を指数型分布族としたときの方程式と一致します (この場合の最尤推定量は確率そのものになります)。しかし、n_i の偏りが大きい場合は、ベルヌーイ分布を使ってスコア法を適用することはできません。そのため、各試行回数に応じて分布が切り替えられるような仕組みが必要になります。

/*
  ExpFamily_Binomial : 一母数指数型分布族(二項分布)
*/
class ExpFamily_Binomial : public ExpFamily_IF
{
  double n_; // 総試行回数 N

public:

  // コンストラクタ
  ExpFamily_Binomial( unsigned int n )
    : n_( n ) {}

  // A(η)
  virtual double a( double eta ) const
  { return( n_ * log( exp( eta ) + 1 ) ); }

  // 期待値 ( E[y] = A'(η) )
  double average( double eta ) const
  { return( n_ * ( 1.0 - 1.0 / ( exp( eta ) + 1.0 ) ) ); }

  // A(η) の導関数の逆関数 ( η = A'^(-1)(y) )
  virtual double aveInv( double mu ) const
  { return( log( mu / ( n_ - mu ) ) ); }

  // A(η) の二階導関数 = 分散 ( V[y] = A''(η) )
  virtual double variance( double eta ) const
  { return( ( isinf( eta ) ) ? 0 : n_ * exp( eta ) / pow( exp( eta ) + 1, 2 ) ); }

  // 属性を表す文字列
  virtual string ident() const
  { return( "Binomial Distribution" ); }
};

/*
  ExpFamily_MultiBinomial : 一母数指数型分布族(複数パラメータの二項分布)

  各 x の総試行回数に応じて式を切り替える
*/
class ExpFamily_MultiBinomial : public ExpFamily_IF
{
  static const double DEFAULT_H_ = 1E-6; // デフォルトの微小値

  vector<double> n_; // 各 x に対応する二項分布の総試行回数

  // 各関数の利用回数(n_.size()の剰余)
  mutable unsigned int aCnt_;      // a()
  mutable unsigned int aveCnt_;    // average()
  mutable unsigned int aveInvCnt_; // aveInv()
  mutable unsigned int varCnt_;    // variance()

  double h_; // 微小値 ( π = 0, 1 のときに計算結果がゼロになるのを防ぐためのしきい値 )

  // π = 0, 1 のときに計算結果がゼロになるのを防ぐためのチェック
  static double checkProb( double p )
  {
    return(
           ( p <= 0 ) ? DEFAULT_H_ :
           ( ( p >= 1 ) ? 1.0 - DEFAULT_H_ : p )
           );
  }

  // η = eta に対する π の値を返す
  static double eta2Pi( double eta )
  { return( checkProb( 1.0 - 1.0 / ( exp( eta ) + 1.0 ) ) ); }

  // π = pi に対する η の値を返す
  static double pi2Eta( double pi )
  { return( log( pi / ( 1.0 - pi ) ) ); }

public:

  // コンストラクタ
  ExpFamily_MultiBinomial( const vector<unsigned int>& n, double h = DEFAULT_H_ )
    : n_( n.size() ), aCnt_( 0 ), aveCnt_( 0 ), aveInvCnt_( 0 ), varCnt_( 0 ), h_( h )
  {
    if ( h_ <= 0 ) {
      cerr << "Specified h [" << h_ << "] must be greater than zero.";
      cerr << " Changed to default value [" << DEFAULT_H_ << "]" << endl;
      h_ = DEFAULT_H_;
    }
    for ( unsigned int i = 0 ; i < n_.size() ; ++i )
      n_[i] = n[i];
  }

  // A(η)
  virtual double a( double eta ) const
  {
    double ans = n_[aCnt_] * log( exp( eta ) + 1.0 );
    aCnt_ = ( aCnt_ + 1 ) % n_.size();
    return( ans );
  }

  // 期待値 ( E[y] = A'(η) )
  virtual double average( double eta ) const
  {
    double ans = n_[aveCnt_] * eta2Pi( eta );
    aveCnt_ = ( aveCnt_ + 1 ) % n_.size();
    return( ans );
  }

  // A(η) の導関数の逆関数 ( η = A'^(-1)(y) )
  virtual double aveInv( double mu ) const
  {
    double eta = pi2Eta( checkProb( mu / n_[aveInvCnt_] ) );
    aveInvCnt_ = ( aveInvCnt_ + 1 ) % n_.size();
    return( eta );
  }

  // A(η) の二階導関数 = 分散 ( V[y] = A''(η) )
  virtual double variance( double eta ) const
  {
    double p = eta2Pi( eta );
    double ans = n_[varCnt_] * p * ( 1.0 - p );
    varCnt_ = ( varCnt_ + 1 ) % n_.size();
    return( ans );
  }

  // 属性を表す文字列
  virtual string ident() const
  { return( "Binomial Distribution (multiple parameters)" ); }
};

/*
  SigmoidFunc : シグモイド型連結関数

  各 x の総試行回数に応じて式を切り替える
*/
class SigmoidFunc : public LinkFunction_IF
{
  static const double DEFAULT_H_ = 1E-6; // デフォルトの微小値

  const LinkFunction_IF& linkFunc_; // 利用する連結関数
  vector<double> n_; // 各 x に対応する二項分布の総試行回数

  // 各関数の利用回数(n_.size()の剰余)
  mutable unsigned int opCnt_;   // operator()()
  mutable unsigned int dfCnt_;   // df()
  mutable unsigned int invfCnt_; // invf()

  double h_; // 微小値 ( π = 0, 1 のときに計算結果がゼロになるのを防ぐためのしきい値 )

  // checkProb : π = 0, 1 のときに計算結果がゼロになるのを防ぐためのチェック
  double checkProb( double p ) const
  {
    return(
           ( p <= 0 ) ? h_ :
           ( ( p >= 1 ) ? 1.0 - h_ : p )
           );
  }

public:

  /*
    コンストラクタ

    const LinkFunction_IF& linkFunc : 利用する連結関数
    const vector<unsigned int>& n : 各 x の総試行回数
    double h : 微小値
  */
  SigmoidFunc( const LinkFunction_IF& linkFunc, const vector<unsigned int>& n, double h = DEFAULT_H_ )
    : linkFunc_( linkFunc ), n_( n.size() ), opCnt_( 0 ), dfCnt_( 0 ), invfCnt_( 0 ), h_( h )
  {
    if ( h_ <= 0 ) {
      cerr << "Specified h [" << h_ << "] must be greater than zero.";
      cerr << " Changed to default value [" << DEFAULT_H_ << "]" << endl;
      h_ = DEFAULT_H_;
    }
    for ( unsigned int i = 0 ; i < n_.size() ; ++i )
      n_[i] = n[i];
  }

  // 連結関数 g(x)
  virtual double operator()( double x ) const
  {
    double d = linkFunc_( checkProb( x / n_[opCnt_] ) );
    opCnt_ = ( opCnt_ + 1 ) % n_.size();
    return( d );
  }

  // 導関数 g'(x)
  virtual double df( double x ) const
  {
    double d = linkFunc_.df( checkProb( x / n_[dfCnt_] ) ) / n_[dfCnt_];
    dfCnt_ = ( dfCnt_ + 1 ) % n_.size();
    return( d );
  }

  // 逆関数 g^-1(y)
  virtual double invf( double y ) const
  {
    double d = checkProb( linkFunc_.invf( y ) ) * n_[invfCnt_];
    invfCnt_ = ( invfCnt_ + 1 ) % n_.size();
    return( d );
  }

  // 属性を表す文字列
  virtual string ident() const
  { return( linkFunc_.ident() ); }
};

/*
  ScoringMethod_Binomial : 二項分布を確率分布としたスコア法

  連結関数として Probit, Logistic, Complementary Log-log を利用することを想定

  const vector< vector<double> >& x : 独立変数(p個のパラメータのベクトルからなるn個のベクトル)
  const vector<unsigned int>& n : 試行回数
  const vector<unsigned int>& y : 成功した数
  vector<double>& a : 求めた係数
  const LinkFunction_IF& g : 連結関数
  bool verbose : 冗長モード(ON/OFF)
  unsigned int maxCount : 反復処理の最大回数
  double threshold : 収束条件(全係数が threshold 以下なら処理終了)
  double h : 微小値 (π = 0, 1 のときに計算結果がゼロになるのを防ぐためのしきい値)

  戻り値 : 係数が得られた ... true ; データ異常・反復処理回数が最大値を超えた ... false
*/
bool ScoringMethod_Binomial( const vector< vector<double> >& x, const vector<unsigned int>& n, const vector<unsigned int>& y,
               vector<double>& a, const LinkFunction_IF& g,
               bool verbose, unsigned int maxCount, double threshold, double h )
{
  if ( &n == 0 ) {
    cerr << "n not defined." << endl;
    return( false );
  }
  if ( &y == 0 ) {
    cerr << "y not defined." << endl;
    return( false );
  }
  if ( n.size() != y.size() ) {
    cerr << "n size (" << n.size() << ") and y size (" << y.size() << ") not matched." << endl;
    return( false );
  }
  ExpFamily_MultiBinomial pdf( n ); // 各 x 用の二項分布
  SigmoidFunc sf( g, n );           // 各 x 用の連結関数

  vector<double> yDbl( n.size() ); // y をdouble型に置き換える
  for ( unsigned int i = 0 ; i < n.size() ; ++i ) {
    if ( n[i] == 0 ) {
      cerr << "n[" << i << "] is zero ( n must be greater than zero )." << endl;
      return( false );
    }
    if ( y[i] > n[i] ) {
      cerr << "y[" << i << "] (" << y[i] << ") is more than n[" << i << "] (" << n[i] << ") ( y must be less than n )." << endl;
      return( false );
    }
    yDbl[i] = y[i];
    if ( y[i] == 0 ) yDbl[i] = h * (double)( n[i] );
    if ( y[i] == n[i] ) yDbl[i] -= h * (double)( n[i] );
  }

  return( ScoringMethod( x, yDbl, a, pdf, sf, verbose, maxCount, threshold ) );
}

ExpFamily_Binomial は二項分布を指数型分布族としたクラスで、n_i が全て等しければこれをスコア法に用いることができます。n_i が各グループに対して異なる場合は、ExpFamily_MultiBinomial クラスを利用します。このクラスの中で、グループごとの試行回数を保持しておき、スコア法の中では各試行回数を順番に使って計算する流れになりますが、スコア法ではメンバ関数の aveInv と variance を一回の計算の中で各グループに対し一回ずつ、合計で独立変数の数 n だけ呼び出しているため、各メンバ関数が呼び出された時にカウンタを 1 ずつ増加させることで切り替えを行っています。スコア法の実装に依存したやり方でありいい方法とは言えませんが、今回は妥協しました。なお、スコア法が反復処理であるため、一回のループでカウンタをリセットする必要があり、n の剰余をカウンタに反映する形で実現しています。

SigmoidFunc は、連結関数オブジェクト ProbitModelFunc, LogitFunc, LoglogFunc を利用して、試行回数を切り替えながらスコア法用のパラメータを求めるためのクラスで。連結関数オブジェクトのメンバ関数も、operator(), df, invf がスコア法の中でそれぞれ一回ずつ呼び出されているので、切り替えの原理は ExpFamily_MultiBinomial クラスと全く同じです。
ProbitModelFunc, LogitFunc, LoglogFunc では、引数や戻り値が π_i であることを前提としています。しかし、スコア法の中では μ_i = n_iπ_i が利用されるため、SigmoidFunc の中で補正を行なっています。具体的には、operator() と df の引数が μ_i であるため、n_i で除算して連結関数へ渡すようにし、さらに

なので、invf で求めた値に n_i を掛ければ μ_i を得ることができます。

ScoringMethod_Binomial 関数は、ExpFamily_MultiBinomial と SigmoidFunc を使って、二項分布を指数型分布族としたスコア法を適用するための専用関数で、内部ではスコア法をそのまま利用しています。スコア法とは異なり、各グループの試行回数 n と成功した回数 y をそのまま別々に渡し、n を使って ExpFamily_MultiBinomial と SigmoidFunc を初期化した上で、それらを使ってスコア法を行います。

サンプル・プログラムを利用して、文献にあったデータの回帰係数を計算してみたいと思います。下記データは、二硫化炭素 CS₂ に 5 時間暴露されたカブトムシの死亡数を、使用したガス濃度ごとに示したものです。

図表 1-4. 二硫化炭素ガス濃度に対するカブトムシの死亡率
i	CS₂ ガス濃度 x_i (log₁₀mgl^-1)	カブトムシの数 n_i	死亡数 y_i	死亡率 p_i=y_i/n_i
1	1.6907	59	6	0.1017
2	1.7242	60	13	0.2167
3	1.7552	62	18	0.2903
4	1.7842	56	28	0.5000
5	1.8113	63	52	0.8254
6	1.8369	59	53	0.8983
7	1.8610	62	61	0.9839
8	1.8839	60	60	1.0000

このデータから、まずは独立変数を CS₂ ガス濃度 x_i、従属変数を死亡率 p_i = y_i / n_i とし、

の連結関数 g(π) にプロビット・モデル、ロジット関数、Complementary Log-log 関数をそれぞれ適用します。また、π が従う指数型分布族を N = 1 のときの二項分布 (つまりベルヌーイ分布) として回帰係数を求めると、結果は以下のようになります (但し、p₈ = 1 になるため、1 に近い値として 0.999999 に置き換えて計算します)。

図表 1-5. 実測値と当てはめ値の比較 (ベルヌーイ分布)
ガス濃度	カブトムシの数	実測値	Probit	Logit	Loglog
死亡率	死亡数	死亡率	死亡数	死亡率	死亡数	死亡率	死亡数
1.6907	59	0.1017	6	0.0573	3.38	0.0590	3.48	0.0947	5.59
1.7242	60	0.2167	13	0.1789	10.74	0.1642	9.85	0.1878	11.27
1.7552	62	0.2903	18	0.3782	23.45	0.3614	22.41	0.3373	20.92
1.7842	56	0.5000	28	0.6024	33.74	0.6035	33.80	0.5413	30.31
1.8113	63	0.8254	52	0.7859	49.51	0.7933	49.98	0.7571	47.70
1.8369	59	0.8983	53	0.9024	53.24	0.9019	53.21	0.9168	54.09
1.8610	62	0.9839	61	0.9615	59.62	0.9544	59.17	0.9854	61.10
1.8839	60	1.0000	60	0.9868	59.21	0.9786	58.72	0.9991	59.95

グラフを見る限り、Complementary Log-log 関数を使った場合が最も適合しているようです。次に、ScoringMethod_Binomial 関数を各連結関数と組み合わせて処理を行った結果を以下に示します。

図表 1-6. 実測値と当てはめ値の比較 (グループ単位の二項分布)
ガス濃度	カブトムシの数	実測値	Probit	Logit	Loglog
死亡率	死亡数	死亡率	死亡数	死亡率	死亡数	死亡率	死亡数
1.6907	59	0.1017	6	0.0569	3.36	0.0586	3.4575	0.0947	5.59
1.7242	60	0.2167	13	0.1787	10.72	0.1640	9.84	0.1880	11.28
1.7552	62	0.2903	18	0.3787	23.48	0.3621	22.45	0.3380	20.95
1.7842	56	0.5000	28	0.6038	33.82	0.6053	33.90	0.5423	30.37
1.8113	63	0.8254	52	0.7875	49.62	0.7952	50.10	0.7584	47.78
1.8369	59	0.8983	53	0.9037	53.32	0.9032	53.29	0.9177	54.14
1.8610	62	0.9839	61	0.9623	59.66	0.9552	59.22	0.9857	61.11
1.8839	60	1.0000	60	0.9871	59.23	0.9790	58.74	0.9991	59.95

ベルヌーイ分布を利用した場合と比較すると非常に近い値となっています。しかし、後半部分のデータを n_i, y_i ともに 1000 倍すると、ScoringMethod_Binomial 関数を利用した場合、結果は以下のように変化します。

図表 1-7. 実測値と当てはめ値の比較 (グループ単位の二項分布)
ガス濃度	カブトムシの数	実測値	Probit	Logit	Loglog
死亡率	死亡数	死亡率	死亡数	死亡率	死亡数	死亡率	死亡数
1.6907	59	0.1017	6	0.0282	1.67	0.0191	1.12	0.1635	9.65
1.7242	60	0.2167	13	0.1283	7.70	0.0802	4.81	0.2804	16.82
1.7552	62	0.2903	18	0.3374	20.92	0.2593	16.08	0.4398	27.27
1.7842	56	0.5000	28	0.5984	33.51	0.5623	31.49	0.6262	35.07
1.8113	6300	0.8254	5200	0.8090	5097	0.8124	5118	0.8009	5045
1.8369	5900	0.8983	5300	0.9285	5478	0.9317	5497	0.9238	5451
1.8610	6200	0.9839	6100	0.9783	6066	0.9757	6049	0.9816	6086
1.8839	6000	1.0000	6000	0.9946	5968	0.9912	5947	0.9977	5986

グラフを見るとはっきりと分かるように、前半部分の当てはめは後半部分に比べて非常に悪くなります。u_j = Σ_i{1→K}( n_i( p_i - π_i )x_ij / g'(π_i)π_i( 1 - π_i ) ) より、n_i が小さい成分については影響度が小さくなることが主な理由と考えられます。元々のデータでは n_i が互いに非常に近い値だったので、その場合はベルヌーイ分布を利用した近似解でも問題はありませんが、そうでないときは一般化した二項分布で計算する必要があることがこの結果からもわかります。

サンプル・プログラムを利用して u_j = 0 ( j = 1, 2, ... p ) を満たす α を求め、y^_i = n_iπ^_i = n_ig^-1(x_i^Tα) と計算すれば、π^_i は π_i の最尤推定量であり、y^_i は y_i に対する当てはめ値なので、これを対数尤度に代入すれば l( α|y ) の最大値 l( α|y )_max = l( a|y ) が得られ、その式は

となって、y_i の実測値と当てはめ値を使えば各モデル式での対数尤度統計量を得ることができます。y_i と n_i - y_i はそれぞれ施行に成功・失敗した回数の実測値であるのに対し、y^_i = n_iπ_i と n_i - y^_i はモデル式から得られた確率密度関数による期待値なので、実測値を o、期待値を e で表せば、上式は

の形で表せることになります。ここでの Σ は、各試行における成功・失敗回数それぞれ計 2K 個の和を取るという意味になります。これを、別の式に置き換えます。

この式は、「ピアソン・カイ二乗統計量 (Pearson Chi-squared Statistic)」と呼ばれ、「χ²-検定」でも登場しています。o, e にそれぞれの値を代入すると、

となります。f(s) = s・log( s / t ) に対して t のまわりのテーラー展開を使って二次項までの近似式を求めると、

となるので、実測値と当てはめ値が非常に近ければ、Χ² 値は、D 値の近似値となります。D が漸近的に自由度 K - p の χ²-分布に従うことから Χ² もそれに従い、D の代わりに Χ² 値を利用することも可能で、D の方が、値の小さな y_i, n_i - y_i の影響をより強く受けるため Χ² 値の方がよい近似値になるそうです。

D 値、Χ² 値それぞれに対して、「残差 (Residual)」が定義されます。まず、D 値に対する残差は「逸脱度残差 (Deviance Residuals)」と呼ばれ、

と定義されます。但し、sign(x) は x の正負を表し、実測値と当てはめ値の大小関係によって決まります。d_i² は D の和の成分そのものなので、D = Σ_i{1→K}( d_i² ) が成り立ちます ( 補足 1 )。また、Χ² 値に対しては「ピアソン残差 (Pearson Residuals)」

で定義され、これも Χ² = Σ_i{1→K}( Χ_i² ) が成り立ちます。試行回数 N が十分大きければ、二項分布 B_N,π(y) に従う確率変数 y に対し、( y - Nπ ) / [ Nπ( 1 - π ) ]^1/2 は漸近的に標準正規分布 N( 0, 1 ) に従います(*1-2)。従って、ピアソン残差 Χ_i は N( 0, 1 ) に漸近的に従い、Χ_i ≈ d_i より逸脱度残差 d_i についてもそれは成り立ちます。よって、それぞれの残差に対して「標準化残差」

を定義することができます。但し、h_i は「てこ比 (Levarage)」を表します。

対数尤度統計量 D は、飽和モデルとの比較を意味するのに対し、全ての π_i が等しいと仮定した場合のモデルである「最小モデル (Minimal Model)」と比較することも行われます。最小モデルは、対数尤度関数において π_i = π (一定) とすればよいので、

なので、∂l / ∂π = 0 のとき π = Σ_i{1→K}( y_i ) / Σ_i{1→K}( n_i ) であり、これを p^ で表せば、l( π|y )_max = l( p^|y ) より

が最小モデルとの対数尤度統計量になります。C は「尤度比カイ二乗統計量 (Likelihood Ratio Chi-squared Statistc)」と呼ばれ、

より第一項は自由度 p の、第二項は自由度 1 の χ²-分布に従うので、もし l( α|y ) - l( π|y ) ≈ 0 ならば、C は自由度 p - 1 の χ²-分布に従います。この場合、最小モデルが対象のモデルとよく当てはまることを意味し、さらには切片以外の係数が全てゼロに非常に近いということになります。逆に、C の値がありえないほど大きいということは、係数がゼロではありえない、つまり傾きが意味のあるものであると判断することができます。

線形重回帰モデルにおいて、飽和モデルと最小モデルの間の対数尤度統計量を D₀、飽和モデルと関心のあるモデル式の間の対数尤度統計量を D としたとき、

が「決定係数」と等しくなることを前章で示しました(*1-3)。これによく似た式として、

を「(マクファデンの) 擬似 R² 値 ( (McFadden's) Pseudo R-Squared)」といい、アメリカの経済学者「ダニエル・マクファデン (Daniel Little McFadden)」によって提唱されています。尤度は 1 以下であることから対数尤度は必ず負値となり、飽和モデルが対数尤度の最大値であったことから考えれば、最小モデルの対数尤度 l( π|y ) は必ず最小となる (つまり、その絶対値は最大となる) ことから、この擬似 R² 値は必ず 1 以下となります。

ここで注意すべき点として、擬似 R² 値の計算において、対数尤度は定数項部分の log_{n_i}C_{y_i} を無視する必要があります。すなわち、この時の尤度は

で表されます。対数尤度統計量は対数尤度の差で表されるので、定数項は打ち消し合って無視することができるのに対し、擬似 R² 値は定数項の有無によって値が変化します。文献や統計ソフトを確認する限り、擬似 R² 値の計算では定数項は除いてあります。しかし、対数尤度の式には定数項は含まれているので、擬似 R² 値の計算以外では定数項も含めて計算をしています。なお、n_i が全て 1 で y_i が 0 か 1 のいずれかである場合は、ベルヌーイ分布に従う Σ_i{1→K}( n_i ) 回の独立試行を繰り返すモデルを表し、このときは定数項はゼロになります。

表 1-8. 各統計量の比較
ガス濃度	カブトムシの数 n_i	実測値 y_i	Probit	Logit	Loglog
死亡率	死亡数	d_i	Χ_i	C_i	d_i	Χ_i	C_i	d_i	Χ_i	C_i
1.6907	59	0.1017	6	1.34	1.48	31.94	1.28	1.41	32.02	0.18	0.18	32.83
1.7242	60	0.2167	13	0.75	0.77	18.55	1.06	1.10	18.27	0.56	0.57	18.67
1.7552	62	0.2903	18	-1.46	-1.44	11.49	-1.20	-1.18	11.85	-0.80	-0.79	12.24
1.7842	56	0.5000	28	-1.57	-1.59	0.03	-1.59	-1.61	-0.01	-0.63	-0.64	1.06
1.8113	63	0.8254	52	0.75	0.73	6.89	0.61	0.59	6.99	1.29	1.24	6.34
1.8369	59	0.8983	53	-0.14	-0.14	12.80	-0.13	-0.13	12.80	-0.52	-0.54	12.67
1.8610	62	0.9839	61	1.00	0.89	25.96	1.25	1.09	25.68	-0.12	-0.12	26.46
1.8839	60	1.0000	60	1.25	0.88	29.38	1.59	1.13	28.88	0.33	0.23	30.10
総計	Σ(n_i)	p^	Σ(y_i)	D	Χ²	C	D	Χ²	C	D	Χ²	C
481	0.6050	291	10.12	9.51	274.08	11.23	10.03	272.97	3.45	3.29	280.76

上の表は、各連結関数を使って求めた当てはめ値を元に、逸脱度残差 d_i、ピアソン残差 Χ_i の他、尤比度カイ二乗統計量の和の成分 C_i を求めた結果です。表の下に総計として、対数尤度統計量 D、ピアソン・カイ二乗統計量 Χ²、尤比度カイ二乗統計量 C を示しています。D, Χ² は Complementary Log-log 関数を使った場合が最も小さくなるため、この関数が最も適合していることが定量的に示されたことになります。これらは自由度 6 のカイ二乗分布に近似的に従うと考えてよいので、上側 5% 点が 12.59 であることから、プロビット・モデルとロジスティック・モデルに対してはあまり当てはまりがよいとは言えません。C は自由度 1 のカイ二乗分布に近似的に従い、求められた値に対する p 値はありえないほど小さいため、死亡率がガス濃度に無関係に全て等しいとする最小モデルは適用できず、傾きは必要であると判断することができます。

表 1-9. 擬似 R² 値の比較
ガス濃度	l_i
Probit	Logit	Loglog	最小モデル
1.6907	-20.30	-20.22	-19.41	-52.24
1.7242	-31.64	-31.92	-31.51	-50.19
1.7552	-38.42	-38.07	-37.67	-49.91
1.7842	-40.05	-40.09	-39.02	-40.08
1.8113	-29.46	-29.36	-30.01	-36.35
1.8369	-19.41	-19.41	-19.54	-32.21
1.8610	-5.62	-5.90	-5.13	-31.58
1.8839	-0.78	-1.27	-0.05	-30.15
計	-185.68	-186.24	-182.34	-322.72
擬似 R² 値	0.4246	0.4229	0.4350

上表は、各連結関数ごとの擬似 R² 値を計算した結果です。l_i は、対数尤度を計算するときの和の各成分を表しており、下側に示した合計は対数尤度そのものを意味します。各連結関数に対する対数尤度と最小モデルの対数尤度を使えば擬似 R² 値が得られ、それは表の最も下側に示されています。但し前述の通り、対数尤度の計算において定数項は含んでいないことに注意して下さい。どの連結関数に対しても R² は 40% を少し超える程度です。

2) ホズマー・レメショウ検定 (Hosmer-Lemeshow Test)

二項分布を利用した回帰分析は、ある独立変数ベクトル x_i に対して二値の観測値が n_i 個あって、試行に成功した個数が y_i としたときに、p_i = y_i / n_i を元にモデル式を決めるという方式をとっていました。しかし、x_i が連続値であれば、全ての試行に対して x_i が異なるというのが通常となって、n_i はほとんど 1 となり、y_i は 0 か 1 のいずれかしかとらなくなります。全ての i に対して n_i = 1 の場合、飽和モデルに対する対数尤度 l( p|y ) の和の成分は必ず 0 になるので l( p|y ) = 0 です。l( p|y ) が χ²-分布に漸近的に従うことを利用していることから分かるように、この場合は D や Χ² 値とカイ二乗値との近似性はあまりよくなく、そのまま用いることはできません。そこで、次のような処理を行います。

グループに分割するときの目安としては、通常 10 個程度のグループにすることが多いようです。また、各グループの試行回数 n_i はできるだけ等しくなるようにします。このようにして求めたピアソン・カイ二乗統計量 Χ² を「ホズマー・レメショウ統計量 (Hosmer-Lemeshow Statistic)」といい、Χ²_HL で表されます。数値実験の結果から、Χ²_HL は自由度を ( グループ数 - 2 ) とするカイ二乗分布に従うことが知られています。

この方法は、二人の生物統計学者「David W. Hosmer」と「Stanley Lemeshow」によって 1980 年に発表されました。

ホズマー・レメショウ検定を行うためのサンプル・プログラムを以下に示します。

/*
  vector用 less 関数オブジェクト
*/
template<class T> struct LessVector
{
  bool operator()( const vector<T>& v1, const vector<T>& v2 )
  { for ( unsigned int i = 0 ; i < v1.size() && i < v2.size() ; ++i )
      if ( v1[i] != v2[i] ) return( v1[i] < v2[i] );
    return( false );
  }
};

/*
  NullCheck : NULLチェック関数

  T& t : 対象の変数へのリファレンス
  string arg : 変数名(出力用)
*/
template<class T> bool NullCheck( T& t, string arg )
{
  if ( &t != 0 ) return( true );

  cerr << arg << " not defined." << endl;
  return( false );
}

/*
  SizeCheck : コンテナクラスのサイズチェック

  T& t : 対象のコンテナクラス
  string arg1 : t の変数名(出力用)
  unsigned int sz : サイズ
  string arg2 : sz の大きさを持つ変数の名称(出力用)
*/
template<class T> bool SizeCheck( T& t, string arg1, unsigned int sz, string arg2 )
{
  if ( t.size() == sz ) return( true );

  cerr << arg1 << " size (" << t.size() << ") and " << arg2 << " size (" << sz << ") not matched." << endl;
  return( false );
}

/*
  CalcPredictiveValue : 予測値の計算

  const vector< vector<double> >& x : 独立変数
  const vector<double>& a : 係数
  vector<double>& y : 計算した予測値
  const LinkFunction_IF& g : 連結関数
*/
void CalcPredictiveValue( const vector< vector<double> >& x, const vector<double>& a, vector<double>& y,
                          const LinkFunction_IF& g )
{
  // NULL のチェック
  if ( ! NullCheck( x, "Independent Variable x" ) ) return;
  if ( ! NullCheck( a, "Coefficient a" ) ) return;
  if ( ! NullCheck( y, "Dependent Variable y" ) ) return;
  if ( ! NullCheck( g, "Link Function g" ) ) return;

  unsigned int n = x.size(); // 独立変数の数
  unsigned int p = a.size(); // パラメータ数

  y.assign( n, 0 );

  for ( unsigned int i = 0 ; i < n ; ++i ) {
    if ( p != x[i].size() ) {
      cerr << "The size of x[" << i << "] is not equal to p ( = " << p << " ). Stop processing." << endl;
      return;
    }
    double d = 0;
    for ( unsigned int j = 0 ; j < p ; ++j )
      d += a[j] * x[i][j];
    if ( &g != 0 )
      y[i] = g.invf( d );
  }
}

/*
  Bin2Cnt : 二値データから計数データに変換する

  vector< vector<double> >& x : 独立変数ベクトル
  const vector<bool>& bin : 二値データ
  vector<unsigned int>& n : 計数データ(試行回数の総数)
  vector<unsigned int>& y : 計数データ(成功回数)
  bool verbose : 冗長モード(ON/OFF)
*/
bool Bin2Cnt( vector< vector<double> >& x, const vector<bool>& bin,
              vector<unsigned int>& n, vector<unsigned int>& y, bool verbose )
{
  // NULL のチェック
  if ( ! NullCheck( x, "Independent Variable x" ) ) return( false );
  if ( ! NullCheck( bin, "Binary Data bin" ) ) return( false );
  if ( ! NullCheck( n, "Total Trial Count n" ) ) return( false );
  if ( ! NullCheck( y, "Success Count y" ) ) return( false );

  unsigned int totalCnt = bin.size(); // データ総数
  if ( totalCnt == 0 ) {
    cerr << "bin has no data." << endl;
    return( false );
  }

  if ( ! SizeCheck( x, "x", totalCnt, "bin" ) ) return( false );

  LessVector<double> lessVec; // ベクトル同士の比較関数オブジェクト
  /*
    mapCnt : 計数用Map
     first ... 独立変数ベクトル x
     second ... 総試行回数 n と 成功回数 y のペア
  */
  typedef map< vector<double>, pair<unsigned int, unsigned int>, LessVector<double> > MapCntType;
  MapCntType mapCnt( lessVec );

  // 計数処理
  for ( unsigned int i = 0 ; i < totalCnt ; ++i ) {
    MapCntType::iterator it = mapCnt.find( x[i] );
    if ( it == mapCnt.end() ) {
      mapCnt[x[i]] = pair<unsigned int, unsigned int>( 1, ( bin[i] ) ? 1 : 0 );
    } else {
      ++( ( it->second ).first );
      if ( bin[i] ) ++( ( it->second ).second );
    }
  }

  // x, n, y に計数結果を登録 (各データはクリアされる)
  x.clear();
  n.clear();
  y.clear();
  for ( MapCntType::iterator it = mapCnt.begin() ;
        it != mapCnt.end() ; ++it ) {
    x.push_back( it->first );
    n.push_back( ( it->second ).first );
    y.push_back( ( it->second ).second );
    if ( verbose ) {
      PrintVector( "x : ", x.back() );
      cout << "n = " << n.back() << " ; y = " << y.back() << endl;
    }
  }

  return( true );
}

/*
  HL_ShowPara : ホズマー・レメショウ・テストでのグループ化結果を出力する

  unsigned int grpNo : グループ番号
  const vector< vector<double> >& x : 独立変数ベクトル
  const vector<double>& pi : 確率の当てはめ値
  multimap<double,unsigned int>::const_iterator s, e : 要素番号保持データの範囲(eは最後の要素の次を指す)
  double so, fo : 成功・失敗回数(観測値)
  double se, fe : 成功・失敗回数(期待値)
*/
void HL_ShowPara( unsigned int grpNo,
                  const vector< vector<double> >& x, const vector<double>& pi,
                  multimap<double,unsigned int>::const_iterator s,
                  multimap<double,unsigned int>::const_iterator e,
                  double so, double fo, double se, double fe )
{
  // イテレータの型定義
  typedef multimap<double,unsigned int>::const_iterator MMapCit;

  cout << "*** Group No. = " << grpNo << " ***" << endl;
  MMapCit back; // 末尾のイテレータを保持する
  for ( MMapCit cit = s ; cit != e ; ++cit ) {
    unsigned int i = cit->second; // 要素番号
    std::ostringstream header;
    header << "x[" << i << "] = ";
    PrintVector( header.str(), x[i] );
    back = cit;
  }
  cout << pi[s->second] << " <= pi <= " << pi[back->second] << endl;
  cout << '\t' << "Obs." << '\t' << "Est." << endl;
  cout << "success\t" << so << '\t' << se << endl;
  cout << "failure\t" << fo << '\t' << fe << endl << endl;
}

/*
  HosmerLemeshowTest : ホズマー・レメショウ・テスト

  const vector< vector<double> >& x : 独立変数(表示のみに使用する)
  const vector<unsigned int>& n : 試行回数
  const vector<unsigned int>& y : y の実測値
  const vector<double>& pi : 確率の当てはめ値
  unsigned int grpCnt : グループ数
  bool verbose : 冗長モード(ON/OFF)

  戻り値 : Χ^2HLが得られた ... true ; 引数ミスなど ... false
*/
bool HosmerLemeshowTest( const vector< vector<double> >& x, const vector<unsigned int>& n,
                         const vector<unsigned int>& y, const vector<double>& pi,
                         unsigned int grpCnt, bool verbose )
{
  // NULL のチェック
  if ( ! NullCheck( x, "Independent Variable x" ) ) return( false );
  if ( ! NullCheck( n, "Total Trial Count n" ) ) return( false );
  if ( ! NullCheck( y, "Success Count y" ) ) return( false );
  if ( ! NullCheck( pi, "Predictive Probability pi" ) ) return( false );

  // 要素数の取得とチェック
  unsigned int sz = x.size();
  if ( ! SizeCheck( n, "Total Trial Count n", sz, "Independent Variable x" ) ) return( false );
  if ( ! SizeCheck( y, "Success Count y", sz, "Independent Variable x" ) ) return( false );
  if ( ! SizeCheck( pi, "Predictive Probability pi", sz, "Independent Variable x" ) ) return( false );

  if ( grpCnt == 0 ) {
    cerr << "group count must be more than zero." << endl;
    return( false );
  }

  // 計数データの合計 totalCnt を求める
  unsigned int totalCnt = 0;
  for ( unsigned int i = 0 ; i < sz ; ++i ) {
    if ( n[i] < y[i] ) {
      cerr << "y[" << i << "] (" << y[i] << ") is more than n[" << i << "] (" << n[i] << ") ( o must be less than n )." << endl;
      return( false );
    }
    totalCnt += n[i];
  }
  if ( totalCnt == 0 ) {
    cerr << "total count must be more than zero." << endl;
    return( false );
  }
  unsigned int grpSz = totalCnt / grpCnt; // 1 グループあたりの要素数の目安

  // pi と要素番号をペアとする map の作成
  typedef multimap<double,unsigned int> MMap;
  MMap mmapIndex;
  for ( unsigned int i = 0 ; i < sz ; ++i )
    mmapIndex.insert( pair<double,unsigned int>( pi[i], i ) );

  double se = 0, so = 0, fe = 0, fo = 0; // グループごとの成功・失敗回数(実測・期待値)
  unsigned int cnt = 0;   // グループの総数
  unsigned int grpNo = 1; // グループ番号
  double chisqHL = 0;     // 求める X^2HL
  MMap::const_iterator cit = mmapIndex.begin();
  MMap::const_iterator preCit = cit; // グループの範囲の開始番号
  while ( cit != mmapIndex.end() ) {
    unsigned int i = cit->second; // 要素番号を取得
    so += y[i];
    fo += n[i] - y[i];
    se += (double)( n[i] ) * pi[i];
    fe += (double)( n[i] ) * ( 1.0 - pi[i] );
    cnt += n[i];
    ++cit;
    if ( ( ( cnt >= grpSz ) && ( grpNo < grpCnt ) ) || // グループの要素数が目安を超えた(かつ最後のグループではない)
         ( ( cnt > 0 ) && cit == mmapIndex.end() ) ) { // 最後のグループ(全データ集計完了)
      if ( verbose ) HL_ShowPara( grpNo, x, pi, preCit, cit, so, fo, se, fe );
      ++grpNo;
      preCit = cit;
      chisqHL += pow( so - se, 2 ) / se;
      chisqHL += pow( fo - fe, 2 ) / fe;
      so = fo = se = fe = 0;
      cnt = 0;
    }
  }
  cout << "Chi Square HL = " << chisqHL << endl;

  return( true );
}

サンプル・プログラムは大きく三つの関数に分かれていて、得られたデータが二値の場合、同じ独立変数ベクトル x を持つデータごとに総試行回数 n と成功回数 y を計算するための Bin2Cnt、求めた係数 a を使って成功確率の当てはめ値 pi を求めるための CalcPredictiveValue、最後に、ホズマー・レメショウ・テストを行うための HosmerLemeshowTest があります。また、ホズマー・レメショウ・テストを行うためには、独立変数ベクトル x・計数データ ( 総試行回数 n と成功回数 y )・連結関数 g を使って回帰係数 a を求めるための ScoringMethod_Binomial が必要になりますが、これは前節でのサンプル・プログラムで紹介しています。ホズマー・レメショウ・テスト自体は HosmerLemeshowTest 関数だけで実行することができて、その他の関数はサポートに用います。独立変数と二値データがある場合、ロジット関数を連結関数としてホズマー・レメショウ・テストを行うのであれば、以下のような流れで処理することになります。

当てはめ値を計算するための関数 CalcPredictiveValue は非常に単純なものですが、様々な場面で利用することができます。実際、前節にてサンプル・データから求められた当てはめ値はこれを利用して計算しています。ここで得られる当てはめ値は確率に対するものなので、n を掛けることで成功確率の当てはめ値が得られます。同様に、

と計算することで成功・失敗回数の実測値と予測値が得られるので、これをいくつかのグループごとに集計すれば、その結果を使って Χ²_HL を求めることができます。グループの分け方としては、各グループの試行回数ができるだけ均等になるようにする方がよいので、引数として渡したグループ数 grpCnt で総試行回数 totalCnt を割って 1 グループあたりの試行回数の目安 grpSz を求め、あるグループに対して試行回数の和が grpSz を超えたら次のグループへ切り替えるようにします。

3) 名義ロジスティック回帰 (Nomial Logistic Regression)

前節において、各独立変数 x に対して二項分布 B_N,π(y) を指数型分布族とし、連結関数 g(π) = x^Tα としてロジット関数などの累積分布関数を利用した一般化線形モデルを紹介しました。各従属変数は成功・失敗のいずれかを表す二値であり、同じグループに対する成功回数が二項分布になることを利用してモデル式を定義したわけですが、「成功」と「失敗」という呼び方は便宜的なもので、先の例では「生存率(死亡率)」を意味することになり、その他に性別(男性・女性)などが考えられます。しかし、この分類が三つ以上に分かれる場合を想定すると、二項分布を利用したモデル式は適用できなくなります。そこで、別の確率分布を利用することを検討してみます。

C 個のカテゴリからなる確率変数があり、それぞれのカテゴリの発生する確率を π_k ( k = 1, 2, ... C ) とします。但し、Σ_k{1→C}( π_k ) = 1 という制約条件が付加されます。これは、C 個のカテゴリの中で必ずどれか一つのみが発生するということを意味します。このとき、合計 N 回の試行に対して各カテゴリが y_k 回発生するときの確率分布は「多項分布 (Multinomial Distribution)」に従い、y = ( y₁, y₂, ... y_C )^T, π = ( π₁, π₂, ... π_C )^T とすれば

より、η( π ) = ( log( π₁ ), log( π₂ ), ... log( π_C ) )、h( y ) = N! / Π_k{1→C}( y_k! ) とすれば上式は

となって、母数を複数とする指数型分布族であることを示すことができます。J = 2 のとき、多項分布は二項分布と等しくなり、このときは一母数指数型分布族となることは前に示した通りです。従って、C 個のカテゴリに対しては (C-1)-母数指数型分布族であることが予想できて、実際にそのようになります。これは、y と π の各要素の和に対する制約条件から、母数の数を一つ減らすことができるためです。

N 個の独立変数に対して π_i = ( π_i1, π_i2, ... π_iC )^T ( i = 1, 2, ... N ) を定義したとき ( 但し Σ_k{1→C}( π_ik ) = 1 )、対数尤度 l は

となります。但し、y_i = ( y_i1, y_i2, ... y_iC ) は i 番目の独立変数に対する各カテゴリの発生回数で、Σ_k{1→C}( y_ik ) = n_i を満たします。飽和モデルに対しては、∂l / ∂π_ik = 0 を満たす π_ik を求めればよいので、

がゼロになるような π_ik を計算すればいいわけですが、π_ik' 全てが独立した変数ではなく制約条件 Σ_k'{1→C}( π_ik' ) = 1 を持つことから、例えば π_i1 = 1 - Σ_k'{2→C}( π_ik' ) と表され、

という結果が得られます。この結果は k = 1 に対しても成り立っているので、k = 1 の場合だけ特別扱いする必要はありません。また、直感的にも理解しやすい結果となっています。

カテゴリが二つの場合、「ロジスティック・モデル」を利用した連結関数は次のような式でした。

π_i2 ≡ π_i はいわば「成功」する確率で、π_i1 ≡ 1 - π_i が「失敗」する確率を意味します。これを三つ以上のカテゴリの場合に拡張して、ρ_ik ≡ π_ik / π_i1 が、連結関数 log( ρ_ik ) = log( π_ik / π_i1 ) によって p 個の独立変数の成分に対して線形で表されるとします。この時、

と表すことができます。逆に、π_i1 = 1 - Σ_k'{2→C}( π_ik' ) より

なので、π_ik は ρ_i2 から ρ_iC までの C - 1 個の変数を持つ関数で、ρ_ik は π_i2 から π_iC までの C - 1 個の変数を持つ関数です。対数尤度 l は

と表されます。ここで dρ_ik / dξ_ik = e^ξ_ik = ρ_ik、∂ξ_ik / ∂α_kj = x_ij より

なので、分母は打ち消されて 1 になり、先ほど得られた式に一致します。つまり、二項分布を利用したロジスティック・モデルは、先ほど求めた多項分布によるモデル式に含まれることがわかります。

で表され、H の jr 行 jc 列目の要素は ∂u_jr / ∂α_jc = ∂²l / ∂α_jr∂α_jc となるのでした。u_kj は ( C - 1 ) x p 個あるので、α と u( α ) は要素数が ( C - 1 ) x p のベクトルで、H は行列数が ( C - 1 ) x p の対称行列です。∂u_kj / ∂α_k'j' を直接計算してみると、k ≠ k' ならば

となります。多項分布において、n_iπ_ikπ_ik' は共分散を、また n_iπ_ik( 1 - π_ik ) は分散を表します。

とし、H の行と列を、k が等しいものどうしを固める形で構成したとき ( つまり、H の ( k - 2 ) x p + j 行 ( k' - 2 ) x p + j' 列の要素を ∂u_kj / ∂α_k'j' = ∂²l / ∂u_kj∂u_k'j' としたとき )、H を ( C - 1 ) x ( C - 1 ) 個の p x p 部分行列からなる分割行列 ( 行列をいくつかの部分行列に区切って表現した行列 ) とみれば、対角成分にあたる部分行列は k = k' となるので後者の成分から構成され、それ以外の部分行列は前者の成分になります。x'_j = ( x_1j, x_2j, ... x_Nj )^T とし、

を対角要素とする対角行列を W_kk' とすれば、∂u_kj / ∂α_k'j' = x'_j^TW_kk'x'_j' で表されます。H を p x p 部分行列 ( ブロック ) に分けたとき、k 行 k' 列めのブロックは X^TW_kk'X で表され、分割行列は

と変形できて、H の p x ( k - 2 ) + j 行目の行ベクトル h_kj^T は

となり、左辺は H を係数行列とする連立方程式なので、これを解くことを α_kj が収束するまで繰り返します。

三つ以上のカテゴリに分類できる場合、多項分布とロジスティック・モデルを利用して確率の推定を行うことができることがこれまでの流れでわかりました。この手法は「多項ロジスティック回帰 (Multinomial Logistic Regression)」と呼ばれる中の一種で、分類されたカテゴリが「名義尺度 (Nominal Scale, Categorical Scale)」の場合に利用されることから「名義ロジスティック回帰 (Nominal logistic regression)」ともいいます。

名義ロジスティック回帰を行うためのサンプル・プログラムを以下に示します。

/*
  SizeCheck_Loop : コンテナクラスのサイズチェック(要素もコンテナクラスの場合)

  C& c : コンテナを要素とするコンテナクラス
  string arg1 : t の変数名(出力用)
  unsigned int sz : サイズ
  string arg2 : sz の大きさを持つ変数の名称(出力用)
*/
template<class C> bool SizeCheck_Loop( C& c, string arg1, unsigned int sz, string arg2 )
{
  std::ostringstream oss;
  unsigned int i = 0;
  for ( typename C::const_iterator it = c.begin() ; it != c.end() ; ++it ) {
    oss << arg1 << "[" << i << "]";
    if ( it->size() != sz ) {
      cerr << oss.str() << " size (" << it->size() << ") and " << arg2 << " size (" << sz << ") not matched." << endl;
      return( false );
    }
    ++i;
  }
  return( true );
}

/*
  MultinomialLogit_CalcDiagMatrix : 多項分布モデル用対角行列計算

  const vector< vector<double> >& pi : 確率の当てはめ値
  const vector<double>& ni : 各独立変数ごとの総数
  vector<double>& w : 対角行列 W_kk'
  unsigned int c : カテゴリ数
  unsigned int n : 独立変数ベクトルの数
*/
void MultinomialLogit_CalcDiagMatrix( const vector< vector<double> >& pi, const vector<double>& ni, vector<double>& w,
                                      unsigned int c, unsigned int n )
{
  vector<double>::iterator it = w.begin();
  for ( unsigned int rk = 0 ; rk < c - 1 ; ++rk ) {
    // 対角成分にあたる分割行列を先に処理する
    for ( unsigned int i = 0 ; i < n ; ++i )
      *it++ = -pi[i][rk + 1] * ni[i] * ( 1.0 - pi[i][rk + 1] );
    // その他の分割行列を処理
    for ( unsigned int ck = rk + 1 ; ck < c - 1 ; ++ck )
      for ( unsigned int i = 0 ; i < n ; ++i )
        *it++ = pi[i][rk + 1] * pi[i][ck + 1] * ni[i];
  }
}

/*
  MultinomialLogit_CalcCoef : 多項分布モデル用係数計算

  LinearEquationSystem<double>& s : 係数行列を求める対象の連立方程式計算用インスタンス
  const vector<double>& w : 対角行列 W_kk'
  const vector< vector<double> >& x : 独立変数
  unsigned int c : カテゴリ数
  unsigned int p : 独立変数ベクトルの要素数
  unsigned int n : 独立変数ベクトルの数
  unsigned int rk, ck : 分割行列の行列番号
  unsigned int rj, cj : 分割行列内の計算対象要素の行列番号
*/
void MultinomialLogit_CalcCoef( LinearEquationSystem<double>& s,
                                const vector<double>& w, const vector< vector<double> >& x,
                                unsigned int c, unsigned int p, unsigned int n,
                                unsigned int rk, unsigned int rj, unsigned int ck, unsigned int cj )
{
  s[rk * p + rj][ck * p + cj] = 0;

  vector<double>::const_iterator cit =
    w.begin() + n * ( ( c * ( c - 1 ) - ( c - rk - 1 ) * ( c - rk ) ) / 2 + ck - rk );

  for ( unsigned int i = 0 ; i < n ; ++i )
    s[rk * p + rj][ck * p + cj] += x[i][rj] * x[i][cj] * *( cit + i );
}

/*
  Multinomial_CalcCoefMatrix : 多項分布モデル用係数行列の計算

  LinearEquationSystem<double>& s : 係数行列を求める対象の連立方程式計算用インスタンス
  const vector<double>& w : 対角行列 W_kk'
  const vector< vector<double> >& x : 独立変数
  unsigned int c : カテゴリ数
  unsigned int p : 独立変数ベクトルの要素数
  unsigned int n : 独立変数ベクトルの数
*/
void MultinomialLogit_CalcCoefMatrix( LinearEquationSystem<double>& s,
                                      const vector<double>& w, const vector< vector<double> >& x,
                                      unsigned int c, unsigned int p, unsigned int n )
{
  for ( unsigned int rk = 0 ; rk < c - 1 ; ++rk ) {
    for ( unsigned int rj = 0 ; rj < p ; ++rj ) {
      // 分割行列もその要素も対角成分
      MultinomialLogit_CalcCoef( s, w, x, c, p, n, rk, rj, rk, rj );
      // 分割行列が対角成分
      for ( unsigned int cj = rj + 1 ; cj < p ; ++cj ) {
        MultinomialLogit_CalcCoef( s, w, x, c, p, n, rk, rj, rk, cj );
        s[rk * p + cj][rk * p + rj] = s[rk * p + rj][rk * p + cj];
      }
    }
    for ( unsigned int ck = rk + 1 ; ck < c - 1 ; ++ck ) {
      for ( unsigned int rj = 0 ; rj < p ; ++rj ) {
        // 分割行列は非対角で要素は対角成分
        MultinomialLogit_CalcCoef( s, w, x, c, p, n, rk, rj, ck, rj );
        s[ck * p + rj][rk * p + rj] = s[rk * p + rj][ck * p + rj];
        // 分割行列もその要素も非対角
        for ( unsigned int cj = rj + 1 ; cj < p ; ++cj ) {
          MultinomialLogit_CalcCoef( s, w, x, c, p, n, rk, rj, ck, cj );
          s[ck * p + cj][rk * p + rj] = s[ck * p + rj][rk * p + cj] =
            s[rk * p + cj][ck * p + rj] = s[rk * p + rj][ck * p + cj];
        }
      }
    }
  }
}

/*
  Multinomial_CalcUkj : 多項分布モデル用 u_kj の計算

  const vector< vector<double> >& x : 独立変数
  const vector< vector<double> >& y : 従属変数(各カテゴリの発生回数)
  const vector<double>& ni : 各独立変数の総数 ( n )
  const vector< vector<double> >& pi : 確率の当てはめ値
  unsigned int k : カテゴリ番号
  unsigned int j : 独立変数の番号

  戻り値 : 計算した u_kj の値
*/
double MultinomialLogit_CalcUkj( const vector< vector<double> >& x, const vector< vector<double> >& y,
                                 const vector<double>& ni, const vector< vector<double> >& pi,
                                 unsigned int k, unsigned int j )
{
  unsigned int n = x.size();
  vector<double> d( n );

  for ( unsigned int i = 0 ; i < n ; ++i )
    d[i] = ( y[i][k + 1] - ni[i] * pi[i][k + 1] ) * x[i][j];

  return( sum( d ) );
}

/*
  MultinomialLogit_CalcRSide : 多項分布モデル用連立方程式の右辺の計算

  LinearEquationSystem<double>& s : 右辺を求める対象の連立方程式計算用インスタンス
  const vector<double>& w : 対角行列 W_kk'
  const vector< vector<double> >& x : 独立変数
  const vector< vector<double> >& y : 従属変数(各カテゴリの発生回数)
  const vector<double>& ni : 各独立変数ごとの総数
  const vector< vector<double> >& pi : 確率の当てはめ値
  unsigned int c : カテゴリ数
  unsigned int p : 独立変数ベクトルの要素数
  unsigned int n : 独立変数ベクトルの数
*/
void MultinomialLogit_CalcRSide( LinearEquationSystem<double>& s, const vector<double>& w,
                                   const vector< vector<double> >& x, const vector< vector<double> >& y,
                                   const vector<double>& ni, const vector< vector<double> >& pi,
                                   unsigned int c, unsigned int p, unsigned int n )
{
  for ( unsigned int k = 0 ; k < c - 1 ; ++k ) {
    for ( unsigned int j = 0 ; j < p ; ++j ) {
      s.ans( k * p + j ) = 0;
      for ( unsigned int k2 = 0 ; k2 < c - 1 ; ++k2 ) {
        unsigned int rk = ( k > k2 ) ? k2 : k; // 分割行列の行番号
        unsigned int ck = ( k > k2 ) ? k : k2; // 分割行列の列番号
        vector<double>::const_iterator cit =
          w.begin() + n * ( ( c * ( c - 1 ) - ( c - rk - 1 ) * ( c - rk ) ) / 2 + ck - rk );
        for ( unsigned int i = 0 ; i < n ; ++i )
          s.ans( k * p + j ) += x[i][j] * log( pi[i][k2 + 1] / pi[i][0] ) * *( cit + i );
      }
      s.ans( k * p + j ) -= MultinomialLogit_CalcUkj( x, y, ni, pi, k, j );
    }
  }
}

/*
  MultinomialLogistic : 名義ロジスティック回帰

  const vector< vector<double> >& x : 独立変数
  const vector< vector<double> >& y : 従属変数(各カテゴリの発生回数)
  vector< vector<double> >& a : 求めた係数
  bool verbose : 冗長モード(ON/OFF)
  unsigned int maxCount : 反復処理の最大回数
  double threshold : 収束条件(全係数が threshold 以下なら処理終了)

  n は独立変数ベクトルの数
  p は独立変数ベクトルの要素数
  c はカテゴリの数

  x は p 個のパラメータのベクトルからなる n 個のベクトル
  y は c 個のパラメータのベクトルからなる n 個のベクトル
  a は p 個のパラメータのベクトルからなる c - 1 個のベクトル

  戻り値 : 係数が得られた ... true ; データ異常・反復処理回数が最大値を超えた ... false
*/
bool MultinomialLogistic( const vector< vector<double> >& x, const vector< vector<double> >& y,
                          vector< vector<double> >& a,
                          bool verbose, unsigned int maxCount, double threshold )
{
  cout << "*** Multinomial Logistic Regression ***" << endl << endl;

  // NULL のチェック
  if ( ! NullCheck( x, "Independent Variable x" ) ) return( false );
  if ( ! NullCheck( y, "Occurrence Count y" ) ) return( false );
  if ( ! NullCheck( a, "Coefficient a" ) ) return( false );

  unsigned int n = x.size(); // 独立変数ベクトル x_i の数
  if ( n == 0 ) {
    cerr << "x has no data." << endl;
    return( false );
  }
  if ( ! SizeCheck( y, "Occurrence Count y", n, "Independent Variable x" ) )
    return( false );

  unsigned int p = x[0].size(); // 独立変数ベクトルの要素数
  if ( ! SizeCheck_Loop( x, "Independent Variable x", p, "Independent Variable x[0]" ) )
    return( false );
  if ( p == 0 ) {
    cerr << "The size of x is zero." << endl;
    return( false );
  }
  unsigned int c = y[0].size(); // カテゴリ数
  if ( ! SizeCheck_Loop( y, "Occurrence Count y", c, "Occurrence Count y[0]" ) )
    return( false );
  if ( c == 0 ) {
    cerr << "The size of categories is zero." << endl;
    return( false );
  }

  cout << "N = " << n << " ; p = " << p << " ; c = " << c << endl << endl;

  if ( verbose ) {
    PrintMatrix( "x = ", x );
    cout << endl;
    PrintMatrix( "y = ", y );
    cout << endl;
  }

  vector<double> ni( n ); // 各独立変数の総数 ( n )
  for ( unsigned int i = 0 ; i < n ; ++i )
    ni[i] = sum( y[i] );

  vector< vector<double> > pi( n, vector<double>( c ) ); // 確率の当てはめ値 ( n x c )
  // piは y_ik / ni_i で初期化
  for ( unsigned int i = 0 ; i < n ; ++i )
    for ( unsigned int k = 0 ; k < c ; ++k )
      pi[i][k] = y[i][k] / ni[i];

  LinearEquationSystem<double> s( ( c - 1 ) * p ); // 連立方程式計算用インスタンス
  vector<double> w( n * c * ( c - 1 ) / 2 ); // 対角行列 W_kk'

  // 係数の初期化
  a.resize( c - 1 );
  for ( unsigned int k = 0 ; k < c - 1 ; ++k )
    a[k].assign( p, 0 );

  bool isMatched;   // 収束したか
  unsigned int cnt; // 計算回数
  for ( cnt = 0 ; cnt < maxCount ; ++cnt ) {

    if ( verbose ) {
      cout << "----- cnt = " << cnt << " -----" << endl << endl;
      PrintMatrix( "pi = ", pi );
      cout << endl;
    }

    // 対角行列の計算 W_kk' の計算
    MultinomialLogit_CalcDiagMatrix( pi, ni, w, c, n );
    // 係数行列の計算
    MultinomialLogit_CalcCoefMatrix( s, w, x, c, p, n );
    // 右辺の計算
    MultinomialLogit_CalcRSide( s, w, x, y, ni, pi, c, p, n );

    if ( verbose ) {
      cout << "Equation System :" << endl;
      s.print();
      cout << endl;
    }

    // 連立方程式の計算
    if ( ! GaussianElimination( s ) ) {
      cerr << "Failed to calculate coefficients." << endl;
      return( false );
    }

    // 各係数が収束しているかを確認する
    isMatched = true;
    for ( unsigned int k = 0 ; k < c - 1 ; ++k ) {
      for ( unsigned int j = 0 ; j < p ; ++j ) {
        if ( fabs( a[k][j] - s.ans( k * p + j ) ) >= threshold )
          isMatched = false;
        a[k][j] = s.ans( k * p + j );
      }
    }

    if ( verbose ) {
      for ( unsigned int k = 0 ; k < c - 1 ; ++k ) {
        std::ostringstream oss;
        oss << "Regression equation : y[" << k << "] = ";
        PrintEquation( oss.str(), a[k] );
      }
      cout << endl;
    }

    if ( isMatched ) break;

    // カテゴリごとの確率計算
    vector<double> rho( c - 1 );
    for ( unsigned int i = 0 ; i < n ; ++i ) {
      for ( unsigned int k = 0 ; k < c - 1 ; ++k ) {
        double xi = a[k][0] * x[i][0];
        for ( unsigned int j = 1 ; j < p ; ++j )
          xi += a[k][j] * x[i][j];
        rho[k] = exp( xi );
      }
      pi[i][0] = 1.0 / ( 1.0 + sum( rho ) );
      for ( unsigned int k = 1 ; k < c ; ++k )
        pi[i][k] = rho[k - 1] * pi[i][0];
    }
  }

  if ( cnt < maxCount ) {
    // 係数行列の再計算(=フィッシャー情報行列)
    MultinomialLogit_CalcDiagMatrix( pi, ni, w, c, n );
    MultinomialLogit_CalcCoefMatrix( s, w, x, c, p, n );
    LinearEquationSystem<double> inv( 0 ); // 連立方程式計算用インスタンスの逆行列
    Inverse( s, inv );

    cout << "Estimated regression equation" << endl << endl;
    for ( unsigned int k = 0 ; k < c - 1 ; ++k ) {
      std::ostringstream oss;
      oss << " y[" << k << "] = ";
      PrintEquation( oss.str(), a[k] );
      cout << "variance of a[" << k << "] = ( " << -inv[k * p][k * p];
      for ( unsigned int j = 1 ; j < p ; ++j )
        cout << ", " << -inv[k * p + j][k * p + j];
      cout << " )" << endl << endl;
    }

    cout << "Estimated probability" << endl;
    PrintMatrix( "pi = ", pi );
    cout << endl;
  } else {
    cout << "Failed to estimate regression coefficient" << endl << endl;
  }

  return( cnt < maxCount );
}

メインの関数は MultinomialLogistic で、非常に長いプログラムとなっています。最初の半分程度は引数のチェックや必要な変数の初期化などを行っている部分で、その後の for によるループ処理が反復計算を行っている最も重要な個所になります。
計算処理の最初で部分行列 W_kk' の要素を MultinomialLogit_CalcDiagMatrix を使って求めています。W_kk' は N 個の対角成分を持ち、H 内にある部分行列の個数は ( C - 1 )² 個ありますが、明らかに W_kk' = W_k'k であり、W_kk' のうち対角成分とその右上、または左下の部分行列だけに対して計算をすれば充分です。従って、対角成分を保持するのに必要な要素数は N x ( C - 1 ) x [ ( C - 1 ) + 1 ] / 2 = NC( C - 1 ) / 2 個でよく、計算量も半分程度減らすことができます。なお、計算は左上の部分行列から開始して右方向へ処理した後、その下の行の対角成分からまた右方向へ計算するという形で行い、最後に右下の対角成分にあたる部分行列を処理します。
次に、左辺の係数行列にあたる H の要素を MultinomialLogit_CalcCoefMatrix で計算します。X^TW_kk'X は明らかに対称行列であり、H 内での W_kk' の配置も対称であることから、X^TW_kk'X の対角成分より右上または左下を計算すれば、X^TW_kk'X の要素全てだけでなく、H 内で対称の位置にある部分行列 X^TW_k'kX の全要素も得られたことになります。係数行列の計算では、H の右上側の部分行列に対し、さらに右上の要素のみについて行い、その値を部分行列内の左下だけでなく、H 内での左下の部分行列に対してもコピーしています。これによって、計算量はほぼ 1 / 4 程度になります。
最後に MultinomialLogit_CalcRSide で右辺の計算を行い、連立方程式を解いたら、その解が収束しているかを確認します。なお、連立方程式の解法には「ガウスの消去法 (Gaussian elimination)」を利用しています。

ここでも文献にあったサンプル・データを使って処理を行なってみます。下記データは、車の安全性や装備の嗜好に関するドライバーへの聞き取り調査結果の中で、エアコンとパワーステアリングをどれだけ重視するかを示したものです。

どれだけ重視するか(反応)は 4 段階評価になっていて、「D : 重要でない」「C : あまり重要でない」「B : 重要」「A : 非常に重要」となっています。但し、表の中では「D : 重要でない」と「C : あまり重要でない」を合算した形にしてあります。このデータに対して、性別と年齢を独立変数ベクトルとし、各反応が独立変数ベクトルそれぞれの中で発生する確率を名義ロジスティック回帰で推定してみたいと思います。ここでは、「D : 重要でない」「C : あまり重要でない」を基準としてこの確率との比率を使います。また、独立変数ベクトルは以下のようにします。

表 3-1. エアコンとパワーステアリングの重要性
性別	年齢(歳)	反応
C/D	B	A	計
女性	18-23	26	12	7	45
24-40	9	21	15	45
>40	5	14	41	60
男性	18-23	40	17	8	65
24-40	17	15	12	44
>40	8	15	18	41
計	105	94	101	300

と表され、π_i1 が、基準となる「D : 重要でない」「C : あまり重要でない」の発生確率、π_i2, π_i3 がそれ以外の「B : 重要」と「A : 非常に重要」に対する発生確率になります。また、18 から 23 歳までの女性は x_i1 = x_i2 = x_i3 = 0 であり、これが他の性別・年齢に対する基準となります。

サンプル・プログラムを使って処理を行った結果を以下に示します (出力を冗長にするため引数 verbose を ON にしています)。

*** Multinomial Logistic Regression ***

N = 6 ; p = 4 ; c = 3

x = ( 1, 0, 0, 0 )
    ( 1, 0, 1, 0 )
    ( 1, 0, 0, 1 )
    ( 1, 1, 0, 0 )
    ( 1, 1, 1, 0 )
    ( 1, 1, 0, 1 )

y = ( 26, 12, 7 )
    ( 9, 21, 15 )
    ( 5, 14, 41 )
    ( 40, 17, 8 )
    ( 17, 15, 12 )
    ( 8, 15, 18 )

----- cnt = 0 -----

pi = ( 0.577778, 0.266667, 0.155556 )
     ( 0.2, 0.466667, 0.333333 )
     ( 0.0833333, 0.233333, 0.683333 )
     ( 0.615385, 0.261538, 0.123077 )
     ( 0.386364, 0.340909, 0.272727 )
     ( 0.195122, 0.365854, 0.439024 )

Equation System :
(-62.6857)x0 + (-31.9524)x1 + (-21.0864)x2 + (-20.2455)x3 + (31.2019)x4 + (12.7686)x5 + (11.0909)x6 + (16.152)x7 = 14.0668
(-31.9524)x0 + (-31.9524)x1 + (-9.88636)x2 + (-9.5122)x3 + (12.7686)x4 + (12.7686)x5 + (4.09091)x6 + (6.58537)x7 = 6.5478
(-21.0864)x0 + (-9.88636)x1 + (-21.0864)x2 + (0)x3 + (11.0909)x4 + (4.09091)x5 + (11.0909)x6 + (0)x7 = -6.10144
(-20.2455)x0 + (-9.5122)x1 + (0)x2 + (-20.2455)x3 + (16.152)x4 + (6.58537)x5 + (0)x6 + (16.152)x7 = 8.43913
(31.2019)x0 + (12.7686)x1 + (11.0909)x2 + (16.152)x3 + (-54.7347)x4 + (-25.8402)x5 + (-18.7273)x6 + (-23.0809)x7 = -2.35318
(12.7686)x0 + (12.7686)x1 + (4.09091)x2 + (6.58537)x3 + (-25.8402)x4 + (-25.8402)x5 + (-8.72727)x6 + (-10.0976)x7 = 7.97945
(11.0909)x0 + (4.09091)x1 + (11.0909)x2 + (0)x3 + (-18.7273)x4 + (-8.72727)x5 + (-18.7273)x6 + (0)x7 = 3.35057
(16.152)x0 + (6.58537)x1 + (0)x2 + (16.152)x3 + (-23.0809)x4 + (-10.0976)x5 + (0)x6 + (-23.0809)x7 = -21.5174

Regression equation : y[0] = -0.594729x0 + -0.381807x1 + 1.13426x2 + 1.58749x3
Regression equation : y[1] = -1.03114x0 + -0.803674x1 + 1.46287x2 + 2.9008x3

----- cnt = 1 -----

pi = ( 0.524023, 0.28911, 0.186867 )
     ( 0.23501, 0.40309, 0.3619 )
     ( 0.098186, 0.264972, 0.636842 )
     ( 0.650933, 0.24515, 0.103917 )
     ( 0.349621, 0.409351, 0.241028 )
     ( 0.174038, 0.32061, 0.505352 )

Equation System :
(-63.3591)x0 + (-31.5974)x1 + (-21.4658)x2 + (-20.6163)x3 + (31.7604)x4 + (12.64)x5 + (10.9058)x6 + (16.7676)x7 = 14.4935
(-31.5974)x0 + (-31.5974)x1 + (-10.6384)x2 + (-8.93059)x3 + (12.64)x4 + (12.64)x5 + (4.34126)x6 + (6.64286)x7 = 7.13147
(-21.4658)x0 + (-10.6384)x1 + (-21.4658)x2 + (0)x3 + (10.9058)x4 + (4.34126)x5 + (10.9058)x6 + (0)x7 = -6.14966
(-20.6163)x0 + (-8.93059)x1 + (0)x2 + (-20.6163)x3 + (16.7676)x4 + (6.64286)x5 + (0)x6 + (16.7676)x7 = 8.99699
(31.7604)x0 + (12.64)x1 + (10.9058)x2 + (16.7676)x3 + (-55.4565)x4 + (-24.3506)x5 + (-18.4408)x6 + (-24.1253)x7 = -4.94817
(12.64)x0 + (12.64)x1 + (4.34126)x2 + (6.64286)x3 + (-24.3506)x4 + (-24.3506)x5 + (-8.04907)x6 + (-10.2488)x7 = 6.37968
(10.9058)x0 + (4.34126)x1 + (10.9058)x2 + (0)x3 + (-18.4408)x4 + (-8.04907)x5 + (-18.4408)x6 + (0)x7 = 2.62442
(16.7676)x0 + (6.64286)x1 + (0)x2 + (16.7676)x3 + (-24.1253)x4 + (-10.2488)x5 + (0)x6 + (-24.1253)x7 = -22.8294

Regression equation : y[0] = -0.590797x0 + -0.388125x1 + 1.12827x2 + 1.5877x3
Regression equation : y[1] = -1.03902x0 + -0.813x1 + 1.47805x2 + 2.91668x3

----- cnt = 2 -----

pi = ( 0.524195, 0.290344, 0.185461 )
     ( 0.234584, 0.40153, 0.363886 )
     ( 0.0975796, 0.264427, 0.637993 )
     ( 0.652471, 0.245143, 0.102386 )
     ( 0.350992, 0.407527, 0.241481 )
     ( 0.174276, 0.32035, 0.505374 )

Equation System :
(-63.3346)x0 + (-31.5786)x1 + (-21.4374)x2 + (-20.5971)x3 + (31.7196)x4 + (12.5992)x5 + (10.905)x6 + (16.7599)x7 = 14.5863
(-31.5786)x0 + (-31.5786)x1 + (-10.6237)x2 + (-8.92676)x3 + (12.5992)x4 + (12.5992)x5 + (4.33004)x6 + (6.63776)x7 = 7.17966
(-21.4374)x0 + (-10.6237)x1 + (-21.4374)x2 + (0)x3 + (10.905)x4 + (4.33004)x5 + (10.905)x6 + (0)x7 = -6.13127
(-20.5971)x0 + (-8.92676)x1 + (0)x2 + (-20.5971)x3 + (16.7599)x4 + (6.63776)x5 + (0)x6 + (16.7599)x7 = 9.00441
(31.7196)x0 + (12.5992)x1 + (10.905)x2 + (16.7599)x3 + (-55.3536)x4 + (-24.2819)x5 + (-18.4757)x6 + (-24.1063)x7 = -5.07926
(12.5992)x0 + (12.5992)x1 + (4.33004)x2 + (6.63776)x3 + (-24.2819)x4 + (-24.2819)x5 + (-8.05939)x6 + (-10.2488)x7 = 6.25696
(10.905)x0 + (4.33004)x1 + (10.905)x2 + (0)x3 + (-18.4757)x4 + (-8.05939)x5 + (-18.4757)x6 + (0)x7 = 2.62152
(16.7599)x0 + (6.63776)x1 + (0)x2 + (16.7599)x3 + (-24.1063)x4 + (-10.2488)x5 + (0)x6 + (-24.1063)x7 = -22.7995

Regression equation : y[0] = -0.590799x0 + -0.388128x1 + 1.12827x2 + 1.58771x3
Regression equation : y[1] = -1.03908x0 + -0.813018x1 + 1.47811x2 + 2.91675x3

Estimated regression equation

 y[0] = -0.590799x0 + -0.388128x1 + 1.12827x2 + 1.58771x3
 variance of a[0] = ( 0.0806426, 0.0903072, 0.116722, 0.162328 )

 y[1] = -1.03908x0 + -0.813018x1 + 1.47811x2 + 2.91675x3
 variance of a[1] = ( 0.109228, 0.103064, 0.160738, 0.178863 )

Estimated probability
pi = ( 0.524195, 0.290344, 0.185461 )
     ( 0.234584, 0.40153, 0.363886 )
     ( 0.0975796, 0.264427, 0.637993 )
     ( 0.652471, 0.245143, 0.102386 )
     ( 0.350992, 0.407527, 0.241481 )
     ( 0.174276, 0.32035, 0.505374 )

反復処理は三回で収束し、最後に推定結果が出力されます。推定確率は処理時に計算しているので、後で計算する必要はなくそのまま続けて出力されます。その結果は次のようになりました。

表 3-2. 名義ロジスティック回帰結果
性別	年齢(歳)	反応	計
1	2	3
推定確率	当てはめ値	実測値	ピアソン残差	推定確率	当てはめ値	実測値	ピアソン残差	推定確率	当てはめ値	実測値	ピアソン残差
女性	18-23	0.5242	23.59	26	0.4965	0.2903	13.07	12	-0.2948	0.1855	8.35	7	-0.4658	45
24-40	0.2346	10.56	9	-0.4790	0.4015	18.07	21	0.6896	0.3639	16.37	15	-0.3398	45
>40	0.0976	5.85	5	-0.3533	0.2644	15.87	14	-0.4684	0.6380	38.28	41	0.4397	60
男性	18-23	0.6525	42.41	40	-0.3702	0.2451	15.93	17	0.2670	0.1024	6.66	8	0.5213	65
24-40	0.3510	15.44	17	0.3960	0.4075	17.93	15	-0.6922	0.2415	10.63	12	0.4218	44
>40	0.1743	7.15	8	0.3197	0.3204	13.13	15	0.5148	0.5054	20.72	18	-0.5976	41
二乗和		0.997		1.597		1.333	3.927

反応の部分は 1 が「D : 重要でない」または「C : あまり重要でない」、2 が「B : 重要」、そして 3 が「A : 非常に重要」になります。各反応に対し、サンプル・プログラムで得られた推定確率とそれによる当てはめ値 ( e )、実測値 ( o )、最後にピアソン残差 ( o - e ) / √e を示してあります。下端の「二乗和」行は、ピアソン残差の二乗和を計算した結果を示しています。右端の列にある「計」は各反応の当てはめ値または実測値の合計ですが、最も右下の数値はピアソン残差の二乗和の合計を表しており、これは「ピアソン・カイ二乗統計量 (Pearson Chi-squared Statistic)」と等しくなります。

に推定確率 π_ik を代入することで求められます。但し、π_ik に依存しない項は定数項 K で表しています。最大モデルは π_ik = y_ik / n_i となることは先述した通りです。また最小モデルは、各カテゴリに対する確率が i に依存しない ( すなわち π_ik = π_k ) と考えればよいので、

より ( 但し、Y_k = Σ_i{1→N}( y_ik ) としています )、π_k による偏微分 ∂l / ∂π_k は

で求められます。ここでも π₁ = 1 - Σ_k'{2→C}( π_k' ) であることから

という結果が得られます。つまり、各カテゴリごとの総数を全ての総数で割った値を使えばよいことになります。

名義ロジスティック・モデルを使って求めた推定確率から計算した対数尤度を l_α とすると、サンプル・データにおける値は -25.37 になります。また、最大モデルに対する対数尤度 l_max = -23.40、最小モデルの場合 l_min = -64.29 となるので、対数尤度統計量 D と尤度比カイ二乗統計量 C は

で求められますが、ここでの対数尤度 l_min, l_α は定数項を除く必要があるので(*3-1)、その値は l_min = -329.27, l_α = -290.35 となって、

という結果になります。飽和モデルでのパラメータ数は、性別 X 年齢(三種類) = 6 個の変数それぞれに対して基準カテゴリを除く二種類のカテゴリ(反応)があったので合計 12 個です。名義ロジスティック・モデルでは係数が 8 つあり、最小モデルでは 2 つになります。従って、D は自由度が 12 - 8 = 4、C は自由度が 8 - 2 = 6 の χ²-分布に漸近的に従い、それぞれの p 値は 0.4144, 9.966E-15 となるので、名義ロジスティック・モデルは最小モデルよりもデータに対して有意に適合しており、飽和モデルと比較してもうまくデータを表せていることが示されています。その反面、擬似 R² 値は 0.1182 と低く、このモデルはデータの全変動の 11.82% 程度しか説明できていないことになります。

確率の比 ρ_ik = π_ik / π_i1 は「オッズ (Odds)」と呼ばれます。ギャンブルなどでもよく使われる用語ですが、その場合は勝った時に支払われる量を意味し、日本の競馬などでは払戻金の倍率を表しています。C = 2 の場合、オッズは π_i / ( 1 - π_i ) となり、その対数はロジット関数そのものになります。
二つのオッズ ρ_ik = π_ik / π_i1 と ρ_lk = π_lk / π_l1 の比を「オッズ比(Odds Ratio)」といいます。ρ_ik と ρ_lk の対数が線形関係

となるので、もし α_k がゼロベクトルに近い場合、オッズ比は 1 に近い値になるはずです。例えば、単純な例として x_i = ( 1, 1 )^T、x_l = ( 1, 0 )^T の場合を考えると、

であり、α_k1 = 0 ならば ρ_ik / ρ_lk = 1 です。独立変数の第一成分が定数項、第二成分が要因の有無を表しているとすれば、オッズ比が 1 ならば要因の有無に対して発生確率が影響しないことを表します。また、オッズ比が 1 より大きければ α_k1 > 0 であり要因によって確率は上がる傾向に、逆に 1 より小さければ確率は小さくなる傾向にあることも読み取ることができます。

サンプルデータにおける独立変数は次のような構成になっていました。

この表から、i = 4 の i = 1 に対するオッズ比は性別に対する影響を、また i = 2, 3 の i = 1 に対するオッズ比は年齢に対する影響をそれぞれ表すことがわかります。これらを計算した結果は次のようになります。

表 3-3. 独立変数ベクトルの構成
i	定数項	性別	年齢 24-40	年齢 >40
1	1	0	0	0
2	1	0	1	0
3	1	0	0	1
4	1	1	0	0
5	1	1	1	0
6	1	1	0	1

表 3-4. オッズ比
反応	B	A
性別	0.678	0.444
年齢 24-40	3.090	4.384
年齢 >40	4.892	18.48

性別による影響については、男性の方が女性よりもエアコンやパワーステアリングを重視しない傾向にあることがわかります。また、年齢に関しては、上昇するほど重要視する傾向が強くなることがはっきりと読み取れます。これらの値は定数項以外の係数 α_kj を使って exp( α_kj ) から求めることもできます。各係数の分散はフィッシャー情報行列から得ることができて、係数が正規分布に漸近的に従うことから 95% 信頼区間は平均 ± 1.96 x 標準誤差(S.E.) を計算すれば求められ、次のような結果になります。

性別による差異については、B に対して区間が 1 を含んでおり、A はそうではないものの、最大値は 1 に近くなっています (実際、信頼度を 99% にすれば 1 を含むようになります)。よって、性別による差異についてははっきりとあると断言することはできません。

表 3-5. オッズ比の信頼区間
反応	B	A
性別	[ 0.376, 1.222 ]	[ 0.236, 0.832 ]
年齢 24-40	[ 1.582, 6.037 ]	[ 1.998, 9.621 ]
年齢 >40	[ 2.221, 10.777 ]	[ 8.067, 42.34 ]

*3-1) 多項分布において、N = 1 で y が二値変数 ( 0 または 1 ) からなるベクトルとすれば、P_1,π( y ) = Π_k( π_k ) であり、対数尤度の定数項はゼロになります。二項分布によるモデルにおいて、指数型分布族をベルヌーイ分布と考えた時、定数項が消滅するのと内容は同じです。ロジスティックモデルの場合、指数型分布族は二項分布や多項分布とはせずに、ベルヌーイ分布や、そのカテゴリが三つ以上になった場合の分布と考えるのが通常のようです。

補足1) ロジスティック回帰の逸脱度残差

ロジスティック・モデルにおける逸脱度残差は次のように表されるのでした。

d_i が実数になるためには、平方根の中身がゼロ以上でなければなりません。ロジスティック・モデルの場合は成り立つことが以下のように証明できます。

y^_i を変数とする以下の関数を定義します。y_i は定数としておきます。

なので、f'( y^_i ) = 0 のとき y^_i = y_i となり、f( y_i ) は極値となります。二階導関数は

であり、y_i ≥ 0 かつ n_i - y_i ≥ 0 ならば常に f''( y^_i ) ≥ 0 となるので、f( y_i ) は極小値になります。最後に、f( y_i ) = 0 より f( y^_i ) ≥ 0 が成り立ちます。

となります。但し、N は独立変数の数、C はカテゴリ数で、i 番目の独立変数の k 番目のカテゴリに対し、y_ik が発生回数、π_ik が発生確率を表します。また、n_i は i 番目の独立変数の総発生回数で Σ_k{1→C}( y_ik ) = n_i を満たします。

飽和モデルの場合は π_ik = y_ik / n_i のとき最大対数尤度となるので、逸脱度残差が計算できるとすれば

で計算できますが、残念ながら log( y_ik / n_i ) ≥ log π_ik が常に成り立つわけではないので名義ロジスティック・モデルの場合は逸脱度残差は計算できない場合があります。

Π_j{1→N}( π_j^x_j( 1 - π_j )^1-x_j )	=	exp( logΠ_j{1→N}( π_j^x_j( 1 - π_j )^1-x_j ) )
	=	exp( Σ_j{1→N}( logπ_j^x_j( 1 - π_j )^1-x_j ) )
	=	exp( Σ_j{1→N}( x_jlogπ_j + ( 1 - x_j )log( 1 - π_j ) ) )
	=	exp( Σ_j{1→N}( x_j[ logπ_j - log( 1 - π_j ) ] + log( 1 - π_j ) ) )
	=	exp( Σ_j{1→N}( x_jlog( π_j / ( 1 - π_j ) ) + log( 1 - π_j ) ) )

Π_i{1→K}( _{n_i}C_{y_i}π_i^y_i( 1 - π_i )^n_i-y_i )	=	exp( logΠ_i{1→K}( _{n_i}C_{y_i}π_i^y_i( 1 - π_i )^n_i-y_i ) )
	=	exp( Σ_i{1→K}( log_{n_i}C_{y_i}π_i^y_i( 1 - π_i )^n_i-y_i ) )
	=	exp( Σ_i{1→K}( log_{n_i}C_{y_i} + y_ilogπ_i + ( n_i - y_i )log( 1 - π_i ) ) )
	=	exp( Σ_i{1→K}( y_ilog( π_i / ( 1 - π_i ) ) + n_ilog( 1 - π_i ) + log_{n_i}C_{y_i} ) )

P(x)	=	0	[ x < a ]
	=	( x - a ) / ( b - a )	[ a ≤ x ≤ b ]
	=	1	[ x > b ]

P(x)	=	( 1 / ( 2πσ² )^1/2 ) ∫{-∞→(x-μ)/σ} exp( -s² / 2 )・σ ds
	=	( 1 / ( 2π )^1/2 ) ∫{-∞→(x-μ)/σ} exp( -s² / 2 ) ds
	=	Φ( ( x - μ ) / σ )

P(x) = ∫{-∞→x^Tα} f(t) dt	=	∫{0→exp( x^Tα )} [ s / ( 1 + s )² ]・[ 1 / s ] ds
	=	∫{0→exp( x^Tα )} ( 1 + s )^-2 ds
	=	[ -( 1 + s )^-1 ]{0→exp( x^Tα )}
	=	1 - 1 / [ 1 + exp( x^Tα ) ]
	=	exp( x^Tα ) / [ 1 + exp( x^Tα ) ]

exp( y・log( π / ( 1 - π ) ) + N・log( 1 - π ) )	=	[ π / ( 1 - π ) ]^y( 1 - π )^N
	=	π^y( 1 - π )^N-y

P(x) = ∫{-∞→x^Tα} f(t) dt	=	∫{0→exp(x^Tα)} [ se^-s ]・[ 1 / s ] ds
	=	∫{0→exp(x^Tα)} e^-s ds
	=	[ -e^-s ]{0→exp(x^Tα)}
	=	1 - exp( -exp( x^Tα ) )

u_j = ∂l / ∂α_j	=	Σ_i{1→K}( x_ijy_i / g'(π_i)π_i( 1 - π_i ) - x_ijn_i / g'(π_i)( 1 - π_i ) )
	=	Σ_i{1→K}( ( y_i - n_iπ_i )x_ij / g'(π_i)π_i( 1 - π_i ) )

l( a\|y )	=	Σ_i{1→K}( y_ilog( π^_i / ( 1 - π^_i ) ) + n_ilog( 1 - π^_i ) + log_{n_i}C_{y_i} )
	=	Σ_i{1→K}( y_ilog( y^_i / ( n_i - y^_i ) ) + n_ilog( ( n_i - y^_i ) / n_i ) + log_{n_i}C_{y_i} )

D = 2[ l( p\|y ) - l( a\|y ) ]	=	2[ Σ_i{1→K}( y_ilog( y_i / ( n_i - y_i ) ) + n_ilog( ( n_i - y_i ) / n_i ) + log_{n_i}C_{y_i} )
		- Σ_i{1→K}( y_ilog( y^_i / ( n_i - y^_i ) ) + n_ilog( ( n_i - y^_i ) / n_i ) + log_{n_i}C_{y_i} ) ]
	=	2Σ_i{1→K}( y_ilog( y_i( n_i - y^_i ) / y^_i( n_i - y_i ) ) + n_ilog( ( n_i - y_i ) / ( n_i - y^_i ) ) )
	=	2Σ_i{1→K}( y_ilog( y_i / y^_i ) - y_ilog( ( n_i - y_i ) / ( n_i - y^_i ) ) + n_ilog( ( n_i - y_i ) / ( n_i - y^_i ) ) )
	=	2Σ_i{1→K}( y_ilog( y_i / y^_i ) + ( n_i - y_i )log( ( n_i - y_i ) / ( n_i - y^_i ) ) )

Χ²	=	Σ_i{1→K}( ( y_i - y^_i )² / y^_i ) + [ ( n_i - y_i ) - ( n_i - y^_i ) ]² / ( n_i - y^_i ) )
	=	Σ_i{1→K}( [ ( n_i - y^_i )( y_i - y^_i )² + y^_i( y^_i - y_i )² ] / y^_i( n_i - y^_i ) )
	=	Σ_i{1→K}( n_i( y_i - y^_i )² / n_iy^_i( 1 - π_i ) )
	=	Σ_i{1→K}( ( y_i - n_iπ_i )² / n_iπ_i( 1 - π_i ) )

2Σ( o・log( o / e ) )	≈	2Σ( ( o - e ) + ( o - e )² / 2e )
	=	Σ( [ e + ( o - e ) ]( o - e ) / e )
	=	Σ( o( o - e ) / e )
	=	Σ_i{1→K}( y_i( y_i - y^_i ) / y^_i +
		( n_i - y_i )[ ( n_i - y_i ) - ( n_i - y^_i ) ] / ( n_i - y^_i ) )
	=	Σ_i{1→K}( y_i( n_i - y^_i )( y_i - y^_i ) / y^_i( n_i - y^_i ) +
		y^_i( n_i - y_i )( y^_i - y_i ) / y^_i( n_i - y^_i ) )
	=	Σ_i{1→K}( [ y_i( n_i - y^_i ) - y^_i( n_i - y_i ) ]( y_i - y^_i ) / y^_i( n_i - y^_i ) )
	=	Σ_i{1→K}( n_i( y_i - y^_i )² / y^_i( n_i - y^_i ) )
	=	Σ_i{1→K}( ( y_i - n_iπ_i )² / n_iπ_i( 1 - π_i ) ) = χ²

∂l / ∂π	=	Σ_i{1→K}( y_i / π( 1 - π ) - n_i / ( 1 - π ) )
	=	Σ_i{1→K}( ( y_i - n_iπ ) / π( 1 - π ) )

C = 2[ l( a\|y ) ) - l( p^\|y ) ]	=	2[ Σ_i{1→K}( y_ilog( y^_i / ( n_i - y^_i ) ) + n_ilog( ( n_i - y^_i ) / n_i ) + log_{n_i}C_{y_i} )
		- Σ_i{1→K}( y_ilog( p^ / ( 1 - p^ ) ) + n_ilog( 1 - p^ ) + log_{n_i}C_{y_i} ) ]
	=	2[ Σ_i{1→K}( y_ilog( y^_i( 1 - p^ ) / p^( n_i - y^_i ) ) + n_ilog( ( n_i - y^_i ) / n_i( 1 - p^ ) ) ) ]
	=	2[ Σ_i{1→K}( -y_ilog( p^( n_i - y^_i ) / y^_i( 1 - p^ ) ) + n_ilog( ( n_i - y^_i ) / n_i( 1 - p^ ) ) ) ]
	=	2[ Σ_i{1→K}( -y_ilog( ( n_ip^ / y^_i )・[ ( n_i - y^_i ) / n_i( 1 - p^ ) ] ) + n_ilog( ( n_i - y^_i ) / n_i( 1 - p^ ) ) ) ]
	=	2[ Σ_i{1→K}( y_ilog( y^_i / n_ip^ ) + ( n_i - y_i )log( ( n_i - y^_i ) / n_i( 1 - p^ ) ) ) ]

C	=	2[ l( a\|y ) - l( p^\|y ) ]
	=	2{ [ l( a\|y ) - l( α\|y ) ] - [ l( p^\|y ) - l( π\|y ) ] - [ l( α\|y ) - l( π\|y ) ] }

x	=	1	(試行が成功)
	=	0	(試行が失敗)

l	=	Σ_i{1→N}( log( P_{n_i,π_i}( y_i ) ) )
	=	Σ_i{1→N}( y_i^Tη( π_i ) + log( h( y_i ) ) )
	=	Σ_i{1→N}( Σ_k{1→C}( y_iklog( π_ik ) ) + log n_i! - Σ_k{1→C}( log y_ik! ) )

∂l / ∂π_ik	=	( ∂ / ∂π_ik )Σ_i'{1→N}( log( P_{n_i',π_i'}( y_i' ) ) )
	=	( ∂ / ∂π_ik )Σ_k'{1→C}( y_ik'log( π_ik' ) )

l	=	Π_i{1→N}( log( P_{n_i,π_i}( y_i ) ) )
	=	Σ_i{1→N}( Σ_k'{1→C}( y_ik'log( π_ik' ) ) + log n_i! - Σ_k'{1→C}( log y_ik'! ) )
	=	Σ_i{1→N}( -y_i1log( 1 + Σ_k'{2→C}( ρ_ik' ) )
		+ Σ_k'{2→C}( y_ik'log( ρ_ik' / [ 1 + Σ_l{2→C}( ρ_il ) ] ) )
		+ log n_i! - Σ_k'{1→C}( log y_ik'! ) )
	=	Σ_i{1→N}( -y_i1log( 1 + Σ_k'{2→C}( ρ_ik' ) )
		+ Σ_k'{2→C}( y_ik'log( ρ_ik' ) - y_ik'log( 1 + Σ_l{2→C}( ρ_il ) ) )
		+ log n_i! - Σ_k'{1→C}( log y_ik'! ) )
	=	Σ_i{1→N}( Σ_k'{2→C}( y_ik'log( ρ_ik' ) ) - n_ilog( 1 + Σ_k'{2→C}( ρ_ik' ) )
		+ log n_i! - Σ_k'{1→C}( log y_ik'! ) )

P_N,π( y )	=	[ N! / Π_k{1→C}( y_k! ) ]exp( log ( Π_k{1→C}( π_k^y_k ) ) )
	=	[ N! / Π_k{1→C}( y_k! ) ]exp( Σ_k{1→C}( y_klog( π_k ) ) )

∂l / ∂π_ik	=	( ∂ / ∂π_ik )[ y_i1log( 1 - Σ_k'{2→C}( π_ik' ) ) + Σ_k'{2→C}( y_ik'log( π_ik' ) ) ]
	=	-y_i1 / ( 1 - Σ_k'{2→C}( π_ik' ) ) + y_ik / π_ik
	=	y_ik / π_ik - y_i1 / π_i1 = 0

Σ_k'{2→C}( ρ_ik' )	=	Σ_k'{2→C}( π_ik' ) / π_i1
	=	( 1 - π_i1 ) / π_i1

u_kj	=	( ∂ / ∂α_kj )Σ_i{1→N}( Σ_k'{2→C}( y_ik'log( ρ_ik' ) ) - n_ilog( 1 + Σ_k'{2→C}( ρ_ik' ) ) )
	=	Σ_i{1→N}( ( y_ik / ρ_ik )( ∂ρ_ik / ∂α_kj )
		- [ n_i / ( 1 + Σ_k'{2→C}( ρ_ik' ) ) ]( ∂ρ_ik / ∂α_kj ) )
	=	Σ_i{1→N}( ( y_ik / ρ_ik )( ∂ρ_ik / ∂α_kj )
		- { n_i / [ 1 + ( 1 - π_i1 ) / π_i1 ] }( ∂ρ_ik / ∂α_kj ) )
	=	Σ_i{1→N}( ( y_ik / ρ_ik - n_iπ_i1 )( ∂ρ_ik / ∂α_kj ) )
	=	Σ_i{1→N}( [ ( y_ik - n_iπ_ik ) / ρ_ik ]( ∂ρ_ik / ∂α_kj ) )

∂ρ_ik / ∂α_kj	=	( dρ_ik / dξ_ik )( ∂ξ_ik / ∂α_kj )
	=	x_ijρ_ik

g'(π_i)	=	( d / dπ_i )log( π_i / ( 1 - π_i ) )
	=	[ ( 1 - π_i ) / π_i ]{ [ ( 1 - π_i ) + π_i ] / ( 1 - π_i )² }
	=	1 / π_i( 1 - π_i )

∂u_kj / ∂α_k'j'	=	( ∂ / ∂α_k'j' )Σ_i{1→N}( ( y_ik - n_iπ_ik )x_ij )
	=	Σ_i{1→N}( -n_ix_ij( ∂π_ik / ∂α_k'j' ) )
	=	Σ_i{1→N}( -n_ix_ij・( ∂ / ∂α_k'j' ){ ρ_ik / [ 1 + Σ_l{2→C}( ρ_il ) ] } )
	=	Σ_i{1→N}( -n_ix_ijρ_ik{ -1 / [ 1 + Σ_l{2→C}( ρ_il ) ]² }( ∂ρ_ik' / ∂α_k'j' ) )
	=	Σ_i{1→N}( n_ix_ijρ_ik[ 1 + ( 1 - π_i1 ) / π_i1 ]²x_ij'ρ_ik' )
	=	Σ_i{1→N}( n_ix_ijx_ij'π_i1²ρ_ikρ_ik' )
	=	Σ_i{1→N}( n_ix_ijx_ij'π_ikπ_ik' )

	( ∂ / ∂α_kj' ){ ρ_ik / [ 1 + Σ_l{2→C}( ρ_il ) ] }
=	{ [ ( 1 + Σ_l{2→C}( ρ_il ) ) - ρ_ik ] / [ 1 + Σ_l{2→C}( ρ_il ) ]² }( ∂ρ_ik / ∂α_kj' )
=	[ ( 1 / π_i1 - π_ik / π_i1 ) / ( 1 / π_i1 )² ]x_ij'ρ_ik
=	π_i1( 1 - π_ik )x_ij'ρ_ik

∂u_kj / ∂α_kj'	=	Σ_i{1→N}( -n_ix_ij[ π_i1( 1 - π_ik )x_ij'ρ_ik ] )
	=	Σ_i{1→N}( -n_ix_ijx_ij'π_i1( 1 - π_ik )ρ_ik )
	=	Σ_i{1→N}( -n_ix_ijx_ij'π_ik( 1 - π_ik ) )

w_ikk'	=	n_iπ_ikπ_ik'	[k ≠ k']
w_ikk	=	-n_iπ_ik( 1 - π_ik )	[k = k']

h_kj^Tα	=	Σ_k'{2→C}( Σ_j'{1→p}( ( ∂u_kj / ∂α_k'j' )α_k'j' ) )
	=	Σ_k'{2→C}( Σ_j'{1→p}( Σ_i{1→N}( w_ikk'x_ijx_ij'α_k'j' ) ) )
	=	Σ_k'{2→C}( Σ_i{1→N}( w_ikk'x_ijΣ_j'{1→p}( x_ij'α_k'j' ) ) )
	=	Σ_k'{2→C}( Σ_i{1→N}( w_ikk'x_ijξ_ik' ) )

x_i1	=	1	(男性)
	=	0	(女性)
x_i2	=	1	(年齢 24-40)
	=	0	(それ以外)
x_i3	=	1	(年齢 >40)
	=	0	(それ以外)

l( π \| y )	=	Σ_i{1→N}( Σ_k{1→C}( y_iklog( π_ik ) ) + log n_i! - Σ_k{1→C}( log y_ik! ) )
	=	Σ_i{1→N}( Σ_k{1→C}( y_iklog( π_ik ) ) ) + K

l( π \| y )	=	Σ_i{1→N}( Σ_k{1→C}( y_iklog( π_k ) ) ) + K
	=	Σ_k{1→C}( Σ_i{1→N}( y_ik )log( π_k ) ) + K
	≡	Σ_k{1→C}( Y_klog( π_k ) ) + K