// poireg.C Version 7.45  Copyleft (c) 杉本典夫
// Time-stamp: <2026-04-29 09:40:38 tac>
//
// Poisson回帰分析。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法  :g++ -O -s -W poireg poireg.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 =            // プログラム使用法表示用
{
  "poireg  Version 7.45         Copyleft(Vo\\oV) 2026 杉本典夫",
  "書式:poireg [<オプション> …]  [<データファイル> …]",
  "機能:最初の項目を目的変数Y、以後の項目を説明変数Xとした",
  "     Poisson回帰分析を行います。",
  "<オプション>",
  "  -b=<β0の初期値>[,<β1の初期値>,…]:偏回帰係数βの初期値を指定します。",
  "  -cc[=<信頼係数>]:信頼係数を指定します。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -eps[=<収束誤差>]:最尤法の収束誤差を指定します。(0<収束誤差)",
  "  -f=<ファイル>:指定したファイルからオプションを読み込みます。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示します。(-h:-h=man)",
  "  -ix=<説明変数の開始番号>[-<終了番号>]x<説明変数の開始番号>[-<終了番号>][,…]",
  "    交互作用項目を作成する2つの説明変数の番号を指定します。",
  "  -ixc[={-|+}]:交互作用項目の説明変数を中心化するかどうか指定します。",
  "    -:中心化しない(オプション省略時)  +:中心化する  -ixc:ixc=+",
  "  -meth[={c|x}]:計算方法を指定します。(-meth:-meth=x)",
  "  -model[={a|l|p}]:モデルを指定します。(-model:-model=a)",
  "  -ngp={ex|lambda|rr}[,…]:Ngraph用データを指定します。",
  "    ex:指数関数  lambda:ln(λ)とy  rr:リスク比",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データを出力するファイル名を指定します。",
  "    オプション省略時:Ngraph用データ出力せず  -ngpname:-ngpname=poireg_ngp.dat",
  "  -res=[-][{def|all|reg|stat|v}][,…]:解析結果の出力項目を指定します。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    reg:回帰分析結果  stat:基礎統計量  v:偏回帰係数初期値",
  "  -rmax=<最大回数>:反復計算の最大回数を指定します。(0≦最大回数)",
  "  -x<説明変数番号>=<説明変数の値>:X軸にしない説明変数の値を指定します。",
  "  -xno=<説明変数の番号>:X軸にする説明変数の番号を指定します。",
  "  -xrange=<最小値>-<最小値>:X軸の最小値と最大値を指定します。"
};

const string Title = "=== Poisson回帰分析(Poisson regression analysis) ==="; // 表題
const string NGPNAMEdef = "poireg_ngp.dat"; // 省略時Ngraphファイル名
constexpr double EPSdef = 1.0e-5; // 収束誤差省略時値
constexpr int RMAXdef = 200;      // 反復計算の省略時最大回数

enum methtype                   // 計算方法
{
  methnon,                      // 指定なし
  methc,                        // ケース単位
  methx,                        // 説明変数単位
  methetc                       // その他
};

