/*
  OrderedLogisticModel : 順序ロジスティックモデル
*/
class OrderedLogisticModel
{
  bool isValid_; // モデルは有効か

protected:

  const vector< vector<double> >& x_; // 独立変数
  const vector< vector<double> >& y_; // 従属変数(各カテゴリの発生回数)

  LinearEquationSystem<double> s_; // スコア法に利用する連立方程式

  vector<double> a_;               // 回帰係数
  vector< vector<double> > rho_;   // ρ値

  // x は p_ 個のパラメータのベクトルからなる n_ 個のベクトル
  // y は c_ 個のパラメータのベクトルからなる n_ 個のベクトル
  unsigned int n_; // 独立変数ベクトルの数
  unsigned int p_; // 独立変数ベクトルの要素数
  unsigned int c_; // カテゴリ数

  // π の初期化 ( y から計算する )
  void initPi( vector< vector<double> >& pi ) const;

  // ρ の初期化 ( π から計算する )
  virtual void initRho() = 0;
  // ρ から π を求める
  virtual void rho2Pi( vector< vector<double> >& pi ) const = 0;

  // 回帰係数の数を返す
  unsigned int coefSize() const
  { return( a_.size() ); }

  // 回帰係数の数をセットする
  void setCoefSize( unsigned int coefSz )
  {
    if ( ! isValid() ) return;
    a_.resize( coefSz );
    s_.resize( coefSz );
  }

public:

  // コンストラクタ
  OrderedLogisticModel( const vector< vector<double> >& x, const vector< vector<double> >& y );

  virtual void coef2Rho() = 0;       // 求めた係数から ρ を計算する
  virtual void calcCoefMatrix() = 0; // 係数行列を計算する
  virtual void calcRSide() = 0;      // 連立方程式の右辺を計算する
  int solve( double threshold );     // 連立方程式を解いて、収束を調べる

  bool isValid() const { return( isValid_ ); } // モデルは有効か

  void printSize() const // 独立変数、その要素数(元の数)、カテゴリ数を出力
  {
    if ( ! isValid() ) return;
    cout << "N = " << n_ << " ; p = " << p_ << " ; c = " << c_ << endl << endl;
  }
  void printX() const // 独立変数を出力
  {
    if ( ! isValid() ) return;
    PrintMatrix( "x = ", x_ ); cout << endl;
  }
  void printY() const // 従属変数(各カテゴリの発生回数)を出力
  {
    if ( ! isValid() ) return;
    PrintMatrix( "y = ", y_ ); cout << endl;
  }
  void printPi() const; // π を出力
  void printRho() const // ρ を出力
  {
    if ( ! isValid() ) return;
    PrintMatrix( "rho = ", rho_ ); cout << endl;
  }
  void printLES() const // 連立方程式を出力
  {
    if ( ! isValid() ) return;
    s_.print(); cout << endl;
  }

  virtual void printVar() = 0; // 回帰係数に対する分散の出力
  virtual void printEquation() const = 0; // 回帰式の出力

  virtual string ident() const = 0; // モデル名
};

/*
  OrderedLogisticModel コンストラクタ

  const vector< vector<double> >& x : 独立変数
  const vector< vector<double> >& y : 従属変数(各カテゴリの発生回数)
  unsigned int coefSz : 回帰係数の数

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

  x は p 個のパラメータのベクトルからなる n 個のベクトル
  y は c 個のパラメータのベクトルからなる n 個のベクトル
*/
OrderedLogisticModel::OrderedLogisticModel( const vector< vector<double> >& x, const vector< vector<double> >& y )
  : isValid_( false ), x_( x ), y_( y ), s_( 0 )
{
  // NULL のチェック
  if ( ! NullCheck( x_, "Independent Variable x" ) ) return;
  if ( ! NullCheck( y_, "Occurrence Count y" ) ) return;

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

  // 独立変数ベクトルの要素数
  p_ = x[0].size();
  if ( ! SizeCheck_Loop( x_, "Independent Variable x", p_, "Independent Variable x[0]" ) )
    return;

  // カテゴリ数
  c_ = y[0].size();
  if ( ! SizeCheck_Loop( y_, "Occurrence Count y", c_, "Occurrence Count y[0]" ) )
    return;
  if ( c_ == 0 ) {
    cerr << "The size of categories is zero." << endl;
    return;
  }

  isValid_ = true;
}

