// logreg.C Version 7.45  Copyleft (c) 杉本典夫
// Time-stamp: <2026-04-29 09:43:29 tac>
//
// ロジスティック(ロジット)回帰分析。原則として前向きデータに適用する。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法  :g++ -O -s -W logreg logreg.C
//                -I/usr/local/include -L/usr/local/lib -ldans -lfun

#include <dans.H>               // dans用ヘッダファイル(fun.H,mat.H含む)

using namespace dans;           // dans名前空間使用(名前空間matvec含む)

const vecstr USAGE =            // プログラム使用法表示用
{
  "logreg  Version 7.45         Copyleft(Vo\\oV) 2026 杉本典夫",
  "書式:logreg [<オプション> …]  [<データファイル> …]",
  "機能:最初の項目を群項目(目的変数Y)、以後の項目を説明変数Xとした",
  "     逐次変数選択法によるロジスティック(ロジット)回帰分析を行います。",
  "<オプション>",
  "  -b=<β0の初期値>[,<β1の初期値>,…]:偏回帰係数βの初期値を指定します。",
  "  -cc[=<信頼係数>]:信頼係数を指定します。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -cin=<取り込み基準>:変数の取り込み基準Cinを指定します。(Cin≧0)",
  "  -comp[={-|+}]:ロジットスコア計算時の欠損を補間します。(-comp:-comp=+)",
  "  -dnsname[=<出力ファイル名>]:DANSデータを出力するファイル名を指定する。",
  "    オプション省略時:DANSデータ出力せず  -dnsname:-dnsname=logreg_dns.d",
  "  -eps[=<収束誤差>]:最尤法の収束誤差を指定します。(0<収束誤差)",
  "  -f=<ファイル>:指定したファイルからオプションを読み込みます。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示します。(-h:-h=man)",
  "  -ix=<説明変数の開始番号>[-<終了番号>]x<説明変数の開始番号>[-<終了番号>][,…]",
  "    交互作用項目を作成する2つの説明変数の番号を指定します。",
  "  -ixc[={-|+}]:交互作用項目の説明変数を中心化するかどうか指定します。",
  "    -:中心化しない(オプション省略時)  +:中心化する  -ixc:ixc=+",
  "  -ngp={lci|log|or|score}[,…]:Ngraph用データを指定します。",
  "    lci:ロジスティック曲線の信頼区間  log:ロジスティック曲線  or:オッズ比",
  "    score:ロジットスコア  オプション省略時:全て指定",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データを出力するファイル名を指定する。",
  "    オプション省略時:Ngraph用データ出力せず  -ngpname:-ngpname=logreg_ngp.dat",
  "  -pin=<取り込み基準>:変数の取り込み基準Pinを指定します。(0＜Pin≦1)",
  "  -r{p|n}[=<コード>]:反応有(rp)または反応無(rn)のコードを指定します。",
  "  -res=[-][{def|all|allx|ltab|selx|stat|v}][,…]:",
  "      解析結果の出力項目を指定します。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    allx:全変数選択結果  ltab:尤度分析表  selx:変数選択結果",
  "    stat:基礎統計量  v:偏回帰係数初期値",
  "  -rmax=<最大回数>:反復計算の最大回数を指定します。(0≦最大回数)",
  "  -x<説明変数番号>=<説明変数の値>:X軸にしない説明変数の値を指定します。",
  "  -xno=<説明変数の番号>:X軸にする説明変数の番号を指定します。",
  "  -xrange=<最小値>-<最小値>:X軸の最小値と最大値を指定します。"
};

const string Title = "=== ロジスティック回帰分析(logistic regression analysis) ==="; // 表題
const string DNSNAMEdef = "logreg_dns.d"; // 省略時DANSファイル名
const string NGPNAMEdef = "logreg_ngp.dat"; // 省略時Ngraphファイル名
constexpr double EPSdef = 1.0e-5; // 収束誤差省略時値
constexpr double CIOdef = 2.0; // χ2乗値取り込み・追い出し基準省略時値
constexpr double PIOdef = 0.1; // p値取り込み・追い出し基準省略時値
constexpr int RMAXdef = 200;   // 反復計算最大回数の省略時値
constexpr int RDdef = 50;      // グラフの描画範囲分割数省略時値

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  lstr blst;                    // 説明変数の初期値リスト
  vecdbl v_bvec;                // 説明変数の初期値ベクトル
  double v_cc;                  // 信頼係数
  double v_cin;                 // χ2乗値の取り込み基準
  bool v_comp;                  // 欠損補間の有無
  string v_dnsname;             // DANSファイル名
  double v_eps;                 // 最尤法の収束誤差
  lstr ixstrlst;                // 交互作用項目の説明変数リスト
  lino v_ixlst;                 // 交互作用項目の説明変数ペアリスト
  bool v_ixc;                   // 交互作用項目の中心化フラッグ
  string v_ngpname;             // Ngraphファイル名
  lstr ngplst;                  // Ngraph出力リスト
  mapsb ngpmap;                 // Ngraph出力の有無マップ
  double v_pin;                 // p値の取り込み基準
  vecstr rvec;                  // 反応有・無のコード
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
  int v_rmax;                   // 最大反復回数
  lstr xlst;                    // 説明変数リスト
  vecdbl xvec;                  // 説明変数ベクトル
  int v_xno;                    // X軸にする説明変数の番号
  pairdbl xrange;               // X軸の最小値と最大値
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  const vecdbl& bvec() const    // 説明変数の初期値ベクトルを返す
  { return v_bvec; }
  double b(size_t idx) const    // 説明変数の初期値を返す
  { return (idx < v_bvec.size()) ? v_bvec[idx] : 0.0; }
  double cc() const             // 信頼係数を返す
  { return v_cc; }
  double cin() const            // χ2乗値の取り込み基準を返す
  { return v_cin; }
  bool comp() const             // 欠損補間の有無を返す
  { return v_comp; }
  const string& dnsname() const // DANSデータファイル名を返す
  { return v_dnsname; }
  double eps() const            // 最尤法の収束誤差を返す
  { return v_eps; }
  const lino& ixlst()    // 交互作用項目の説明変数番号ペアリストを返す
  { return v_ixlst; }
  bool ixc()                   // 交互作用項目の中心化フラッグを返す
  { return v_ixc; }
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  bool ngp(const string& code); // Ngraph出力フラッグを返す
  double pin() const            // p値の取り込み基準を返す
  { return v_pin; }
  const string& rn() const      // 反応無のコードを返す
  { return rvec[0]; }
  const string& rp() const      // 反応有のコードを返す
  { return rvec[1]; }
  bool res(const string& code); // 分析結果出力フラッグを返す
  int rmax() const              // 反復計算の最大回数を返す
  { return v_rmax; }
  double x(size_t xno) const    // 説明変数の値を返す
  { return (0 < xno && xno <= xvec.size()) ? xvec[xno - 1] : BLANKf; }
  int xno() const               // X軸にする説明変数の番号を返す
  { return v_xno; }
  double xmin() const           // X軸の最小値を返す
  { return xrange.first; }
  double xmax() const           // X軸の最大値を返す
  { return xrange.second; }
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check(size_t ni);         // オプションのチェックとセット
  int set_ngp();                // -ngpオプションのセット
  int set_res();                // -resオプションのセット
  int set_ix(int nx);           // -ixオプションのセット
  int set_x(int nx);            // -xオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), blst(), v_bvec(), v_cc(CCdef), v_cin(CIOdef), v_comp(false), v_dnsname(), ixstrlst(), v_ixlst(), v_ixc(false), v_eps(EPSdef), v_ngpname(), ngplst(), ngpmap(), v_pin(1.0), rvec(2), reslst(), resmap(), v_rmax(RMAXdef), xlst(), xvec(), v_xno(0), xrange(BLANKf, BLANKf)
// コンストラクタ
{
  if (c_pname == nullptr) return; // 初期設定せず
  v_pname = getname(c_pname);     // パス名なしのプログラム名取得
  v_path = c_pname;
  strsize_t idx = v_path.find(v_pname);
  if (idx != string::npos) v_path.erase(idx); // プログラム名削除
  rvec[1] = "1";
}

bool option::ngp(const string& code) // Ngraph出力フラッグを返す
{
  mapsb_citr ngpitr = ngpmap.find(code);
  if (ngpitr == ngpmap.end()) return false; // code無し！
  return ngpmap[code];                      // true/falseを返す
}