enum modeltype                  // モデル
{
  modelnon,                     // 指定なし
  modela,                       // 最適モデルを自動選択
  modell,                       // λ(出現件数)モデル
  modelp,                       // p=λ/n(出現率)モデル
  modeletc                      // その他
};

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  lstr blst;                    // 説明変数の初期値リスト
  vecdbl v_bvec;                // 説明変数の初期値ベクトル
  double v_cc;                  // 信頼係数
  double v_eps;                 // 最尤法の収束誤差
  lstr ixstrlst;                // 交互作用項目の説明変数リスト
  lino v_ixlst;                 // 交互作用項目の説明変数ペアリスト
  bool v_ixc;                   // 交互作用項目の中心化フラッグ
  methtype v_meth;              // 計算方法
  modeltype v_model;            // モデル
  string v_ngpname;             // Ngraphファイル名
  lstr ngplst;                  // Ngraph出力リスト
  mapsb ngpmap;                 // Ngraph出力の有無マップ
  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 eps() const            // 最尤法の収束誤差を返す
  { return v_eps; }
  const lino& ixlst()    // 交互作用項目の説明変数番号ペアリストを返す
  { return v_ixlst; }
  bool ixc()                    // 交互作用項目の中心化フラッグを返す
  { return v_ixc; }
  methtype meth() const         // 計算方法を返す
  { return v_meth; }
  modeltype model() const       // モデルを返す
  { return v_model; }
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  bool ngp(const string& code); // Ngraph出力フラッグを返す
  bool res(const string& code); // 分析結果出力フラッグを返す
  int rmax() const              // 反復計算の最大回数を返す
  { return v_rmax; }
  double x(size_t no) const     // 説明変数の値を返す
  { return (0 < no && no <= xvec.size()) ? xvec[no - 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_eps(EPSdef), ixstrlst(), v_ixlst(), v_ixc(false), v_meth(methx), v_model(modela), v_ngpname(), ngplst(), ngpmap(), 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); // プログラム名削除
}

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] == "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] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  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] == "meth")    // 計算方法指定
  {
    if (opt[1].empty() || opt[1] == "x") v_meth = methx;
    else if (opt[1] == "c") v_meth = methc;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "model")   // モデル指定
  {
    if (opt[1].empty() || opt[1] == "a") v_model = modela;
    else if (opt[1] == "l") v_model = modell;
    else if (opt[1] == "p") v_model = modelp;
    else err = badopt(pname(), optstr, USAGE);
  }
  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] == "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(), "項目が2つ以上必要です！");
  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["ex"] = ngpmap["lambda"] = ngpmap["rr"] = true;
  if (ngplst.size() > 0)
  {
    if (v_ngpname == "") return disperr(pname(), "-ngpnameオプションの指定が必要！");
    ngpmap["ex"] = ngpmap["lambda"] = ngpmap["rr"] = 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["reg"] = resmap["stat"] = true;
  resmap["v"] = false;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["reg"] = 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["reg"] = resmap["stat"] = true; }
      else if (resstr == "all") // 全て出力
      {
        resmap["def"] = resmap["reg"] = 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        // ケースデータの説明変数比較用関数オブジェクト
{
public:
  data_less()                   // コンストラクタ
  {}
  ~data_less()                  // デストラクタ
  {}
  bool operator()(datac& cdata1, datac& cdata2) const; // ケースデータを比較する
};

bool data_less::operator()(datac& cdata1, datac& cdata2) const
// ケースデータの説明変数を比較する
{
  size_t ni = cdata1.size();
  int comp = 0;
  for (size_t i = 1; comp == 0 && i < ni; ++i)
  {
    double fdata1 = cdata1.fdata(i), fdata2 = cdata2.fdata(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 v_mevec;               // 平均値ベクトル
  vecdbl v_vrvec;               // 分散ベクトル
  vec<double> b0v;              // 偏回帰係数初期値ベクトル
  vec<double> bv;               // 偏回帰係数ベクトル
  vec<double> vbv;              // 偏回帰係数の分散ベクトル
  double v_lb;                  // ポアソンモデルの対数尤度
  double v_l0;                  // 切片モデルの対数尤度
  double v_lf;                  // 飽和モデルの対数尤度
  double v_st;                  // 全体平方和
  double v_sb;                  // 回帰平方和
  double v_slof;                // ズレ平方和
  int v_nx;                     // 説明変数の数
  int v_nxix;                   // 交互作用項目を含んだ説明変数の数
  int v_nr;                     // 合計個数
  int v_n0;                     // 0件の例数
  int v_n;                      // 例数
  int v_dft;                    // 全体自由度
  modeltype v_modelt;           // モデルの種類
  int v_nrep;                   // 反復計算回数
  restype v_resvar;             // 処理結果の種類
  int set_mean(int& nr, int& n0, vecdbl& mxvec, ldata& datalst) const;
  // 説明変数の平均値ベクトルをセットする
  void set_initb(vec<double>& b0v, vec<double>& vbv, option& opt) const;
  // 重回帰分析→判別分析による偏回帰係数初期値のセット
  int cal_poilf(int& nt, double& lf, const vecdbl& mevec, option& opt, ldata& datalst) const;
  // 飽和モデルのグループ数と例数と尤度を求める
  restype newton_poi(int& nrep, vec<double>& bv, vec<double>& vbv, double& lb, modeltype modelt, const vecdbl& mevec, int rmax, option& opt, ldata& datalst) const;
  // Newton法によるポアソン回帰の計算(最尤法)
  restype cal_lgh(vec<double>& gv, mat<double>& hmat, double& lb, const vec<double>& bv, modeltype modelt, const vecdbl& mevec, option& opt, ldata& datalst) const;
  // 傾斜ベクトルgとHessian行列Hを計算する
  bool data_equal(datac& cdata1, datac& cdata2) 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; }
  double me(int idx) const      // 平均値を返す
  { return (0 <= idx && idx < np()) ? v_mevec[idx] : 0.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; }
  double l0() const             // 切片モデルの対数尤度を返す
  { return v_l0; }
  double lf() const             // 飽和モデルの対数尤度を返す
  { return v_lf; }
  double st() const             // 全体平方和を返す
  { return v_st; }
  double sb() const             // 回帰平方和を返す
  { return v_sb; }
  double slof() const           // ズレ平方和を返す
  { return v_slof; }
  double r2() const             // 疑似重寄与率を返す
  { return (st() > 0.0) ? sb() / st() : 0.0; }
  int nx() const                // 説明変数の数を返す
  { return v_nx; }
  int nxix() const              // 交互作用項目を含んだ説明変数の数を返す
  { return v_nxix; }
  int np() const                // パラメータ数を返す
  { return nxix() + 1; }
  int nr() const                // 合計個数を返す
  { return v_nr; }
  int n0() const                // 0件の例数を返す
  { return v_n0; }
  int n() const                 // 例数を返す
  { return v_n; }
  int dft() const               // 全体自由度を返す
  { return v_dft; }
  int dfb() const               // 回帰自由度を返す
  { return np() - 1; }
  int dflof() const             // ズレ自由度を返す
  { return dft() - dfb(); }
  modeltype modelt() const      // モデルの種類を返す
  { return v_modelt; }
  int nrep() const              // 反復計算回数を返す
  { return v_nrep; }
  restype resvar() const        // 処理結果の種類を返す
  { return v_resvar; }
  void cal(option& opt, ldata& datalst); // データを集計して各種統計量を求める
  void poireg(option& opt, ldata& datalst); // ポアソン回帰分析実行
  double lambda(option& opt, datac& cdata) const; // λまたはλ/n推測値を返す
};

varclass::varclass(size_t ni, option& opt) : v_smat(), v_mevec(), v_vrvec(), b0v(), bv(), vbv(), v_lb(0.0), v_l0(0.0), v_lf(0.0), v_st(0.0), v_sb(0.0), v_slof(0.0), v_nx(0), v_nxix(0), v_nr(0), v_n0(0), v_n(0), v_dft(0), v_modelt(modela), v_nrep(0), v_resvar(resnon)
// コンストラクタ
{
  v_nx = ni - 1;
  v_nxix = v_nx + (opt.ixlst()).size();
  int dim = v_nxix + 2;
  v_smat.resize(dim, dim);
  v_mevec.resize(dim, 0.0);  v_vrvec.resize(dim, 0.0);
  b0v.resize(dim - 1);  bv.resize(dim - 1);  vbv.resize(dim - 1);
}

void varclass::cal(option& opt, ldata& datalst)
// データを集計して各種統計量を求める
{
  vecdbl mxvec(nx() + 1, 0.0);
  v_n = set_mean(v_nr, v_n0, mxvec, datalst); // 説明変数の平均値ベクトルをセットする
  if (v_n < np() || v_nr < 1)   // ポアソン回帰式計算不可能
  {
    v_resvar = reserr;
    return;
  }
  int dim = np() + 1;
  vecdbl fdatavec(dim, 0.0);
  for (auto& cdatai : datalst)
  {
    if (cdatai.fdata(0) != BLANKf) // 有効なケース
    {
      setixdata(fdatavec, opt.ixlst(), opt.ixc(), mxvec, 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]);
      }
    }
  }
  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);
  for (int i = 0; i < dim; ++i)
  {
    v_mevec[i] = v_smat(0, i);
    v_vrvec[i] = (v_n > 1) ? v_smat(i, i) / (double)(v_n - 1) : 0.0;
  }
  v_resvar = resok;
}