/*
  OrderedLogisticModel::initPi : π の初期化 ( y から計算する )

  vector< vector<double> >& pi : 計算結果を格納するベクトル
*/
void OrderedLogisticModel::initPi( vector< vector<double> >& pi ) const
{
  if ( ! isValid() ) return;

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

  pi.assign( 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];
}

/*
  ProportionalOddsModel::solve : 連立方程式を解いて、収束を調べる

  double threshold : 収束条件(全係数が threshold 以下なら処理終了)

  戻り値 : >0 収束した ; =0 収束しない ; <0 計算に失敗
*/
int OrderedLogisticModel::solve( double threshold )
{
  if ( ! isValid() ) return( -1 );

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

  // 各係数が収束しているかを確認する
  bool isMatched = true;
  for ( unsigned int i = 0 ; i < coefSize() ; ++i ) {
    if ( fabs( a_[i] - s_.ans( i ) ) >= threshold )
      isMatched = false;
    a_[i] = s_.ans( i );
  }

  return( ( isMatched ) ? 1 : 0 );
}

/*
  OrderedLogisticModel::printPi : π の出力
*/
void OrderedLogisticModel::printPi() const
{
  if ( ! isValid() ) return;

  vector< vector<double> > pi; // π の計算結果

  rho2Pi( pi );
  PrintMatrix( "pi = ", pi );
  cout << endl;
}

/*------------------------------------------------------------------------------------------------
  CumulativeLogitModel_Base : 累積ロジットモデルの基底クラス
  ------------------------------------------------------------------------------------------------*/
class CumulativeLogitModel_Base : public OrderedLogisticModel
{
protected:

  // ρ の初期化 ( π から計算する )
  virtual void initRho();
  // ρ から π を求める
  virtual void rho2Pi( vector< vector<double> >& pi ) const;

public:

  // コンストラクタ
  CumulativeLogitModel_Base( const vector< vector<double> >& x, const vector< vector<double> >& y )
    : OrderedLogisticModel( x, y )
    { if ( isValid() ) initRho(); }
};

/*
  CumulativeLogitModel_Base::initRho : ρ の初期化 ( y から計算する )
*/
void CumulativeLogitModel_Base::initRho()
{
  if ( ! isValid() ) return;

  // π の計算 ( y から求める )
  vector< vector<double> > pi;
  initPi( pi );

  rho_.assign( n_, vector<double>( c_ - 1 ) );
  for ( unsigned int i = 0 ; i < n_ ; ++i ) {
    double sumPi = 0; // i までの累積確率
    for ( unsigned int k = 0 ; k < c_ - 1 ; ++k ) {
      sumPi += pi[i][k];
      if ( sumPi == 0 ) sumPi = 1E-9; // ゼロ除算防止
      if ( sumPi > 1 ) sumPi = 1;
      rho_[i][k] = ( 1.0 - sumPi ) / sumPi;
    }
  }
}

/*
  CumulativeLogitModel_Base::rho2Pi : ρ から π を求める

  vector< vector<double> >& pi : 計算結果を格納するベクトル
*/
void CumulativeLogitModel_Base::rho2Pi( vector< vector<double> >& pi ) const
{
  if ( ! isValid() ) return;

  pi.assign( n_, vector<double>( c_ ) ); // 確率の当てはめ値 ( n x c )
  for ( unsigned int i = 0 ; i < n_ ; ++i ) {
    pi[i][0] = 1.0 / ( 1.0 + rho_[i][0] );
    for ( unsigned int k = 1 ; k < c_ - 1 ; ++k )
      pi[i][k] = 1.0 / ( 1.0 + rho_[i][k] ) - 1.0 / ( 1.0 + rho_[i][k - 1] );
    if ( c_ > 1 )
      pi[i][c_ - 1] = 1.0 - 1.0 / ( 1.0 + rho_[i][c_ - 2] );
  }
}