bool option::res(const string& code) // 分析結果出力フラッグを返す
{
  mapsb_citr resitr = resmap.find(code);
  if (resitr == resmap.end()) return false; // code無し！
  return resmap[code];                      // true/falseを返す
}

void option::set(const string& optstr, char devc)
// オプション指定によってオプションを設定する
{
  string opt[2];                // オプションのopt[0]:項目,opt[1]:内容
  if (getopt(opt, optstr, devc) < 1) return; // オプション指定なし
  int err = 0;
  if (opt[0] == "b")            // 偏回帰係数初期値の指定
  {
    if (strsplit(blst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "cc")      // 信頼係数の指定
  {
    v_cc = opt[1].empty() ? CCdef : strtodf(opt[1]);
    if (v_cc < 0.0 || 1.0 <= v_cc) err = badopt(pname(), optstr, USAGE, 0.0, 1.0);
  }
  else if (opt[0] == "cin")     // χ2乗値の取り込み基準指定
  {
    v_cin = opt[1].empty() ? CIOdef : strtodf(opt[1]);
    if (v_cin < 0.0) err = badopt(pname(), optstr, USAGE, 0.0, FDmax);
  }
  else if (opt[0] == "comp")    // 欠損補完の有無の指定
  {
    if (opt[1].empty() || opt[1] == "+") v_comp = true;
    else if (opt[1] == "-") v_comp = false;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "dnsname") // DANSファイル名指定
  {
    if (opt[1].empty()) v_dnsname = DNSNAMEdef;
    else v_dnsname = opt[1];
    if (v_dnsname.size() < 3) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "eps")     // 最尤法の収束誤差の指定
  {
    v_eps = opt[1].empty() ? EPSdef : strtodf(opt[1]);
    if (v_eps < FDeps) err = badopt(pname(), optstr, USAGE, FDeps, FDmax);
  }
  else if (opt[0] == "f")       // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "ix")      // 交互作用項目の説明変数ペア指定
  {
    if (strsplit(ixstrlst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "ixc")     // 交互作用項目の説明変数の中心化指定
  {
    if (opt[1].empty() || opt[1] == "+" ) v_ixc = true;
    else if (opt[1] == "-") v_ixc = false;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else if (opt[0] == "ngp")     // Ngraph出力項目指定
  {
    if (strsplit(ngplst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "ngpname")     // Ngraphファイル名指定
  {
    if (opt[1].empty()) v_ngpname = NGPNAMEdef;
    else v_ngpname = opt[1];
    if (v_ngpname.size() < 3) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "pin")     // p値の取り込み基準指定
  {
    v_pin = opt[1].empty() ? PIOdef : strtodf(opt[1]);
    if (v_pin <= 0.0 || 1.0 < v_pin) err = badopt(pname(), optstr, USAGE, 0.0, 1.0);
  }
  else if (opt[0] == "rn")      // 反応有のコード指定
  { rvec[0] = opt[1]; }
  else if (opt[0] == "rp")      // 反応無のコード指定
  { rvec[1] = opt[1]; }
  else if (opt[0] == "res")     // 分析結果の出力指定
  {
    if (strsplit(reslst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "rmax")    // 反復計算の最大回数指定
  {
    char* ptr = nullptr;
    v_rmax = strtol(opt[1].c_str(), &ptr, 0);
    if (v_rmax < 0) err = badopt(pname(), optstr, USAGE, 0.0, (double)INT_MAX);
  }
  else if (opt[0] == "xno")     // X軸にする説明変数の番号指定
  {
    char* ptr = nullptr;
    v_xno = strtol(opt[1].c_str(), &ptr, 0);
    if (v_xno < 1) err = badopt(pname(), optstr, USAGE, 1.0, (double)INT_MAX);
  }
  else if (opt[0] == "xrange")  // X軸の最小値と最大値の指定
  {
    if (getrange(xrange.first, xrange.second, opt[1]) == 0) // 範囲指定無
    { err = badopt(pname(), optstr, USAGE, -FDmax, FDmax); }
    else if (xrange.first >= xrange.second) // 範囲指定不適切
    { err = badopt(pname(), optstr, USAGE, -FDmax, FDmax); }
  }
  else if (opt[0][0] == 'x')    // 説明変数の値指定
  {
    if (strsplit(xlst, optstr, CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else                          // 間違ったオプション指定
  { err = nonopt(pname(), optstr, USAGE); }
  if (err > 0) throw help_except("ヘルプ表示");
  if (err < 0) throw option_error(string("オプション指定違反→") + optstr);
}

int option::check(size_t ni)    // オプションのセットとチェック
{
  if (ni < 2) return disperr(pname(), "項目数が足りません！");
  if (rvec[0].empty() && rvec[1].empty()) return disperr(pname(), "反応有・無のコードが指定されていません！");
  int nx = (int)ni - 1;         // 交互作用項目を含まない説明変数の数
  if (nx < v_xno) return disperr(pname(), "X軸にする説明変数の番号が間違っています！");
  if (set_ngp() < 0) return -2; // -ngpオプションのセット
  if (set_res() < 0) return -3; // -resオプションのセット
  int nxix = set_ix(nx);        // -ixオプションのセット
  if (nxix < 0) return -4;      // 交互作用項目を含んだ説明変数の数
  if (set_x(nx) < 0) return -5;   // -xオプションのセット
  if (blst.empty()) return 0;     // 偏回帰係数初期値の指定なし
  int np = nxix + 1;              // パラメータ数
  if ((int)(blst.size()) > np) return disperr(pname(), "偏回帰係数初期値の値が多すぎます！");
  v_bvec.resize(np, BLANKf);
  int j = 0;
  for (const string& bstr : blst)
  {
    double fdata = strtodf(bstr);
    if (fdata == BLANKf) return disperr(pname(), "初期値指定が間違っています→", bstr);
    v_bvec[j] = fdata;
    ++j;
  }
  return 0;                     // エラーなし
}

int option::set_ngp()           // -ngpオプションのセット
{
  ngpmap["lci"] = ngpmap["log"] = ngpmap["or"] = ngpmap["score"] = true;
  if (ngplst.size() > 0)
  {
    if (v_ngpname == "") return disperr(pname(), "-ngpnameオプションの指定が必要！");
    ngpmap["lci"] = ngpmap["log"] = ngpmap["or"] = ngpmap["score"] = false;
    for (const string& ngpstr : ngplst)
    {
      mapsb_citr ngpitr = ngpmap.find(ngpstr);
      if (ngpitr == ngpmap.end())
      { return disperr(pname(), "-ngp=" + ngpstr + "オプション指定が不適切！"); }
      else
      { ngpmap[ngpstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["ltab"] = resmap["selx"] = resmap["stat"] = true;
  resmap["allx"] = resmap["v"] = false;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["ltab"] = resmap["selx"] = resmap["stat"] = false;
    for (const string& resstr : reslst)
    {
      mapsb_citr resitr = resmap.find(resstr);
      if (resstr == "-")           // 全て出力せず
      { resmap["res"] = false; }
      else if (resstr[0] == '-' && resstr.size() > 1) // 特定の項目を出力せず
      {
        resitr = resmap.find(resstr.substr(1));
        if (resitr == resmap.end())
        { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
        else
        { resmap[resstr.substr(1)] = false; }
      }
      else if (resstr == "def")    // 省略時出力
      {
        resmap["def"] = resmap["cc"] = resmap["ltab"] = resmap["selx"] = resmap["stat"] = true;
      }
      else if (resstr == "all")    // 全て出力
      {
        resmap["def"] = resmap["allx"] = resmap["ltab"] = resmap["selx"] = resmap["stat"] = resmap["v"] = true;
      }
      else if (resitr == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

int option::set_ix(int nx)      // -ixオプションのセット
{
  if (ixstrlst.empty()) return nx;
  int ixno = nx + 1;            // 交互作用項目の最初の番号
  for (const string& ixstr : ixstrlst)
  {
    string opt[2];
    if (getopt(opt, ixstr, 'x') > 0)
    {
      int low0 = 0, up0 = 0, low1 = 0, up1 = 0;
      if (getrange(low0, up0, opt[0]) < 1)
      { return disperr(pname(), "1番目の説明変数番号が間違っています→", ixstr); }
      else
      {
        if (getrange(low1, up1, opt[1]) < 1)
        { return disperr(pname(), "2番目の説明変数番号が間違っています→", ixstr); }
        else
        {
          if (low0 < 1 || low0 > up0) return disperr(pname(), "1番目の説明変数番号が間違っています→", ixstr);
          if (low1 < 1 || low1 > up1) return disperr(pname(), "2番目の説明変数番号が間違っています→", ixstr);
          if (ixno <= up0 || ixno <= up1) return disperr(pname(), "説明変数番号が間違っています→", ixstr);
          for (int j = low1; j <= up1; ++j) // 2番目の番号
          {
            for (int i = low0; i <= up0; ++i) // 1番目の番号
            {
              pairint ixpi {i, j};
              v_ixlst.push_back(ixpi);
            }
          }
          ixno += (up0 - low0 + 1) * (up1 - low1 + 1);
        }
      }
    }
  }
  return ixno - 1;              // 交互作用項目を含めた説明変数の数
}

int option::set_x(int nx)       // -xオプションのセット
{
  if (xlst.empty()) return 0;   // 説明変数指定なし
  if ((int)(xlst.size()) > nx) return disperr(pname(), "説明変数の値が多すぎます！");
  xvec.resize(nx, BLANKf);
  for (const string& xstr : xlst)
  {
    string opt[2];
    if (getopt(opt, xstr, OPTdev) > 0)
    {
      char* ptr = nullptr;
      int xno = strtol(opt[0].substr(1).c_str(), &ptr, 0);
      if (xno < 1 || nx < xno || xno == v_xno) return disperr(pname(), "説明変数指定の項目番号が間違っています→", xstr);
      double fdata = strtodf(opt[1]);
      if (fdata == BLANKf) return disperr(pname(), "説明変数指定の値が間違っています→", xstr);
      xvec[xno - 1] = fdata;
    }
  }
  return (int)(xlst.size());    // 説明変数指定の数を返す
}

class data_less                 // ケースデータの選択済み説明変数比較用関数オブジェクト
{
private:
  option *optptr;               // オプションへのポインター
  const vecint *checkvecptr;    // 変数選択ベクトルへのポインター
  const vecdbl *mevecptr;       // 平均値ベクトルへのポインター
  int ni;                       // 項目数
public:
  data_less(option *optptr0, int ni0, const vecint *checkvecptr0, const vecdbl *mevecptr0) : optptr(optptr0), ni(ni0), checkvecptr(checkvecptr0), mevecptr(mevecptr0)
  // コンストラクタ
  {}
  ~data_less()                  // デストラクタ
  {}
  bool operator()(datac& cdata1, datac& cdata2) const; // ケースデータを比較する
};

bool data_less::operator()(datac& cdata1, datac& cdata2) const
// ケースデータの選択済み説明変数を比較する
{
  size_t dim = checkvecptr->size();
  size_t nxix = dim - 2;
  vecdbl fdatavec1(dim, 0.0), fdatavec2(dim, 0.0);
  setixdata(fdatavec1, optptr->ixlst(), optptr->ixc(), *mevecptr, cdata1);
  setixdata(fdatavec2, optptr->ixlst(), optptr->ixc(), *mevecptr, cdata2);
  int comp = 0;
  for (size_t i = 1; comp == 0 && i <= nxix; ++i)
  {
    if (checkvecptr->at(i))     // 選択済み説明変数だけ比較する
    {
      double fdata1 = fdatavec1[i], fdata2 = fdatavec2[i];
      if (fdata1 == BLANKf)
      {
        if (fdata2 != BLANKf) comp = 1;
      } else if (fdata2 == BLANKf)
      { comp = -1; }
      else
      { comp = sgndiff(fdata1, fdata2); }
    }
  }
  return (comp < 0) ? true : false; // 小さいときだけtrueを返す
}

class varclass                 // 変数クラス
{
private:
  mat<double> v_smat;           // 積和行列
  vecdbl midvec;                // 反応有無別平均値の中央値ベクトル
  vecdbl v_mevec;               // 平均値ベクトル
  vecdbl v_vrvec;               // 分散ベクトル
  vecint v_checkvec;            // 変数選択用ベクトル
  vec<double> b0v;              // 偏回帰係数初期値ベクトル
  vec<double> bv;               // 偏回帰係数ベクトル
  vec<double> vbv;              // 偏回帰係数の分散ベクトル
  double v_lb;                  // 対数尤度
  int v_nx;                     // 説明変数の数
  int v_nxix;                   // 交互作用項目を含んだ説明変数の数
  int v_n;                      // 例数数
  int v_nrp;                    // 反応有の例数
  int v_nrep;                   // 反復計算回数
  restype v_resvar;             // 処理結果の種類
  int set_mean(int& nr, vecdbl& mevec, option& opt, ldata& datalst) const;
  // 説明変数の平均値ベクトルをセットする
  restype cal_csqr(vecdbl& cvec, const vecdbl& mevec, const vec<double>& b0v, option& opt, ldata& datalst) const;
  // 変数選択用χ2乗値を計算する
  restype cal_lgh(vec<double>& gv, mat<double>& hmat, double& lb, const vecint& checkvec, const vec<double>& bv, const vecdbl& mevec, option opt, ldata& datalst) const;
  // 傾斜ベクトルgとHessian行列Hを計算する
  int sel_invar(vecint& checkvec, vec<double>& btv, vec<double>& vbtv, const vecdbl& cvec, const vec<double>& bv) const;
  // 取り込み変数の選択と各種設定
  int sel_outvar(double cin, double pin, vecint& checkvec, vec<double>& bv, vec<double>& vbv, vec<double>& btv, vec<double>& vbtv) const;
  // 係数ベクトルの代入と追い出し変数の選択
  restype newton_log(int& nrep, vec<double>& bv, vec<double>& vbv, double& lb, const vecdbl& mevec, const vecint& checkvec, int rmax, option& opt, ldata& datalst) const;
  // Newton法によるロジットモデルの計算(最尤法)
  void set_initb(vec<double>& b0v, vec<double>& vbv, option& opt) const;
  // 重回帰分析→判別分析による偏回帰係数初期値のセット
public:
  varclass(size_t ni, option& opt); // コンストラクタ
  ~varclass()                   // デストラクタ
  {}
  const mat<double>& smat() const // 積和行列を返す
  { return v_smat; }
  const vecdbl& mevec() const   // 平均値ベクトルを返す
  { return v_mevec; }
  const vecdbl& vrvec() const   // 分散ベクトルを返す
  { return v_vrvec; }
  const vecint& checkvec() const   // 変数選択用を返す
  { return v_checkvec; }
  double mid(int idx) const     // 反応有無別平均値の中央値を返す
  { return (0 <= idx && idx <= np()) ? midvec[idx] : 0.0; }
  double me(int idx) const      // 平均値を返す
  { return (0 <= idx && idx <= np()) ? v_mevec[idx] : 0.0; }
  int check(int idx) const      // 変数選択フラッグを返す
  { return (0 <= idx && idx < np()) ? v_checkvec[idx] : 0; }
  double b0(int idx) const      // 偏回帰係数初期値を返す
  { return (0 <= idx && idx < np()) ? b0v(idx) : 0.0; }
  double b(int idx) const       // 偏回帰係数を返す
  { return (0 <= idx && idx < np()) ? bv(idx) : 0.0; }
  double vb(int idx) const      // 偏回帰係数の分散を返す
  { return (0 <= idx && idx < np()) ? vbv(idx) : 0.0; }
  double bst(int idx) const     // 標準偏回帰係数を返す
  { return (0 <= idx && idx < np()) ? b(idx) * sqrt(v_vrvec[idx]) : 0.0; }
  double lb() const             // 対数尤度を返す
  { return v_lb; }
  int nx() const                // 説明変数の数を返す
  { return v_nx; }
  int nxix() const              // 交互作用項目を含んだ説明変数の数を返す
  { return v_nxix; }
  int np() const                // パラメータ数を返す
  { return nxix() + 1; }
  int n() const                 // 例数を返す
  { return v_n; }
  int nrp() const               // 反応有の例数を返す
  { return v_nrp; }
  int nrn() const               // 反応無の例数を返す
  { return n() - nrp(); }
  int dfb() const               // 回帰自由度を返す
  { return v_checkvec[np()] - 1; }
  int nrep() const              // 反復計算回数を返す
  { return v_nrep; }
  restype resvar() const        // 処理結果の種類を返す
  { return v_resvar; }
  void cal(option& opt, ldata& datalst); // データを集計して各種統計量を求める
  void stepsel(option& opt, ldata& datalst); // 増減法による変数選択
  void allsel(option& opt, ldata& datalst);  // 全変数選択法
  double lscore(double& z, option& opt, datac& cdata) const;
  // ロジットスコアを計算してyを返す
};

varclass::varclass(size_t ni, option& opt) : v_smat(), midvec(), v_mevec(), v_vrvec(), v_checkvec(), b0v(), bv(), vbv(), v_lb(0.0), v_nx(0), v_nxix(0), v_n(0), v_nrp(0), v_nrep(0), v_resvar(resnon)
// コンストラクタ
{
  v_nx = ni - 1;
  v_nxix = v_nx + (opt.ixlst()).size();
  size_t dim = v_nxix + 2;
  v_smat.resize(dim, dim);
  midvec.resize(dim, 0.0);  v_mevec.resize(dim, 0.0);  v_vrvec.resize(dim, 0.0);
  v_checkvec.resize(dim, 0);
  b0v.resize(dim -1);  bv.resize(dim -1);  vbv.resize(dim -1);
}

void varclass::cal(option& opt, ldata& datalst)
// データを集計して各種統計量を求める
{
  vecdbl mevec(nx() + 1);       // 交互作用項計算用平均値ベクトル
  v_n = set_mean(v_nrp, mevec, opt, datalst);
  if (v_n < np() || v_nrp < 1 || v_nrp >= v_n) // ロジスティック関数計算不可能
  {
    v_resvar = reserr;
    return;
  }
  int dim = np() + 1;
  vecdbl mevec0(dim, 0.0), mevec1(dim, 0.0);
  vecdbl fdatavec(dim, 0.0);
  for (auto& cdatai : datalst)
  {
    if (cdatai.sdata(0) == "0" || cdatai.sdata(0) == "1") // 有効なケース
    {
      setixdata(fdatavec, opt.ixlst(), opt.ixc(), mevec, cdatai);
      for (int i = 0; i < dim; ++i)
      {
        double x = fdatavec[i];
        for (int j = i; j < dim; ++j) v_smat.add(i, j, x * fdatavec[j]);
        if (cdatai.sdata(0) == "0") mevec0[i] += x;
        else mevec1[i] += x;
      }
    }
  }
  for (int j = 1; j < dim; ++j) // 反応の有無別平均値とその中央値の計算
  {
    mevec0[j] /= mevec0[0];  mevec1[j] /= mevec1[0];
    midvec[j] = (mevec0[j] + mevec1[j]) * 0.5;
  }
  for (int i = 0; i < np(); ++i)
  {
    for (int j = i + 1; j < dim; ++j) v_smat.assign(j, i, v_smat(i, j));
  }
  v_smat.sweep(0);              // 単純積和行列→積和行列
  int dft = v_n - 1;
  for (int i = 0; i < dim; ++i)
  {
    v_mevec[i] = v_smat(0, i);
    v_vrvec[i] = (dft > 0) ? v_smat(i, i) / (double)dft : 0.0;
  }
  v_resvar = resok;
}

int varclass::set_mean(int& nrp, vecdbl& mevec, option& opt, ldata& datalst) const
// 説明変数の平均値ベクトルをセットする
{
  int ni = nx() + 1;            // 項目数
  vecdbl fdatavec(ni, 0.0);
  int n = 0;
  nrp = 0;
  for (auto& cdatai : datalst)
  {
    if (cdatai.sdata(0) != BLANKs)
    {
      double y = -1.0;
      if (opt.rn().empty())     // 反応無しコード指定なし
      { y = (cdatai.sdata(0) == opt.rp()) ? 1.0 : 0.0; }
      else if (cdatai.sdata(0) == opt.rn())
      { y = 0.0; }
      else if (opt.rp().empty() || cdatai.sdata(0) == opt.rp())
      { y = 1.0; }
      if (y < 0.0)              // 指定コード以外
      { cdatai.assign(0, BLANKs); } // 集計から除外する
      else
      {
        bool blank = false;
        int j = 1;
        while (blank == false && j < ni)
        {
          fdatavec[j] = cdatai.fdata(j);
          if (fdatavec[j] == BLANKf) blank = true;
          else ++j;
        }
        if (blank)
        { cdatai.assign(0, y - 2.0); } // 欠損補完のためにコードを-2/-1にする
        else
        {
          ++n;
          if (y == 1.0) ++nrp;  // 反応有の例数を数える
          cdatai.assign(0, y);  // 後の集計のためにコードを0/1にする
          fdatavec[0] = y;
          for (j = 0; j < ni; ++j) mevec[j] += fdatavec[j];
        }
      }
    }
  }
  for (int j = 0; j < ni; ++j) mevec[j] /= (double)n;
  return n;                     // 有効な例数を返す
}

void varclass::stepsel(option& opt, ldata& datalst)
// 増減法による変数選択
{
  v_checkvec[0] = 1;  v_checkvec[np()] = 1; // 切片は無条件に取り込み
  vecdbl cvec(np(), 0.0);
  vec<double> btv(np()), vbtv(np());
  int iin = 0, iout = 0, nrep = 0;
  restype resnew = resnon;
  do {
    if (cal_csqr(cvec, v_mevec, bv, opt, datalst) == reserr)
    { resnew = reserr; }
    else
    {
      iin = sel_invar(v_checkvec, btv, vbtv, cvec, bv);
      if (iin > 0)                // 取り込み変数あり
      {
        iout = 0;
        int i = 0;
        do {
          resnew = newton_log(nrep, btv, vbtv, v_lb, v_mevec, v_checkvec, opt.rmax(), opt, datalst);
          if (nrep > v_nrep) v_nrep = nrep;
          i = sel_outvar(opt.cin(), opt.pin(), v_checkvec, bv, vbv, btv, vbtv);
          if (i > 0) iout = i;    // 最後に追い出した変数番号を記憶
        } while (resnew != reserr && i > 0); // 追い出し変数が無くなるまで続行
      }
    }
  } while (resnew != reserr && iin > 0 && iin != iout && v_checkvec[np()] < np());
  v_resvar = resnew;
}

restype varclass::cal_csqr(vecdbl& cvec, const vecdbl& mevec, const vec<double>& b0v, option& opt, ldata& datalst) const
// 変数選択用χ2乗値を計算する
{
  int dim = np() + 1;
  vec<double> bv(b0v), gv(np());
  mat<double> hmat(np());
  double lb = 0.0;
  vecint checkvec(dim, 1);
  checkvec[np()] = np();
  if (cal_lgh(gv, hmat, lb, checkvec, bv, mevec, opt, datalst) == reserr) return reserr;
  hmat.inv();
  bv -= hmat * gv;
  for (int i = 0; i < np(); ++i) cvec[i] = (hmat(i, i) < 0.0) ? bv(i) * bv(i) / (-hmat(i, i)) : 0.0;
  return resok;
}

restype varclass::cal_lgh(vec<double>& gv, mat<double>& hmat, double& lb, const vecint& checkvec, const vec<double>& bv, const vecdbl& mevec, option opt, ldata& datalst) const
// 傾斜ベクトルgとHessian行列Hを計算する
{
  constexpr double Emax = 308.0; // 指数最大値
  int dim = np() + 1;
  int npt = checkvec[np()];
  double sbxy = 0.0, sei = 0.0;
  vec<double> sxyv(npt), sxv(npt), xyv(npt + 1);
  xyv.assign(0, 1.0);  gv = 0.0;  hmat = 0.0;  lb = 0.0;
  vecdbl fdatavec(dim, 0.0);
  restype reslgh = resnon;
  for (auto& cdatai : datalst)
  {
    if (cdatai.sdata(0) == "0" || cdatai.sdata(0) == "1")
    {
      setixdata(fdatavec, opt.ixlst(), opt.ixc(), mevec, cdatai);
      double ei = bv(0);
      int j = 1;
      for (int i = 1; i < np(); ++i)
      {
        if (checkvec[i])
        {
          xyv.assign(j, fdatavec[i]);
          ei += xyv(j) * bv(j);
          ++j;
        }
      }
      if (fabs(ei) > Emax)        // 指数最大値以上
      { reslgh = reserr; }
      else
      {
        if (cdatai.sdata(0) == "0") xyv.assign(npt, 0.0);
        else { xyv.assign(npt, 1.0);  sbxy += ei; }
        ei = exp(ei);
        sei += log(1.0 + ei);
        double wi = ei / (1.0 + ei), wi2 = ei / (1.0 + ei) / (1.0 + ei);
        for (j = 0; j < npt; ++j)
        {
          sxyv.add(j, xyv(j) * xyv(npt));
          sxv.add(j, xyv(j) * wi);
          for (int k = j; k < npt; ++k) hmat.add(j, k, -xyv(j) * xyv(k) * wi2);
        }
      }
    }
  }
  if (reslgh == reserr) return reslgh; // 計算エラー
  lb = sbxy - sei;  gv = sxyv - sxv;
  for (int i = 0; i < npt - 1; ++i)
  {
    for (int j = i + 1; j < npt; ++j) hmat.assign(j, i, hmat(i, j));
  }
  return resok;
}

int varclass::sel_invar(vecint& checkvec, vec<double>& btv, vec<double>& vbtv, const vecdbl& cvec, const vec<double>& bv) const
// 取り込み変数の選択と各種設定
{
  double cmax = 0.0;
  int imax = 0;
  for (int i = 1; i < np(); ++i)
  {
    if (checkvec[i] == 0 && cvec[i] > cmax) { cmax = cvec[i];  imax = i; }
  }
  if (imax < 1) return imax;    // 取り込み変数なし
  checkvec[imax] = 1;  checkvec[np()] += 1;
  btv.resize(checkvec[np()]);  vbtv.resize(checkvec[np()]);
  btv = 0.0;
  return imax;                  // 変数番号を返す
}

int varclass::sel_outvar(double cin, double pin, vecint& checkvec, vec<double>& bv, vec<double>& vbv, vec<double>& btv, vec<double>& vbtv) const
// 係数ベクトルの代入と追い出し変数の選択
{
  double cmin = cin, pmax = pin;
  int imin = 0, j = 1;
  bool csel = (pin < 1.0) ? false : true; // χ2乗値による変数選択
  bv.assign(0, btv(0));  vbv.assign(0, vbtv(0));
  for (int i = 1; i < np(); ++i)
  {
    if (checkvec[i])
    {
      bv.assign(i, btv(j));  vbv.assign(i, vbtv(j));
      double csqr = (vbtv(j) > 0.0) ? btv(j) * btv(j) / vbtv(j) : 0.0;
      if (csel)
      {
        if (csqr <= cmin) { cmin = csqr;  imin = i; }
      }
      else
      {
        double pr = ctop(csqr, 1);
        if (pr >= pmax) { pmax = pr;  imin = i; }
      }
      ++j;
    }
  }
  if (imin < 1) return imin;    // 追い出し変数なし
  checkvec[imin] = 0;  checkvec[np()] -= 1;
  btv.resize(checkvec[np()]);  vbtv.resize(checkvec[np()]);
  btv = 0.0;
  return imin;                  // 変数番号を返す
}

restype varclass::newton_log(int& nrep, vec<double>& bv, vec<double>& vbv, double& lb, const vecdbl& mevec, const vecint& checkvec, int rmax, option& opt, ldata& datalst) const
// Newton法によるロジットモデルの計算(最尤法)
{
  int npt = checkvec[np()];      // 一時的パラメータ数
  vec<double> gv(npt);
  mat<double> hmat(npt);
  double lold = 0.0;
  restype resnew = resnon;
  nrep = 0;
  while (resnew == resnon && cal_lgh(gv, hmat, lb, checkvec, bv, mevec, opt, datalst) != reserr)
  {
    ++nrep;
    hmat.inv();
    if (fabs(lb - lold) < opt.eps()) resnew = resok;
    else if (rmax > 0 && nrep >= rmax) resnew = resrmax;
    else { lold = lb;  bv -= hmat * gv; }
  }
  if (resnew == resnon) return reserr; // 発散！
  for (int i = 0; i < npt; ++i) vbv.assign(i, -hmat(i, i));
  return resnew;                // Newton法成功
}

void varclass::allsel(option& opt, ldata& datalst)
// 全変数選択法
{
  set_initb(b0v, vbv, opt);     // 全変数選択時の初期値設定
  bv = b0v;
  v_checkvec.assign(np() + 1, 1); // 全変数選択
  v_checkvec[np()] = np();
  v_resvar = newton_log(v_nrep, bv, vbv, v_lb, v_mevec, v_checkvec, opt.rmax(), opt, datalst);
  if (v_resvar == reserr)       // Newton法が発散した時
  {
    bv = b0v;                   // 初期値を解にする
    v_resvar = newton_log(v_nrep, bv, vbv, v_lb, v_mevec, v_checkvec, 1, opt, datalst);
  }
}

void varclass::set_initb(vec<double>& b0v, vec<double>& vbv, option& opt) const
// 重回帰分析→判別分析による偏回帰係数初期値のセット
{
  mat<double> smat(v_smat);
  double ym = smat(0, np()), st = smat(np(), np());
  for (int i = 1; i < np(); ++i) smat.sweep(i);
  double w = 1.0, sr = smat(np(), np());
  if (sr > 0.0) w = (n() - 2) / sr;
  int dfr = n() - np();
  double vr = (dfr > 0) ? sr / (double)dfr : 0.0;
  for (int i = 1; i < np(); ++i)
  {
    double bi = smat(i, np());
    if (fabs(bi) < FDeps) bi = 0.0;
    b0v.assign(i, bi * w);
    double vbi = vr * smat(i, i);
    vbv.assign(i, vbi * w * w);
  }
  b0v.assign(0, (smat(0, np()) - ym + (ym - 0.5) * (st - sr) / st) * w - log((1.0 - ym) / ym));
  // (1-π)/πの補正
  if (opt.bvec().empty()) return; // 初期値指定なし
  for (int i = 0; i < np(); ++i)
  {
    if (opt.b(i) != BLANKf) b0v.assign(i, opt.b(i));
  }
}

double varclass::lscore(double& z, option& opt, datac& cdata) const
// ロジットスコアを計算してyを返す
{
  double y = cdata.fdata(0);
  if (y == BLANKf) return BLANKf; // 目的変数欠損のため計算できず
  int dim = np() + 1;
  vecdbl fdatavec(dim, 0.0);
  setixdata(fdatavec, opt.ixlst(), opt.ixc(), v_mevec, cdata);
  z = b(0);
  bool blank = false;
  for (int i = 1; blank == false && i < np(); ++i)
  {
    if (fdatavec[i] == BLANKf)
    {
      if (opt.comp())           // 欠損を平均値で補完する場合
      {
        if (v_checkvec[i]) z += b(i) * midvec[i];
      }
      else
      { blank = true; }
    }
    else
    {
      if (v_checkvec[i]) z += b(i) * fdatavec[i];
    }
  }
  if (y < 0.0) y += 2.0;        // データに欠損がある症例の目的変数
  return blank ? BLANKf : y;    // 計算できたらyを返す
}

class logregclass               // logregメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout, option& opt, const varclass& var);
  // 解析結果のタイトルとデータ識別情報を出力する
  int set_itemix(idd& tid, option& opt, int ni) const;
  // 交互作用項目名をセットする
  void put_stat(ostream& fout, option& opt, const varclass& var) const;
  // 統計量を出力する
  void put_logreg(ostream& fout, option& opt, double cin, double pin, const varclass& var);
  // ロジスティック回帰分析結果を出力する
  void put_brci(ostream& fout, option& opt, const varclass& var) const;
  // 偏回帰係数とオッズ比の信頼区間を出力する
  void put_logtab(ostream& fout, option& opt, const varclass& var);
  // ロジスティック回帰分析の分散分析表を出力する
  int cal_loglf(double& lf, option& opt, const varclass& var);
  // 飽和モデルの尤度を求める
  bool data_equal(option& opt, const vecint& checkvec, const vecdbl& mevec, datac& cdata1, datac& cdata2) const;
  // ケースデータの選択済み説明変数を比較する
  void put_dns(option& opt, const varclass& var);
  // データファイル形式のロジットスコアと反応有の確率を出力する
  void plot_logngp(option& opt, double cin, double pin, const varclass& var);
  // Ngraph用データ出力
  void plot_ngpid(ostream& fout, option& opt, double cin, double pin, const varclass& var);
  // Ngraph用データのタイトルとデータ識別情報を出力する
  void plot_log(ostream& fout, option& opt, const varclass& var);
  // Ngraph用ロジスティック曲線を出力する
  double cal_range(double& xmin, double& xmax, option& opt, datac& cdata, const varclass& var);
  // ロジスティック曲線のxの最小値と最大値を計算する
  void plot_score(ostream& fout, option& opt, const varclass& var);
  // Ngraph用ロジットスコアを出力する
  void plot_lci(ostream& fout, option& opt, const varclass& var);
  // Ngraph用ロジスティック曲線の信頼区間を出力する
  void plot_or(ostream& fout, option& opt, const varclass& var) const;
  // Ngraph用にORとその信頼区間を出力する
public:
  logregclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~logregclass()                // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み累積する
  void put(option& opt);                     // 集計結果を出力する
};

void logregclass::put_resid(ostream& fout, option& opt, const varclass& var)
// 解析結果のタイトルとデータ識別情報を出力する
{
  puttitle(fout, Title);
  fout << "\n";
  id.put_title(fout);
  id.put_stra(fout);
  idd tid(id);
  tid.set(id.title(), id.ni() + (opt.ixlst()).size()); // 項目数を変更する
  set_itemix(tid, opt, id.ni());
  fout << "目的変数y   :" << tid.item(0) << "\n";
  for (size_t i = 1; i < tid.ni(); ++i)
  {
    fout << "説明変数x" << setw(3) << i;
    fout << ":" << tid.item(i);
    if (opt.comp()) fout << " 欠損補完用平均値=" << var.mid(i);
    fout << "\n";
  }
}

int logregclass::set_itemix(idd& tid, option& opt, int ni) const
// 交互作用項目名をセットする
{
  ostringstream ostr;
  int j = ni;
  for (const auto& ixi : opt.ixlst())
  {
    ostr.str("");
    if (opt.ixc())              // 偏差変換
    {
      if (ixi.first < ni) ostr << tid.item(ixi.first) << "偏差 x ";
      else ostr << tid.item(ixi.first) << " x ";
      if (ixi.second < ni) ostr << tid.item(ixi.second) << "偏差";
      else ostr << tid.item(ixi.second);
    }
    else
    { ostr << tid.item(ixi.first) << " x " << tid.item(ixi.second); }
    tid.assign_item(j, ostr.str());
    ++j;
  }
  return ni - j;                // セットした交互作用項目名数を返す
}

void logregclass::put_stat(ostream& fout, option& opt, const varclass& var) const
// 統計量を出力する
{
  putmsdvec(fout, var.n(), var.mevec(), var.vrvec(), var.nxix(), 1);
  fout << "反応有:コード=";
  if (opt.rp().empty()) fout << opt.rn() << "以外";
  else fout << opt.rp();
  fout << " 例数=" << var.nrp();
  fout << "  反応無:コード=";
  if (opt.rn().empty()) fout << opt.rp() << "以外";
  else fout << opt.rn();
  fout << " 例数=" << var.nrn() << "\n";
  putmat(fout, "・相関行列(correlation matrix)", matr, var.n(), var.smat(), var.nxix(), 1);
}

void logregclass::put_logreg(ostream& fout, option& opt, double cin, double pin, const varclass& var)
// ロジスティック回帰分析結果を出力する
{
  if (cin <= 0.0 && pin >= 1.0)
  {
    fout << "\n・全変数を選択した結果  反復計算回数=" << var.nrep();
    if (var.resvar() == reserr)
    { fout << ":Newton法発散のため初期値を使用\n"; }
    else if (var.resvar() == resrmax)
    {
      if (var.nrep() == 1) fout << ":初期値を使用\n";
      else fout << ":最大反復回数\n";
    }
    else
    { fout << "\n"; }
  }
  else
  {
    fout << "\n・前進的変数増減法(stepwise forward selection method)による変数選択結果\n";
    if (pin < 1.0)
    {
      fout << "  取り込み基準:p値＜" << pin;
      fout << "  追い出し基準:p値≧" << pin;
    }
    else
    {
      fout << "  取り込み基準:χ2乗値＞" << cin;
      fout << "  追い出し基準:χ2乗値≦" << cin << "\n";
    }
    fout << "  反復計算回数=" << var.nrep();
    if (var.resvar() == reserr)
    { fout << ":Newton法発散！\n"; }
    else if (var.resvar() == resrmax)
    {
      if (var.nrep() == 1) fout << ":初期値を使用\n";
      else fout << ":最大反復回数\n";
    }
    else
    { fout << "\n"; }
  }
  fout << "  ロジットモデル:p=1/{1+exp(-β0-Σβj・xj)}\n";
  fout << "  p:y=";
  if (opt.rp().empty()) fout << opt.rn() << " 以外";
  else fout << opt.rp();
  fout << "(反応有)の確率  β0:切片  βj:変数xjの偏回帰係数\n";
  if (opt.res("v"))
  {
    fout << "  偏回帰係数初期値:";
    for (int i = 0; i < var.np(); ++i)
    {
      if (i > 0)
      {
        if (i % 4 == 0) fout << "\n" << setw(19) << " ";
        else fout << " ";
      }
      fout << "β" << i << "=" << var.b0(i);
    }
    fout << "\n";
  }
  fout << "                                            標準                   有意確率\n";
  fout << "変数  偏回帰係数   標準誤差    オッズ比  偏回帰係数  Waldのχ2乗       p値\n";
  constexpr int LENline = 78;   // 区切り線の長さ
  putline(fout, LENline);
  double csqr = (var.vb(0) > 0.0) ? var.b(0) * var.b(0) / var.vb(0) : 0.0;
  ldbl pr = ctop(csqr, 1);
  fout << "切片" << setw(12) << var.b(0) << setw(11) << sqrt(var.vb(0));
  fout << setw(36) << csqr << setw(12) << pr << pind(pr) << "\n";
  for (int i = 1; i < var.np(); ++i)
  {
    if (var.check(i))
    {
      csqr = (var.vb(i) > 0.0) ? var.b(i) * var.b(i) / var.vb(i) : 0.0;
      pr = ctop(csqr, 1);
      fout << "x" << setw(3) << i << setw(12) << var.b(i);
      fout << setw(11) << sqrt(var.vb(i)) << setw(12) << exp(var.b(i));
      fout << setw(12) << var.bst(i) << setw(12) << csqr;
      fout << setw(12) << pr << pind(pr) << "\n";
    }
  }
  putline(fout, LENline);
  put_brci(fout, opt, var);
  if (opt.res("ltab")) put_logtab(fout, opt, var);
}

void logregclass::put_brci(ostream& fout, option& opt, const varclass& var) const
// 偏回帰係数とオッズ比の信頼区間を出力する
{
  fout << "\n変数  偏回帰係数" << setw(5) << opt.cc() * 100.0;
  fout << "%CI下限        上限    オッズ比" << setw(5) << opt.cc() * 100.0;
  fout << "%CI下限        上限\n";
  constexpr int LENline = 76;   // 区切り線の長さ
  putline(fout, LENline);
  double z = ptoz(1.0 - opt.cc());
  double ci = z * sqrt(var.vb(0));
  double bl = var.b(0) - ci, bu = var.b(0) + ci;
  fout << "切片" << setw(12) << var.b(0) << setw(12) << bl;
  fout << setw(12) << bu << "\n";
  for (int i = 1; i < var.np(); ++i)
  {
    if (var.check(i))
    {
      ci = z * sqrt(var.vb(i));
      bl = var.b(i) - ci;  bu = var.b(i) + ci;
      fout << "x" << setw(3) << i << setw(12) << var.b(i);
      fout << setw(12) << bl << setw(12) << bu;
      fout << setw(12) << exp(var.b(i)) << setw(12) << exp(bl);
      fout << setw(12) << exp(bu) << "\n";
    }
  }
  putline(fout, LENline);
}

void logregclass::put_logtab(ostream& fout, option& opt, const varclass& var)
// ロジスティック回帰分析の分散分析表を出力する
{
  int dft = var.n() - 1;
  double lf = 0.0;
  dft -= cal_loglf(lf, opt, var);
  int dfl = dft - var.dfb();
  double l0 = (double)(var.nrp()) * log((double)(var.nrp()) / (double)(var.n())) + (double)(var.nrn()) * log((double)(var.nrn()) / (double)var.n());
  double r2 = 0.0, lbdif = var.lb() - l0;
  if (sgndiff(lf, l0) > 0) r2 = lbdif / (lf - l0);
  fout << "対数尤度:L(モデル)=" << var.lb() << "  L(切片)=" << l0 << "  飽和Lf=" << lf << "\n";
  fout << "擬似寄与率R2乗=" << r2 << "  AIC(赤池の情報量基準)=" << 2.0 * ((double)(var.dfb()) + 1.0 - var.lb()) << "\n";
  const string Head = "回帰とズレの検定";
  int w = strclen(Head, CHARwide);
  int npad = Head.size() - w;
  fout << "\n" << setw((68 + w) / 2 + npad) << Head;
  fout << "\n要因         -(対数尤度比)    自由度        χ2乗値   有意確率p値\n";
  constexpr int LENline = 68;
  putline(fout, LENline);
  double csqr = 2.0 * lbdif;
  ldbl pr = ctop(csqr, var.dfb());
  fout << "回帰" << setw(22) << lbdif << setw(10) << var.dfb();
  fout << setw(14) << csqr << setw(14) << pr << pind(pr) << "\n";
  lbdif = lf - var.lb();
  csqr = 2.0 * lbdif;
  if (dfl < 1 || csqr < FDeps) csqr = 0.0;
  pr = ctop(csqr, dfl);
  fout << "ズレ(LOF)" << setw(17) << lbdif << setw(10) << dfl;
  fout << setw(14) << csqr << setw(14) << pr << pind(pr) << "\n";
  putline(fout, LENline);
  fout << "全体" << setw(22) << lf - l0 << setw(10) << dft << endl;
}

int logregclass::cal_loglf(double& lf, option& opt, const varclass& var)
// 飽和モデルの尤度を求める
{
  vecint checkvec(var.checkvec());
  vecdbl mevec(var.mevec());
  datalst.sort(data_less(&opt, id.ni(), &checkvec, &mevec)); // 選択済み説明変数順にソートする
  int nlf = 0;
  lf = 0.0;
   ldata_itr iitr = datalst.begin();
   while (iitr != datalst.end())
  {
    if (iitr->sdata(0) != "0" && iitr->sdata(0) != "1")
    { ++iitr; }
    else
    {
      int ni = 1, ri = 0;
      if (iitr->sdata(0) == "1") ri = 1; // 反応有
      bool equal = true;
      ldata_itr jitr = iitr;
      ++jitr;
      while (equal && jitr != datalst.end())
      {
        if (jitr->sdata(0) != "0" && jitr->sdata(0) != "1")
        { ++jitr; }
        else
        {
          equal = data_equal(opt, checkvec, mevec, *iitr, *jitr);
          if (equal)
          {
            if (jitr->sdata(0) == "1") ++ri;
            ++ni;  ++jitr;
          }
        }
      }
      if (ni > 1)
      {
        nlf += ni - 1;
        if (ri > 0 && ri < ni) lf += (double)ri * log((double)ri / (double)ni) + (double)(ni - ri) * log((double)(ni - ri) / (double)ni);
      }
      iitr = jitr;
    }
  }
  return nlf;                   // 飽和モデルの重複数を返す(自由度減少数)
}

bool logregclass::data_equal(option& opt, const vecint& checkvec, const vecdbl& mevec, datac& cdata1, datac& cdata2) const
// ケースデータの選択済み説明変数を比較する
{
  size_t dim = checkvec.size();
  size_t nxix = dim - 2;
  vecdbl fdatavec1(dim, 0.0), fdatavec2(dim, 0.0);
  setixdata(fdatavec1, opt.ixlst(), opt.ixc(), mevec, cdata1);
  setixdata(fdatavec2, opt.ixlst(), opt.ixc(), mevec, cdata2);
  bool equal = true;
  for (size_t i = 1; equal && i <= nxix; ++i)
  {
    if (checkvec[i])            // 選択済み説明変数だけ比較する
    {
      double fdata1 = fdatavec1[i], fdata2 = fdatavec2[i];
      if (fdata1 == BLANKf)
      {
        if (fdata2 != BLANKf) equal = false;
      } else if (fdata2 == BLANKf)
      { equal = false; }
      else
      {
        if (sgndiff(fdata1, fdata2) != 0) equal = false;
      }
    }
  }
  return equal;                 // true:一致, false:不一致
}

void logregclass::put_dns(option& opt, const varclass& var)
// データファイル形式のロジットスコアと反応有の確率を出力する
{
  ofstream fout;
  if (ofopen(fout, (opt.dnsname()).c_str()) == false) return; // ファイルオープン失敗
  idd tid(id.title(), 3);
  tid.assign_item(0, id.item(0));
  if (opt.cin() > 0.0)
  {
    tid.assign_item(1, "ロジットスコア(変数選択法)");
    tid.assign_item(2, "反応有の確率(変数選択法)");
  }
  else
  {
    tid.assign_item(1, "ロジットスコア(全変数選択)");
    tid.assign_item(2, "反応有の確率(全変数選択)");
  }
  id.put_title(fout, formtag);
  id.put_stra(fout, formtag);
  tid.put_item(fout, formtag);
  datac cdata(3);
  for (auto& cdatai : datalst)
  {
    double z = BLANKf;
    double y = var.lscore(z, opt, cdatai);
    cdata.assign(0, y);
    if (y == BLANKf)
    { cdata.assign(1, BLANKs);  cdata.assign(2, BLANKs); }
    else
    { cdata.assign(1, z);  cdata.assign(2, 1.0 / (1.0 + exp(-z))); }
    cdata.put(fout);
  }
}

void logregclass::plot_logngp(option& opt, double cin, double pin, const varclass& var)
// Ngraph用データ出力
{
  ofstream fout;
  if (ofopen(fout, (opt.ngpname()).c_str()) == false) return; // ファイルオープン失敗
  bool delml = false;
  if (opt.ngp("log"))
  {
    fout << "-s1 -dline -w80 -c255,0,0 -zy-0.1,1.1,0.2\n";
    plot_ngpid(fout, opt, cin, pin, var);
    plot_log(fout, opt, var);
    delml = true;
  }
  if (opt.ngp("lci"))
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline -w80 -c255,0,0 -zy-0.1,1.1,0.2\n";
      delml = true;
    }
    plot_ngpid(fout, opt, cin, pin, var);
    plot_lci(fout, opt, var);
  }
  if (opt.ngp("score"))
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dmark,00 -o300 -zy-0.1,1.1,0.2\n";
      delml = true;
    }
    plot_ngpid(fout, opt, cin, pin, var);
    plot_score(fout, opt, var);
  }
  if (opt.ngp("or"))
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dmark,00 -w80 -o400 -zx0,"<< id.ni() << ",1 -zy0.1,100,10 -eylog\n";
      fout << "%CIを描く時は-s1 -derrorbar_y -w80\n";
    }
    plot_ngpid(fout, opt, cin, pin, var);
    plot_or(fout, opt, var);
  }
}

void logregclass::plot_ngpid(ostream& fout, option& opt, double cin, double pin, const varclass& var)
// Ngraph用データのタイトルとデータ識別情報を出力する
{
  fout << "%";  puttitle(fout, Title);
  fout << "%\n";
  id.put_title(fout, formngp);
  id.put_stra(fout, formngp);
  idd tid(id);
  tid.set(id.title(), id.ni() + (opt.ixlst()).size()); // 項目数を変更する
  set_itemix(tid, opt, id.ni());
  fout << "%目的変数y   :" << tid.item(0) << "\n";
  for (size_t i = 1; i < tid.ni(); ++i)
  {
    if (var.check(i))
    {
      fout << "%説明変数x" << setw(3) << i << ":" << tid.item(i);
      if (opt.comp()) fout << " 欠損補完用平均値=" << var.mid(i);
      fout << "\n";
    }
  }
  if (pin < 1.0) fout << "%取り込み基準:p値＜Pin=" << pin << "\n";
  else fout << "%取り込み基準:Waldのχ2乗値＞Cin=" << cin << "\n";
  if (opt.xno() > 0 && var.check(opt.xno())) // X軸にする説明変数指定あり
  {
    fout << "%X軸にする説明変数:" << tid.item(opt.xno()) << "\n";
    for (size_t i = 1; i < tid.ni(); ++i)
    {
      if (var.check(i))
      {
        if (opt.x(i) != BLANKf) fout << "%X軸にしない説明変数:" << tid.item(i) << " 指定値=" << opt.x(i) << "\n";
      }
    }
  }
}

void logregclass::plot_log(ostream& fout, option& opt, const varclass& var)
// Ngraph用ロジスティック曲線を出力する
{
  fout << "%ロジスティック曲線\n";
  datac cdata(id.ni());
  cdata.assign(0, var.b(0));    // 切片
  for (size_t i = 1; i < id.ni(); ++i)
  {
    if (opt.x(i) == BLANKf) cdata.assign(i, var.mid(i));
    else cdata.assign(i, opt.x(i));
  }
  double xmin = 0.0, xmax = 0.0;
  double xrange = cal_range(xmin, xmax, opt, cdata, var);
  double dx = xrange / (double)RDdef, z = 0.0, pr = 0.0;
  double x = xmin;
  for (int i = 0; i <= RDdef; ++i, x += dx)
  {
    z = 0.0;
    if (opt.xno() > 0 && var.check(opt.xno()))
    {
      cdata.assign(opt.xno(), x);
      var.lscore(z, opt, cdata);
    }
    else
    { z = x; }
    pr = 1.0 / (1.0 + exp(-z));
    fout << x << " " << pr << "\n";
  }
}

double logregclass::cal_range(double& xmin, double& xmax, option& opt, datac& cdata, const varclass& var)
// ロジスティック曲線のxの最小値と最大値を計算する
{
  if (opt.xmin() != BLANKf && opt.xmax() != BLANKf) // X軸の範囲指定あり
  {
    xmin = opt.xmin();  xmax = opt.xmax();
    return xmax - xmin;         // 範囲を返す
  }
  double pr = 0.001;            // 最小確率
  double zmin = log(pr / (1 - pr));
  pr = 1.0 - pr;                // 最大確率
  double zmax = log(pr / (1 - pr));
  if (opt.xno() > 0 && var.check(opt.xno())) // 変化説明変数が指定されている時
  {
    cdata.assign(opt.xno(), 0.0);
    double z = 0.0;
    var.lscore(z, opt, cdata);
    xmin = (zmin - z) / var.b(opt.xno());
    xmax = (zmax - z) / var.b(opt.xno());
    if (var.b(opt.xno()) < 0.0) // 偏回帰係数が負の時
    {
      z = xmin;
      xmin = xmax;
      xmax = z;
    }
  }
  else
  { xmin = zmin;  xmax = zmax; }
  return xmax - xmin;           // 範囲を返す
}

void logregclass::plot_score(ostream& fout, option& opt, const varclass& var)
// Ngraph用ロジットスコアを出力する
{
  fout << "%ロジットスコア\n";
  for (auto& cdatai : datalst)
  {
    double x = BLANKf, z = BLANKf;
    double y = var.lscore(z, opt, cdatai);
    if (opt.xno() < 1) x = z;
    else x = cdatai.fdata(opt.xno());
    if (x != BLANKf && y != BLANKf) fout << x << " " << y << "\n";
  }
}

void logregclass::plot_lci(ostream& fout, option& opt, const varclass& var)
// Ngraph用ロジスティック曲線の信頼区間を出力する
{
  datac cdata(id.ni());
  cdata.assign(0, var.b(0));    // 切片
  for (size_t i = 1; i < id.ni(); ++i)
  {
    if (opt.x(i) == BLANKf) cdata.assign(i, var.mid(i));
    else cdata.assign(i, opt.x(i));
  }
  double xmin = 0.0, xmax = 0.0;
  double xrange = cal_range(xmin, xmax, opt, cdata, var);
  double dx = xrange / (double)RDdef;
  double t = ptoz(1.0 - opt.cc()), z = 0.0, pr = 0.0, ci = 0.0;
  double x = xmin;
  fout << "%ロジスティック曲線の" << opt.cc() * 100.0 << "%信頼区間下限\n";
  for (int i = 0; i <= RDdef; ++i, x += dx)
  {
    if (opt.xno() > 0 && var.check(opt.xno()))
    {
      cdata.assign(opt.xno(), x);
      var.lscore(z, opt, cdata);
    }
    else
    { z = x; }
    pr = 1.0 / (1.0 + exp(-z));
    double prl = pr - t * sqrt(pr * (1.0 - pr) / (double)(var.n()));
    if (sgn(prl) <= 0) prl = 0.0;
    fout << x << " " << prl << "\n";
  }
  fout << NGPdelm2 << "\n";
  fout << "%ロジスティック曲線の" << opt.cc() * 100.0 << "%信頼区間上限\n";
  x = xmin;
  for (int i = 0; i <= RDdef; ++i, x += dx)
  {
    if (opt.xno() > 0 && var.check(opt.xno()))
    {
      cdata.assign(opt.xno(), x);
      var.lscore(z, opt, cdata);
    }
    else
    { z = x; }
    pr = 1.0 / (1.0 + exp(-z));
    double pru = pr + t * sqrt(pr * (1.0 - pr) / (double)(var.n()));
    if (pru > 1.0) pru = 1.0;
    fout << x << " " << pru << "\n";
  }
}

void logregclass::plot_or(ostream& fout, option& opt, const varclass& var) const
// Ngraph用にORとその信頼区間を出力する
{
  fout << "%オッズ比と" << opt.cc() * 100.0 << "%信頼区間\n";
  double z = ptoz(1.0 - opt.cc());
  for (int i = 1; i < var.np(); ++i)
  {
    if (var.check(i))
    {
      double ci = z * sqrt(var.vb(i));
      double ori = exp(var.b(i));
      double orl = exp(var.b(i) - ci), oru = exp(var.b(i) + ci);
      fout << i << " " << ori << " " << oru - ori << " " << orl - ori << "\n";
    }
  }
}

int logregclass::get(const char* c_fname, option& opt)
// データを読み込み保存・累積する
{
  if (c_fname == nullptr) fname.erase();
  else fname = c_fname;
  int nc = fgetdata(fname, id, datalst);
  if (nc < 0) throw runtime_error("実行時エラー！");
  if (opt.check(id.ni()) != 0) throw option_error("実行時オプションエラー！");
  return nc;
}

void logregclass::put(option& opt) // 集計結果を出力する
{
  varclass var(id.ni(), opt);   // 変数オブジェクト
  var.cal(opt, datalst);        // 各種の計算実行
  if (opt.res("res"))
  {
    put_resid(cout, opt, var);
    if (opt.res("stat")) put_stat(cout, opt, var);
    if (var.resvar() == reserr)
    { cout << "※ロジスティック関数計算不可能！\n\n" << endl; }
  }
  if (var.resvar() == reserr) return;
  bool putdns = false;
  if (opt.cin() > 0.0 || (0.0 < opt.pin() && opt.pin() < 1.0)) // 変数選択
  {
    var.stepsel(opt, datalst);
    if (opt.res("res") && opt.res("selx")) put_logreg(cout, opt, opt.cin(), opt.pin(), var);
    if (var.resvar() != reserr)
    {
      if (opt.dnsname() != "") put_dns(opt, var);
      if (opt.ngpname() != "") plot_logngp(opt, opt.cin(), opt.pin(), var);
      putdns = true;
      if (var.dfb() + 1 >= var.np() && opt.res("allx") == false) // 全変数を選択した時
      {
        if (opt.res("res")) cout << "\n" << endl;
        return;
      }
    }
  }
  var.allsel(opt, datalst);
  if (opt.res("res"))
  {
    if (opt.res("allx") || opt.res("def")) put_logreg(cout, opt, 0.0, 1.0, var);
    cout << "\n" << endl;
  }
  if (putdns || var.resvar() == reserr) return;
  if (opt.dnsname() != "") put_dns(opt, var);
  if (opt.ngpname() != "") plot_logngp(opt, 0.0, 1.0, var);
}

int main(int argc, char* argv[]) // メインルーチン
{
  logregclass logreg;           // メインクラス変数
  option opt(argv[0]);          // オプション
  try
  {
    opt.setdef();               // オプションの初期設定
    bool fok = true;            // ファイル処理フラッグ
    for (int i = 1; i <= argc; ++i)
    {
      if (i >= argc)
      {
        if (fok) logreg.get(nullptr, opt); // 標準入力から読み込む
      }
      else if (*(argv[i]) == CHRopt)
      { opt.set(argv[i] + 1);  fok = true; }
      else
      { logreg.get(argv[i], opt);  fok = false; }
    }
    logreg.put(opt);
  }
  catch (help_except& helpex)   // ヘルプファイル表示
  { return 0; }
  catch (option_error& opterr)  // オプション指定エラー
  { return 1; }
  catch (runtime_error& rterr)  // 実行時エラー
  {
    disperr(opt.pname(), rterr.what());
    return 2;
  }
  catch (...)                   // 予期せぬ例外
  {
    disperr(opt.pname(), "内容不明なエラー発生！");
    return 255;
  }
  return 0;
}