int varclass::set_mean(int& nr, int& n0, vecdbl& mxvec, ldata& datalst) const
// 説明変数の平均値ベクトルをセットする
{
  size_t nmx = mxvec.size();
  vecdbl fdatavec(nmx, 0.0);
  int n = 0;
  nr = 0;  n0 = 0;
  for (auto& cdatai : datalst)
  {
    if (cdatai.fdata(0) != BLANKf) // 目的変数有り
    {
      bool blank = false;
      for (size_t j = 1; blank == false && j < nmx; ++j)
      {
        fdatavec[j] = cdatai.fdata(j);
        if (fdatavec[j] == BLANKf) blank = true;
      }
      if (blank)
      { cdatai.assign(0, BLANKf); } // 後の集計のために目的変数を欠損にする
      else
      {
        ++n;
        fdatavec[0] = cdatai.fdata(0);
        if (fdatavec[0] <= 0.0) ++n0; // 目的変数が0の例数
        nr += (int)(fdatavec[0]); // 合計個数を数える
        for (size_t j = 0; j < nmx; ++j) mxvec[j] += fdatavec[j];
      }
    }
  }
  for (size_t j = 0; j < nmx; ++j) mxvec[j] /= (double)n;
  return n;                     // 例数を返す
}

void varclass::poireg(option& opt, ldata& datalst)
// ポアソン回帰分析実行
{
  if (n0() < 1) v_modelt = modell; // 発生件数0がない時は発生件数モデル
  else if (opt.model() == modela) v_modelt = modelp;
  else v_modelt = opt.model();
  set_initb(b0v, vbv, opt);
  bv = b0v;
  v_resvar = newton_poi(v_nrep, bv, vbv, v_lb, v_modelt, v_mevec, opt.rmax(), opt, datalst);
  if (v_resvar == reserr)    // Newton法が発散した時は初期値を解にする
  {
    bv = b0v;
    v_resvar = newton_poi(v_nrep, bv, vbv, v_lb, v_modelt, v_mevec, 1, opt, datalst);
  }
  int nt = 0;
  int m = cal_poilf(nt, v_lf, v_mevec, opt, datalst);
  v_dft = m - 1;
  vec<double> gv(np()), bv(np());
  mat<double> hmat(np());
  double lambda0 = (v_modelt == modelp) ? (double)(nr()) / (double)nt : (double)(nr()) / (double)m;
  bv.assign(0, log(lambda0));
  cal_lgh(gv, hmat, v_l0, bv, v_modelt, v_mevec, opt, datalst);
  v_st = v_lf - v_l0;
  if (v_dft < 1 || v_st < FDeps) v_st = 0.0;
  v_sb = v_lb - v_l0;
  if (dfb() < 1 || v_sb < FDeps) v_sb = 0.0;
  v_slof = v_st - v_sb;
  if (dflof() < 1 || v_slof < FDeps) v_slof = 0.0;
}

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