/*------------------------------------------------------------------------------------------------
  ProportionalOddsModel : 比例オッズモデル
  ------------------------------------------------------------------------------------------------*/
class ProportionalOddsModel : public CumulativeLogitModel_Base
{
  void calcH11(); // H11 の計算
  void calcH22(); // H22 の計算
  void calcH12(); // H12(H21) の計算

public:

  // コンストラクタ
  ProportionalOddsModel( const vector< vector<double> >& x, const vector< vector<double> >& y )
    : CumulativeLogitModel_Base( x, y )
    { setCoefSize( c_ + p_ - 1 ); }

  virtual void coef2Rho();       // 求めた係数から ρ を計算する
  virtual void calcCoefMatrix(); // 係数行列を計算する
  virtual void calcRSide();      // 連立方程式の右辺を計算する

  virtual void printVar();            // 回帰係数に対する分散の出力
  virtual void printEquation() const; // 回帰式の出力

  // モデル名
  virtual string ident() const
  { return( "Proportional Odds Model" ); }
};

/*
  ProportionalOddsModel::coef2Rho : 求めた係数から ρ を計算する
*/
void ProportionalOddsModel::coef2Rho()
{
  if ( ! isValid() ) return;

  for ( unsigned int i = 0 ; i < n_ ; ++i ) {
    for ( unsigned int k = 0 ; k < c_ - 1 ; ++k ) {
      double xi = a_[k];
      for ( unsigned int j = 0 ; j < p_ ; ++j )
        xi += a_[c_ - 1 + j] * x_[i][j];
      rho_[i][k] = exp( xi );
    }
  }
}

/*
  ProportionalOddsModel::calcH11 : H11 の計算
*/
void ProportionalOddsModel::calcH11()
{
  for ( unsigned int k = 0 ; k < c_ - 1 ; ++k ) {
    // 対角成分を計算
    s_[k][k] = 0;
    for ( unsigned int i = 0 ; i < n_ ; ++i ) {
      s_[k][k] -= ( y_[i][k] + y_[i][k + 1] ) * rho_[i][k] / pow( 1.0 + rho_[i][k], 2 );
      if ( k + 1 < c_ - 1 )
        s_[k][k] -= y_[i][k + 1] * rho_[i][k] * rho_[i][k + 1] / pow( rho_[i][k] - rho_[i][k + 1], 2 );
      if ( k > 0 )
        s_[k][k] -= y_[i][k] * rho_[i][k] * rho_[i][k - 1] / pow( rho_[i][k] - rho_[i][k - 1], 2 );
    }
    if ( k + 1 == c_ - 1 ) break;
    // 対角成分の右隣を計算
    s_[k][k + 1] = 0;
    for ( unsigned int i = 0 ; i < n_ ; ++i )
      s_[k][k + 1] += y_[i][k + 1] * rho_[i][k] * rho_[i][k + 1] / pow( rho_[i][k] - rho_[i][k + 1], 2 );
    s_[k + 1][k] = s_[k][k + 1]; // 対称な成分にコピー
    // 残りはゼロで初期化
    for ( unsigned int k2 = k + 2 ; k2 < c_ - 1 ; ++k2 )
      s_[k][k2] = s_[k2][k] = 0;
  }
}

/*
  ProportionalOddsModel::calcH22 : H22 の計算
*/
void ProportionalOddsModel::calcH22()
{
  // k ごとの和の部分を先に計算
  vector<double> sum( n_, 0 );
  for ( unsigned int i = 0 ; i < n_ ; ++i )
    for ( unsigned int k = 0 ; k < c_ - 1 ; ++k )
      sum[i] += ( y_[i][k] + y_[i][k + 1] ) * rho_[i][k] / pow( 1.0 + rho_[i][k], 2 );

  for ( unsigned int sr = c_ - 1 ; sr < coefSize() ; ++sr ) { // sr は係数行列の行番号
    unsigned int jr = sr - ( c_ - 1 ); // x_ に対する行番号
    for ( unsigned int sc = sr ; sc < coefSize() ; ++sc ) { // sc は係数行列の列番号(対角成分からスタート)
      unsigned int jc = sc - ( c_ - 1 ); // x_ に対する列番号
      s_[sr][sc] = 0;
      for ( unsigned int i = 0 ; i < n_ ; ++i )
        s_[sr][sc] -= x_[i][jr] * x_[i][jc] * sum[i];
      // 対角成分(左下)はコピー
      if ( sc != sr )
        s_[sc][sr] = s_[sr][sc];
    }
  }
}