int varclass::cal_poilf(int& nt, double& lf, const vecdbl& mevec, option& opt, ldata& datalst) const
// 飽和モデルのグループ数と例数と尤度を求める
{
  int m = 0;
  nt = 0;  lf = 0.0;
  ldata_itr iitr = datalst.begin();
  while (iitr != datalst.end())
  {
    if (iitr->fdata(0) == BLANKf)
    { ++iitr; }
    else
    {
      int ri = (int)(iitr->fdata(0)); // 反応有の個数
      int ni = (ri > 0) ? ri : 1;   // 全体の例数
      ldata_itr jitr = iitr;
      ++jitr;
      if (opt.meth() != methc)  // 説明変数単位の計算方法
      {
        bool equal = true;
        while (equal && jitr != datalst.end())
        {
          if (jitr->fdata(0) == BLANKf)
          { ++jitr; }
          else
          {
            equal = data_equal(*iitr, *jitr);
            if (equal)
            {
              int rj = int(jitr->fdata(0));
              ri += rj;
              ni += (rj > 0) ? rj : 1;
              ++jitr;
            }
          }
        }
      }
      double ri_f = (double)ri, ni_f = (double)ni;
      double lambdai = ri_f;
      if (ri > 0)               // 発生件数がある時だけ尤度計算
      { lf += ri_f * log(lambdai) - lgamma(ri_f + 1.0) - lambdai; }
      nt += ni;  ++m;
      iitr = jitr;
    }
  }
  return m;                     // グループ数を返す
}