/*
  ProportionalOddsModel::calcH12 : H12(H21) の計算
*/
void ProportionalOddsModel::calcH12()
{
  // 計算は右上側(H12)に対して行う
  for ( unsigned int k = 0 ; k < c_ - 1 ; ++k ) {
    for ( unsigned int sc = c_ - 1 ; sc < coefSize() ; ++sc ) {
      unsigned int j = sc - ( c_ - 1 ); // x_ に対する列番号
      s_[k][sc] = 0;
      for ( unsigned int i = 0 ; i < n_ ; ++i )
        s_[k][sc] -= x_[i][j] * ( y_[i][k] + y_[i][k + 1] ) * rho_[i][k] / pow( 1.0 + rho_[i][k], 2 );
      s_[sc][k] = s_[k][sc]; // 対角成分(左下)はコピー
    }
  }
}

/*
  ProportionalOddsModel::calcCoefMatrix : 係数行列を計算する
*/
void ProportionalOddsModel::calcCoefMatrix()
{
  if ( ! isValid() ) return;

  calcH11(); // H11 の計算
  calcH22(); // H22 の計算
  calcH12(); // H12(H21) の計算
}

/*
  ProportionalOddsModel::calcRSide : 連立方程式の右辺を計算する
*/
void ProportionalOddsModel::calcRSide()
{
  if ( ! isValid() ) return;

  vector<double> u_sum( n_, 0 ); // 対数尤度の偏微分 u の計算に必要な項
  vector<double> h_sum( n_, 0 ); // 係数行列と係数の積 Ha の計算に必要な項
  for ( unsigned int i = 0 ; i < n_ ; ++i ) {
    for ( unsigned int k = 0 ; k < c_ - 1 ; ++k ) {
      u_sum[i] += ( y_[i][k] + y_[i][k + 1] ) / ( 1.0 + rho_[i][k] ) - y_[i][k];
      h_sum[i] += ( y_[i][k] + y_[i][k + 1] ) * rho_[i][k] * log( rho_[i][k] )
        / pow( 1.0 + rho_[i][k], 2 );
    }
  }

  vector<double> u( coefSize(), 0 ); // u の値
  vector<double> h( coefSize(), 0 ); // Ha の値
  for ( unsigned int k = 0 ; k < c_ - 1 ; ++k ) {
    for ( unsigned int i = 0 ; i < n_ ; ++i ) {
      u[k] += ( y_[i][k] + y_[i][k + 1] ) / ( 1.0 + rho_[i][k] );
      h[k] -= ( y_[i][k] + y_[i][k + 1] ) * rho_[i][k] * log( rho_[i][k] ) / pow( 1.0 + rho_[i][k], 2 );
      if ( k + 1 < c_ - 1 ) {
        u[k] += y_[i][k + 1] * rho_[i][k + 1] / ( rho_[i][k] - rho_[i][k + 1] );
        h[k] += y_[i][k + 1] * rho_[i][k] * rho_[i][k + 1]
          * ( log( rho_[i][k + 1] ) - log( rho_[i][k] ) ) / pow( rho_[i][k] - rho_[i][k + 1], 2 );
      }
      if ( k > 0 ) {
        u[k] += y_[i][k] * rho_[i][k - 1] / ( rho_[i][k] - rho_[i][k - 1] );
        h[k] += y_[i][k] * rho_[i][k] * rho_[i][k - 1]
          * ( log( rho_[i][k - 1] ) - log( rho_[i][k] ) ) / pow( rho_[i][k] - rho_[i][k - 1], 2 );
      } else {
        u[k] -= y_[i][k];
      }
    }
  }
  for ( unsigned int j = 0 ; j < p_ ; ++j ) {
    for ( unsigned int i = 0 ; i < n_ ; ++i ) {
      u[j + c_ - 1] += x_[i][j] * u_sum[i];
      h[j + c_ - 1] -= x_[i][j] * h_sum[i];
    }
  }
  for ( unsigned int k = 0 ; k < coefSize() ; ++k )
    s_.ans( k ) = h[k] - u[k]; // 右辺 = Ha - u
}

/*
  ProportionalOddsModel::printVar : 回帰係数に対する分散の出力
*/
void ProportionalOddsModel::printVar()
{
  if ( ! isValid() ) return;

  // 係数行列の再計算(=フィッシャー情報行列)
  calcCoefMatrix();

  LinearEquationSystem<double> inv( 0 ); // 連立方程式計算用インスタンスの逆行列
  Inverse( s_, inv );

  cout << "variance of b = ( " << -inv[0][0];
  for ( unsigned int k = 1 ; k < c_ - 1 ; ++k )
    cout << ", " << -inv[k][k];
  cout << " )" << endl;

  cout << "variance of a = ( " << -inv[c_ - 1][c_ - 1];
  for ( unsigned int j = c_ ; j < coefSize() ; ++j )
    cout << ", " << -inv[j][j];
  cout << " )" << endl << endl;
}

/*
  ProportionalOddsModel::printEquation : 回帰式の出力
*/
void ProportionalOddsModel::printEquation() const
{
  if ( ! isValid() ) return;

  vector<double> a2( a_.begin() + c_ - 1, a_.end() );
  for ( unsigned int k = 0 ; k < c_ - 1 ; ++k ) {
    std::ostringstream oss;
    oss << "Regression equation : y[" << k << "] = " << a_[k] << " + ";
    PrintEquation( oss.str(), a2 );
  }
  cout << endl;
}

/*------------------------------------------------------------------------------------------------
  CumulativeLogitModel : 累積ロジットモデル
  ------------------------------------------------------------------------------------------------*/
class CumulativeLogitModel : public CumulativeLogitModel_Base
{
  void calcCoefDiagBlock( unsigned int k );
  void calcCoefNonDiagBlock( unsigned int k );

public:

  // コンストラクタ
  CumulativeLogitModel( const vector< vector<double> >& x, const vector< vector<double> >& y )
    : CumulativeLogitModel_Base( x, y )
    { setCoefSize( ( c_ - 1 ) * p_ ); }

  virtual void coef2Rho();       // 求めた係数から ρ を計算する
  virtual void calcCoefMatrix(); // 係数行列を計算する
  virtual void calcRSide();      // 連立方程式の右辺を計算する

  virtual void printVar();            // 回帰係数に対する分散の出力
  virtual void printEquation() const; // 回帰式の出力

  // モデル名
  virtual string ident() const
  { return( "Cumulative Logit Model" ); }
};

/*
  CumulativeLogitModel::coef2Rho : 求めた係数から ρ を計算する
*/
void CumulativeLogitModel::coef2Rho()
{
  if ( ! isValid() ) return;

  for ( unsigned int i = 0 ; i < n_ ; ++i ) {
    for ( unsigned int k = 0 ; k < c_ - 1 ; ++k ) {
      double xi = a_[k * p_] * x_[i][0];
      for ( unsigned int j = 1 ; j < p_ ; ++j )
        xi += a_[k * p_ + j] * x_[i][j];
      rho_[i][k] = exp( xi );
    }
  }
}

/*
  CumulativeLogitModel::calcCoefDiagBlock : 係数行列の対角部分行列を計算する
*/
void CumulativeLogitModel::calcCoefDiagBlock( unsigned int k )
{
  // k に依存する部分を先に計算
  vector<double> sum( n_, 0 );
  for ( unsigned int i = 0 ; i < n_ ; ++i ) {
    sum[i] += ( y_[i][k] + y_[i][k + 1] ) * rho_[i][k] / pow( 1.0 + rho_[i][k], 2 );
    if ( k + 1 < c_ - 1 )
      sum[i] += y_[i][k + 1] * rho_[i][k] * rho_[i][k + 1] / pow( rho_[i][k] - rho_[i][k + 1], 2 );
    if ( k > 0 )
      sum[i] += y_[i][k] * rho_[i][k] * rho_[i][k - 1] / pow( rho_[i][k] - rho_[i][k - 1], 2 );
  }

  for ( unsigned int jr = 0 ; jr < p_ ; ++jr ) {
    for ( unsigned int jc = 0 ; jc < p_ ; ++jc ) {
      s_[k * p_ + jr][k * p_ + jc] = 0;
      for ( unsigned int i = 0 ; i < n_ ; ++i )
        s_[k * p_ + jr][k * p_ + jc] -= x_[i][jr] * x_[i][jc] * sum[i];
    }
  }
}