restype varclass::newton_poi(int& nrep, vec<double>& bv, vec<double>& vbv, double& lb, modeltype modelt, const vecdbl& mevec, int rmax, option& opt, ldata& datalst) const
// Newton法によるポアソン回帰の計算(最尤法)
{
  vec<double> gv(np());
  mat<double> hmat(np());
  double lold = 0.0;
  restype resnew = resnon;
  nrep = 0;
  datalst.sort(data_less());    // 説明変数順にソートする
  while (resnew == resnon && cal_lgh(gv, hmat, lb, bv, modelt, 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 < np(); ++i) vbv.assign(i, -hmat(i, i));
  return resnew;                // Newton法成功
}

restype varclass::cal_lgh(vec<double>& gv, mat<double>& hmat, double& lb, const vec<double>& bv, modeltype modelt, const vecdbl& mevec, option& opt, ldata& datalst) const
// 傾斜ベクトルgとHessian行列Hを計算する
{
  constexpr double Emax = 308.0; // 指数最大値
  vec<double> xv(np());
  xv.assign(0, 1.0);  gv = 0.0;  hmat = 0.0;  lb = 0.0;
  vecdbl fdatavec(np() + 1, 0.0);
  restype reslgh = resnon;
  ldata_itr iitr = datalst.begin();
  while (reslgh == resnon && iitr != datalst.end())
  {
    if (iitr->fdata(0) == BLANKf)
    { ++iitr; }
    else
    {
      int ri = (int)(iitr->fdata(0)); // 反応有の個数
      int ni = (ri > 0) ? ri : 1;   // 全体の例数
      ldata_itr jitr = iitr;
      ++jitr;
      if (opt.meth() != methc)  // 説明変数単位の計算方法
      {
        bool equal = true;
        while (equal && jitr != datalst.end())
        {
          if (jitr->fdata(0) == BLANKf)
          { ++jitr; }
          else
          {
            equal = data_equal(*iitr, *jitr);
            if (equal)
            {
              int rj = (int)(jitr->fdata(0));
              ri += rj;
              ni += (rj > 0) ? rj : 1;
              ++jitr;
            }
          }
        }
      }
      setixdata(fdatavec, opt.ixlst(), opt.ixc(), mevec, *iitr);
      for (int i = 1; i < np(); ++i) xv.assign(i, fdatavec[i]);
      double sbxi = inp(bv, xv);
      if (fabs(sbxi) > Emax)    // 指数最大値以上
      { reslgh = reserr; }
      else
      {
        double wi = exp(sbxi), ri_f = (double)ri, ni_f = (double)ni;
        if (modelt == modelp)   // 出現率モデル
        {
          double nwi = ni_f * wi;
          lb += ri_f * (log(ni_f) + sbxi) - lgamma(ri_f + 1.0) - nwi;
          for (int j = 0; j < np(); ++j)
          {
            gv.add(j, (ri_f - nwi) * xv(j));
            for (int k = j; k < np(); ++k) hmat.add(j, k, -nwi * xv(j) * xv(k));
          }
        }
        else                    // 出現例数モデル
        {
          lb += ri_f * sbxi - lgamma(ri_f + 1.0) - wi; // ri=0も含む
          for (int j = 0; j < np(); ++j)
          {
            gv.add(j, (ri_f - wi) * xv(j));
            for (int k = j; k < np(); ++k) hmat.add(j, k, -wi * xv(j) * xv(k));
          }
        }
      }
      iitr = jitr;
    }
  }
  if (reslgh == reserr) return reslgh;
  for (int i = 0; i < np() - 1; ++i)
  {
    for (int j = i + 1; j < np(); ++j) hmat.assign(j, i, hmat(i, j));
  }
  return resok;
}

bool varclass::data_equal(datac& cdata1, datac& cdata2) const
// ケースデータの説明変数を比較する
{
  size_t nx = cdata1.size() - 1;
  bool equal = true;
  for (size_t i = 1; equal && i <= nx; ++i)
  {
    double fdata1 = cdata1.fdata(i), fdata2 = cdata2.fdata(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:不一致
}

double varclass::lambda(option& opt, datac& cdata) const
// λまたはλ/n推測値を返す
{
  if (cdata.fdata(0) == BLANKf) return BLANKf; // 目的変数欠損のため計算できず
  vecdbl fdatavec(np() + 1, 0.0);
  setixdata(fdatavec, opt.ixlst(), opt.ixc(), v_mevec, cdata);
  double bx = b(0);
  for (int i = 1; i < np(); ++i) bx += b(i) * fdatavec[i];
  return exp(bx);               // λまたはλ/n推測値を返す
}

class poiregclass               // poiregメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout, option& opt);
  // 解析結果のタイトルとデータ識別情報を出力する
  int set_itemix(idd& tid, option& opt) const;
  // 交互作用項目名をセットする
  void put_poireg(ostream& fout, option& opt, const varclass& var) const;
  // Poisson回帰分析結果を出力する
  void put_brci(ostream& fout, option& opt, const varclass& var) const;
  // 偏回帰係数とリスク比の信頼区間を出力する
  void put_poitab(ostream& fout, option& opt, const varclass& var) const;
  // Poisson回帰分析の分散分析表を出力する
  void plot_poingp(option& opt, const varclass& var);
  // Ngraph用データ出力
  void plot_ngpid(ofstream& fout, option& opt);
  // Ngraph用データのタイトルとデータ識別情報を出力する
  void plot_ex(ofstream& fout, option& opt, const varclass& var);
  // Ngraph用に指数曲線を出力する
  double search_xrange(double& xmin, double& xmax, option& opt, const varclass& var);
  // x軸の最小値と最大値を検索して範囲を返す
  void plot_lambda(ofstream& fout, option& opt, const varclass& var);
  // Ngraph用にln(λ)とyを出力する
  void plot_rr(ofstream& fout, option& opt, const varclass& var) const;
  // Ngraph用にリスク比とその信頼区間を出力する
public:
  poiregclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~poiregclass()                // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み累積する
  void put(option& opt);                     // 集計結果を出力する
};

void poiregclass::put_resid(ostream& fout, option& opt)
// 解析結果のタイトルとデータ識別情報を出力する
{
  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);
  fout << "目的変数y   :" << tid.item(0) << "\n";
  for (size_t j = 1; j < tid.ni(); ++j)
  {
    fout << "説明変数x" << setw(3) << j;
    fout << ":" << tid.item(j) << "\n";
  }
}

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

void poiregclass::put_poireg(ostream& fout, option& opt, const varclass& var) const
// Poisson回帰分析結果を出力する
{
  constexpr int LENline = 78;   // 区切り線の長さ
  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"; }
  fout << "  計算方法:";
  if (opt.meth() == methc) fout << "ケース単位(SASと同じ)\n";
  else fout << "説明変数単位\n";
  fout << "  モデル:";
  if (var.modelt() == modelp) fout << "p=λ/n=exp(β0+Σβj・xj)  n:例数";
  else fout << "λ=exp(β0+Σβj・xj)";
  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 << "変数  偏回帰係数   標準誤差     exp(bj)  偏回帰係数   Waldのχ^2       p値\n";
  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)
  {
    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);
  put_poitab(fout, opt, var);
}

void poiregclass::put_brci(ostream& fout, option& opt, const varclass& var) const
// 偏回帰係数とリスク比の信頼区間を出力する
{
  fout << "\n変数  偏回帰係数" << setw(5) << opt.cc() * 100.0;
  fout << "%CI下限        上限     exp(bj)" << 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)
  {
    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 poiregclass::put_poitab(ostream& fout, option& opt, const varclass& var) const
// Poisson回帰分析の分散分析表を出力する
{
  fout << "対数尤度:L(モデル)=" << var.lb() << "  L(切片)=" << var.l0() << "  飽和Lf=" << var.lf() << "\n";
  fout << "擬似寄与率R2乗=" << var.r2() << "  AIC(赤池の情報量基準)=" << 2.0 * (var.np() - 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 * var.sb();
  ldbl pr = ctop(csqr, var.dfb());
  fout << "回帰" << setw(22) << var.sb() << setw(10) << var.dfb() << setw(14) << csqr;
  fout << setw(14) << pr << pind(pr) << "\n";
  csqr = (var.dflof() < 1) ? 0.0 : 2.0 * var.slof();
  pr = ctop(csqr, var.dflof());
  fout << "ズレ(LOF)" << setw(17) << var.slof() << setw(10) << var.dflof();
  fout << setw(14) << csqr << setw(14) << pr << pind(pr) << "\n";
  putline(fout, LENline);
  fout << "全体" << setw(22) << var.st() << setw(10) << var.dft() << endl;
}

void poiregclass::plot_poingp(option& opt, const varclass& var)
// Ngraph用データ出力
{
  ofstream fout;
  if (ofopen(fout, (opt.ngpname()).c_str()) == false) return; // ファイルオープン失敗
  bool delml = false;
  if (opt.ngp("ex"))
  {
    fout << "-s1 -dline -w80 -c255,0,0\n";
    plot_ngpid(fout, opt);
    plot_ex(fout, opt, var);
    delml = true;
  }
  if (opt.ngp("lambda"))
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dmark,00 -o300\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_lambda(fout, opt, var);
  }
  if (opt.ngp("rr"))
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dmark,00 -w80 -o400 -zx0,"<< id.ni() << ",1 -zy0.1,100,10 -eypoi\n";
      fout << "%CIを描く時は-s1 -derrorbar_y -w80\n";
    }
    plot_ngpid(fout, opt);
    plot_rr(fout, opt, var);
  }
}