/*
  CumulativeLogitModel::calcCoefNonDiagBlock : 係数行列の非対角部分行列(対角成分の右側)を計算する
*/
void CumulativeLogitModel::calcCoefNonDiagBlock( unsigned int k )
{
  // k に依存する部分を先に計算
  vector<double> sum( n_, 0 );
  for ( unsigned int i = 0 ; i < n_ ; ++i )
    sum[i] += y_[i][k + 1] * rho_[i][k] * rho_[i][k + 1] / pow( rho_[i][k] - rho_[i][k + 1], 2 );

  for ( unsigned int jr = 0 ; jr < p_ ; ++jr ) {
    for ( unsigned int jc = 0 ; jc < p_ ; ++jc ) {
      s_[k * p_ + jr][( k + 1 ) * p_ + jc] = 0;
      for ( unsigned int i = 0 ; i < n_ ; ++i )
        s_[k * p_ + jr][( k + 1 ) * p_ + jc] += x_[i][jr] * x_[i][jc] * sum[i];
      // 対称成分にコピー
      s_[( k + 1 ) * p_ + jc][k * p_ + jr] = s_[k * p_ + jr][( k + 1 ) * p_ + jc];
    }
  }
}

/*
  CumulativeLogitModel::calcCoefMatrix : 係数行列を計算する
*/
void CumulativeLogitModel::calcCoefMatrix()
{
  if ( ! isValid() ) return;

  for ( unsigned int k = 0 ; k < c_ - 1 ; ++k ) {
    calcCoefDiagBlock( k );
    if ( k + 1 < c_ - 1 ) calcCoefNonDiagBlock( k );
  }
}

/*
  CumulativeLogitModel::calcRSide : 連立方程式の右辺を計算する
*/
void CumulativeLogitModel::calcRSide()
{
  if ( ! isValid() ) return;

  vector< vector<double> > u_sum( n_, vector<double>( c_ - 1 ) ); // 対数尤度の偏微分 u の計算に必要な項
  vector< vector<double> > h_sum( n_, vector<double>( c_ - 1 ) ); // 係数行列と係数の積 Ha の計算に必要な項
  // k に依存する部分を先に計算しておく
  for ( unsigned int k = 0 ; k < c_ - 1 ; ++k ) {
    for ( unsigned int i = 0 ; i < n_ ; ++i ) {
      h_sum[i][k] = -( y_[i][k] + y_[i][k + 1] ) * rho_[i][k] * log( rho_[i][k] )
        / pow( 1.0 + rho_[i][k], 2 );
      u_sum[i][k] = ( y_[i][k] + y_[i][k + 1] ) / ( 1.0 + rho_[i][k] );
      if ( k > 0 ) {
        h_sum[i][k] += y_[i][k] * rho_[i][k] * rho_[i][k - 1]
          * ( log( rho_[i][k - 1] ) - log( rho_[i][k] ) ) / pow( rho_[i][k] - rho_[i][k - 1], 2 );
        u_sum[i][k] += y_[i][k] * rho_[i][k - 1] / ( rho_[i][k] - rho_[i][k - 1] );
      } else {
        u_sum[i][k] -= y_[i][k];
      }
      if ( k + 1 < c_ - 1 ) {
        h_sum[i][k] += y_[i][k + 1] * rho_[i][k] * rho_[i][k + 1]
          * ( log( rho_[i][k + 1] ) - log( rho_[i][k] ) ) / pow( rho_[i][k] - rho_[i][k + 1], 2 );
        u_sum[i][k] += y_[i][k + 1] * rho_[i][k + 1] / ( rho_[i][k] - rho_[i][k + 1] );
      }
    }
  }

  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 i = 0 ; i < n_ ; ++i )
        s_.ans( k * p_ + j ) += x_[i][j] * ( h_sum[i][k] - u_sum[i][k] );
    }
  }
}