void poiregclass::plot_ngpid(ofstream& fout, option& opt)
// 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);
  fout << "%目的変数y   :" << tid.item(0) << "\n";
  for (size_t i = 1; i < tid.ni(); ++i)
  { fout << "%説明変数x" << setw(3) << i << ":" << tid.item(i) << "\n"; }
  if (opt.xno() > 0)            // X軸にする説明変数指定あり
  {
    fout << "%X軸にする説明変数:" << tid.item(opt.xno()) << "\n";
    for (size_t i = 1; i < tid.ni(); ++i)
    {
      if (opt.x(i) != BLANKf) fout << "%X軸にしない説明変数:" << tid.item(i) << " 指定値=" << opt.x(i) << "\n";
    }
  }
}

void poiregclass::plot_ex(ofstream& fout, option& opt, const varclass& var)
// Ngraph用指数曲線を出力する
{
  fout << "%指数曲線:";
  if (var.modelt() == modelp) fout << "p=λ/n=exp(b0+Σbj・xj)";
  else fout << "λ=exp(b0+Σbj・xj)";
  fout << "　計算方法:";
  if (opt.meth() == methc) fout << "ケース単位\n";
  else fout << "説明変数単位\n";
  constexpr int RDdef = 50;     // 曲線の分解数
  double xmin = 0.0, xmax = 0.0;
  double xrange = search_xrange(xmin, xmax, opt, var);
  double dx = xrange / (double)RDdef;
  datac cdata(id.ni());
  cdata.assign(0, 0.0);
  for (size_t i = 1; i < id.ni(); ++i)
  {
    if (opt.x(i) == BLANKf) cdata.assign(i, var.me(i));
    else cdata.assign(i, opt.x(i));
  }
  double x = xmin;
  for (int i = 0; i <= RDdef; ++i, x += dx)
  {
    double lambda = BLANKf;
    if (opt.xno() < 1)
    { lambda = exp(x); }
    else
    {
      cdata.assign(opt.xno(), x);
      lambda = var.lambda(opt, cdata);
    }
    if (lambda != BLANKf) fout << x << " " << lambda << "\n";
  }
}