/*
  CumulativeLogitModel::printVar : 回帰係数に対する分散の出力
*/
void CumulativeLogitModel::printVar()
{
  if ( ! isValid() ) return;

  // 係数行列の再計算(=フィッシャー情報行列)
  calcCoefMatrix();

  LinearEquationSystem<double> inv( 0 ); // 連立方程式計算用インスタンスの逆行列
  Inverse( s_, inv );

  for ( unsigned int k = 0 ; k < c_ - 1 ; ++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;
  }
  cout << endl;
}

/*
  CumulativeLogitModel::printEquation : 回帰式の出力
*/
void CumulativeLogitModel::printEquation() const
{
  if ( ! isValid() ) return;

  for ( unsigned int k = 0 ; k < c_ - 1 ; ++k ) {
    vector<double> ak( a_.begin() + k * p_, a_.begin() + ( k + 1 ) * p_ );
    std::ostringstream oss;
    oss << " y[" << k << "] = ";
    PrintEquation( oss.str(), ak );
  }
  cout << endl;
}

/*
  OrderedLogisticRegression : 順序ロジスティック回帰

  OrderedLogisticModel& model : 順序ロジスティックモデル
  bool verbose : 冗長モード(ON/OFF)
  unsigned int maxCount : 反復処理の最大回数
  double threshold : 収束条件(全係数が threshold 以下なら処理終了)

  戻り値 : 係数が得られた ... true ; データ異常・反復処理回数が最大値を超えた ... false
*/
bool OrderedLogisticRegression( OrderedLogisticModel& model,
                                bool verbose, unsigned int maxCount, double threshold )
{
  cout << "*** " << model.ident() << " ***" << endl << endl;

  if ( ! model.isValid() ) {
    cerr << "The model is invalid." << endl;
    return( false );
  }

  model.printSize();

  if ( verbose ) {
    model.printX();
    model.printY();
  }

  unsigned int cnt; // 計算回数
  for ( cnt = 0 ; cnt < maxCount ; ++cnt ) {

    if ( verbose ) {
      cout << "----- cnt = " << cnt << " -----" << endl << endl;
      model.printPi();
      model.printRho();
    }

    // 係数行列の計算
    model.calcCoefMatrix();
    // 右辺の計算
    model.calcRSide();

    if ( verbose ) {
      cout << "Equation System :" << endl;
      model.printLES();
    }

    // 連立方程式の計算
    int res = model.solve( threshold );
    if ( res < 0 ) return( false );

    if ( verbose )
      model.printEquation();

    if ( res > 0 ) break;

    // カテゴリごとの確率計算
    model.coef2Rho();
  }

  if ( cnt < maxCount ) {
    cout << "Estimated regression equation" << endl << endl;
    model.printEquation();

    model.printVar();

    cout << "Estimated probability" << endl;
    model.printPi();
  } else {
    cout << "Failed to estimate regression coefficient" << endl << endl;
  }

  return( cnt < maxCount );
}

vector< vector<double> > x;  // 定数項を含まない独立変数
vector< vector<double> > x2; // (必要なら)定数項を含む独立変数
vector< vector<double> > y;  // 従属変数(各カテゴリの発生回数)

/*
  x[i], x2[i] が独立変数ベクトルで、x[i][j], x2[i][j] がその要素になるように初期化
  但し、x2[i] は(必要なら)定数項を含むのに対し、x[i] は定数項を含めないことに注意
  同様に、y[i] が従属変数ベクトルで、y[i][j] がその要素になるように初期化
*/

// 比例オッズモデル
ProportionalOddsModel proModel( x, y );
OrderedLogisticRegression( proModel, true );

// 累積ロジットモデル
CumulativeLogitModel cumModel( x2, y );
OrderedLogisticRegression( cumModel, true );