double poiregclass::search_xrange(double& xmin, double& xmax, option& opt, const varclass& var)
// x軸の最小値と最大値を検索して範囲を返す
{
  if (opt.xmin() != BLANKf && opt.xmax() != BLANKf) // X軸の範囲指定あり
  {
    xmin = opt.xmin();  xmax = opt.xmax();
    return xmax - xmin;         // 範囲を返す
  }
  xmin = FDmax;  xmax = -FDmax;
  for (auto& cdatai : datalst)
  {
    double x = BLANKf;
    if (opt.xno() > 0)
    { x = cdatai.fdata(opt.xno()); }
    else
    {
      double lambda = var.lambda(opt, cdatai);
      x = (lambda == BLANKf) ? BLANKf : log(lambda);
    }
    if (x != BLANKf)
    {
      if (x < xmin) xmin = x;
      if (x > xmax) xmax = x;
    }
  }
  double xrange = xmax - xmin, xmg = 0.1;
  xmin -= xrange * xmg;  xmax += xrange * xmg;
  return xmax - xmin;           // 範囲を返す
}

void poiregclass::plot_lambda(ofstream& fout, option& opt, const varclass& var)
// Ngraph用にln(λ)とyを出力する
{
  if (var.modelt() == modelp) fout << "%ln(λ/n)とy";
  else fout << "%ln(λ)とy";
  fout << "　計算方法:";
  if (opt.meth() == methc) fout << "ケース単位\n";
  else fout << "説明変数単位\n";
  for (auto& cdatai : datalst)
  {
    double x = BLANKf, y = cdatai.fdata(0);
    if (opt.xno() > 0)
    { x = cdatai.fdata(opt.xno()); }
    else
    {
      double lambda = var.lambda(opt, cdatai);
      x = (lambda == BLANKf) ? BLANKf : log(lambda);
    }
    if (x != BLANKf && y != BLANKf) fout << x << " " << y << "\n";
  }
}

void poiregclass::plot_rr(ofstream& fout, option& opt, const varclass& var) const
// Ngraph用にリスク比とその信頼区間を出力する
{
  fout << "%リスク比と信頼区間　モデル:";
  if (var.modelt() == modelp) fout << "p=λ/n";
  else fout << "λ";
  fout << "　計算方法:";
  if (opt.meth() == methc) fout << "ケース単位\n";
  else fout << "説明変数単位\n";
  double z = ptoz(1.0 - opt.cc());
  for (int i = 1; i < var.np(); ++i)
  {
    double ci = z * sqrt(var.vb(i));
    double rri = exp(var.b(i));
    double rrl = exp(var.b(i) - ci), rru = exp(var.b(i) + ci);
    fout << i << " " << rri << " " << rru - rri << " " << rrl - rri << "\n";
  }
}

int poiregclass::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 poiregclass::put(option& opt) // 集計結果を出力する
{
  varclass var(id.ni(), opt);   // 変数オブジェクト
  var.cal(opt, datalst);        // 各種の計算実行
  if (opt.res("res"))
  {
    put_resid(cout, opt);
    if (opt.res("stat"))
    {
      putmsdvec(cout, var.n(), var.mevec(), var.vrvec(), var.nxix(), 1);
      putmat(cout, "・相関行列(correlation matrix)", matr, var.n(), var.smat(), var.nxix(), 1);
    }
    if (var.resvar() == reserr)
    { cout << "※ポアソン回帰式計算不可能！\n\n" << endl; }
  }
  if (var.resvar() == reserr) return;
  var.poireg(opt, datalst);
  if (opt.res("res"))
  {
    if (opt.res("reg")) put_poireg(cout, opt, var);
    cout << "\n" << endl;
  }
  if (var.resvar() != reserr && opt.ngpname() != "") plot_poingp(opt, var);
}

int main(int argc, char* argv[]) // メインルーチン
{
  poiregclass poireg;           // メインクラス変数
  option opt(argv[0]);          // オプション
  try
  {
    opt.setdef();               // オプションの初期設定
    bool fok = true;            // ファイル処理フラッグ
    for (int i = 1; i <= argc; ++i)
    {
      if (i >= argc)
      {
        if (fok) poireg.get(nullptr, opt); // 標準入力から読み込む
      }
      else if (*(argv[i]) == CHRopt)
      { opt.set(argv[i] + 1);  fok = true; }
      else
      { poireg.get(argv[i], opt);  fok = false; }
    }
    poireg.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;
}
