// dla.C Version 7.44  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-29 08:55:25 tac>
//
// 微分ロジット分析を行う。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -W dla dla.C
//               -I/usr/local/include -L/usr/local/lib -ldans -lfun

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

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

const vecstr USAGE =            // プログラム使用法表示用
{
  "dla  Version 7.44            Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:dla  [<オプション> …]  [<データファイル> …]",
  "機能:最初の項目を説明変数、次の項目を目的変数として微分ロジット分析を行う。",
  "<オプション>", 
  "  -cc[=<信頼係数>]:信頼係数を指定する。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -dnsname[=<出力ファイル名>]:DANSデータを出力するファイル名を指定する。",
  "    オプション省略時:DANSデータ出力せず  -dnsname:-dnsname=dla_dns.d",
  "    出力項目:説明変数,目的変数,目的変数推測値,目的変数ロジット,ロジット推測値",
  "  -f=<ファイル名>:指定されたファイルからオプションを読み込む。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示する。(-h:-h=man)",
  "  -ngp={logc|regc|regd|regl}[,…]:Ngraph用データの指定。",
  "    logc:xとyとロジスティック曲線  regc:yとdyと2次回帰曲線",
  "    regd:xとdyと微分ロジスティック曲線  regl:xとlyとロジット直線",
  "    オプション省略時:全て指定",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データを出力するファイル名を指定する。",
  "    オプション省略時:Ngraph用データ出力せず  -ngpname:-ngpname=dla_ngp.dat",
  "  -outd={app|comp|exc}:範囲外データの処理方法を指定する。",
  "    app:近似法  comp:反復補完法(オプション省略時)  exc:除外",
  "  -res=[-][{def|all|regc|regl|stat|y[{Mh|m|<yの値>}]}][,…]",
  "      :解析結果の出力項目を指定する。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    regc:yとdyの2次回帰分析結果  regl:ロジット回帰分析結果  stat:基礎統計量",
  "    y[{Mh|m|<yの値>}]:xの逆推定値",
  "  -xrange=<最小値>-<最小値>:X軸の最小値と最大値を指定する。"
};

const string Title = "=== 微分ロジット分析(differential logit analysis) ==="; // 表題
const string DNSNAMEdef = "dla_dns.d"; // 省略時DANSファイル名
const string NGPNAMEdef = "dla_ngp.dat"; // 省略時Ngraphファイル名
constexpr int RDdef = 50;       // グラフの描画範囲分割数省略時値
constexpr int Pdef = 2;         // 回帰式の次数
constexpr double EPSly = 0.001;   // 範囲外ロジットの計算用微小値
constexpr int RMAXdef = 200;      // 省略時反復計算最大回数

enum outdtype                   // 範囲外データ処理方法の種類
{
  outdnon,                      // なし
  outdapp,                      // 近似法
  outdcomp,                     // 繰り返し補完法
  outdexc,                      // 除外
  outdetc                       // その他
};

class option                    // オプション集合体
{
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  string v_dnsname;             // DANSファイル名
  string v_ngpname;             // Ngraphファイル名
  lstr ngplst;                  // Ngraph出力リスト
  mapsb ngpmap;                 // NgraphS出力の有無マップ
  double v_cc;                  // 信頼係数
  outdtype v_outd;              // 範囲外データの処理方法
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
  double v_y;                   // xを逆算するyの値
  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 string& dnsname() const // DANSデータファイル名を返す
  { return v_dnsname; }
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  bool ngp(const string& code); // Ngraph出力フラッグを返す
  double cc() const             // 信頼係数を返す
  { return v_cc; }
  outdtype outd() const         // 範囲外データの処理方法を返す
  { return v_outd; }
  bool res(const string& code); // 分析結果出力フラッグを返す
  double xmin() const           // X軸の最小値を返す
  { return xrange.first; }
  double xmax() const           // X軸の最大値を返す
  { return xrange.second; }
  double y() const              // xを逆算するyの値
  { return v_y; }
  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オプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_dnsname(), v_ngpname(), ngplst(), ngpmap(), v_cc(CCdef), v_outd(outdcomp), reslst(), resmap(), xrange(BLANKf, BLANKf), v_y(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 ngppos = ngpmap.find(code);
  if (ngppos == ngpmap.end()) return false; // code無し！
  return ngpmap[code];                      // true/falseを返す
}

bool option::res(const string& code) // 分析結果出力フラッグを返す
{
  mapsb_citr respos = resmap.find(code);
  if (respos == 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] == "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] == "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] == "f")       // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  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] == "outd")       // 範囲外データの処理方法指定
  {
    if (opt[1].empty() || opt[1] == "comp") v_outd = outdcomp;
    else if (opt[1] == "app") v_outd = outdapp;
    else if (opt[1] == "exc") v_outd = outdexc;
    else 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] == "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                          // 間違ったオプション指定
  { 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個！");
  if (set_ngp() < 0) return -2; // -ngpオプションのセット
  if (set_res() < 0) return -3; // -resオプションのセット
  return 0;                     // エラーなし
}

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

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["stat"] = resmap["regc"] = resmap["regl"] = resmap["y"] = true;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["stat"] = resmap["regc"] = resmap["regl"] = resmap["y"] = false;
    for (const string& resstr : reslst)
    {
      mapsb_citr respos = resmap.find(resstr);
      if (resstr == "-")        // 全て出力せず
      { resmap["res"] = false; }
      else if (resstr[0] == '-' && resstr.size() > 1) // 特定の項目を出力せず
      {
        respos = resmap.find(resstr.substr(1));
        if (respos == resmap.end())
        { return disperr(pname(), "-res=" + resstr + "オプション指定が不適当！"); }
        else
        { resmap[resstr.substr(1)] = false; }
      }
      else if (resstr == "def" || resstr == "all") // 省略時出力=全て出力
      {
        resmap["def"] = resmap["stat"] = resmap["regc"] = resmap["regl"] = resmap["y"] = true;
      }
      else if (resstr == "yMh")    // vy=BLANKf:M/2に置き換える
      { resmap["y"] = true; }
      else if (resstr[0] == 'y' && resstr.size() > 1)
      {
        if (resstr == "ym") v_y = DBL_MAX; // 平均値に置き換える
        else v_y = strtodf(resstr.substr(1));
        if (v_y <= 0.0) return disperr(pname(), "-res=" + resstr + "オプション指定が不適当！");
        resmap["y"] = true;
      }
      else if (respos == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適当！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

class data_less                 // ケースデータの比較用関数オブジェクト
{
private:
  int kidx;                     // ソートキー項目のインデックス
public:
  data_less(int kidx0) : kidx(kidx0) // コンストラクタ
  {}
  ~data_less()                  // デストラクタ
  {}
  bool operator()(datac& cdata1, datac& cdata2) const; // ケースデータを比較する
};

bool data_less::operator()(datac& cdata1, datac& cdata2) const
// ケースデータを比較する
{
  int comp = 0;
  double fdata1 = cdata1.fdata(kidx);
  double fdata2 = cdata2.fdata(kidx);
  if (fdata1 < fdata2) comp = -1;
  else if (fdata1 > fdata2) comp = 1;
  return (comp < 0) ? true : false; // 小さいときだけtrueを返す
}

class varclass                  // 変数クラス
{
private:
  pairdbl xrange;               // xの最小値と最大値
  pairdbl yrange;               // yの最小値と最大値
  statcor<double> statcxy;      // xとyの相関基礎統計量
  statcor<double> statcydy;     // yとdyの相関基礎統計量
  statcor<double> statcxly;     // xとlyの相関基礎統計量
  int v_np;                     // 微分回帰式の次数
  int v_noutd;                  // 範囲外データ数
  int v_nr;                     // 反復計算回数
  mat<double> smat;             // yとdyの微分回帰分析用積和行列
  double v_y00;                 // y0初期値
  double v_c0;                  // c初期値
  double v_m0;                  // M初期値
  double v_y0;                  // y0
  double v_c;                   // c
  double v_m;                   // M
  restype v_resregc;            // 微分回帰分析の結果フラッグ
  restype v_resvar;             // 処理結果のフラッグ
  int cal_smat(ldata& datalst, option& opt);
  // 微分回帰分析用に積和行列を計算する
  restype cal_y0(double& y0, double c, double m, ldata& datalst, option& opt);
  // y0をロジット分析によって求める
  restype cal_reglog(ldata& datalst, option& opt);
  // 微分ロジット分析によってパラメータを推定する(Mは更新しない)
  restype cal_replog(ldata& datalst, option& opt);
  // 反復微分ロジット分析によってパラメータを更新する
public:
  varclass() : xrange(FDmax, -FDmax), yrange(FDmax, -FDmax), statcxy(), statcydy(), statcxly(), v_np(Pdef), v_noutd(0), v_nr(0), smat(), v_y00(0.0), v_c0(0.0), v_m0(0.0), v_y0(0.0), v_c(0.0), v_m(0.0), v_resregc(resnon), v_resvar(resnon)
  // デフォールトコンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  void cal(ldata& datalst, option& opt);
  // 微分ロジット分析を実行する
  restype cal_mci(double& mlow, double& mup, double cc0) const;
  // Mの信頼区間を計算する(Newton法)
  double cal_dyest(double y) const; // dy推定値を返す
  double cal_dycr(double y, double x, bool plim) const; // dyの限界幅を返す
  int np() const                // 微分回帰式の次数を返す
  { return v_np; }
  int noutd() const             // 範囲外データ数を返す
  { return v_noutd; }
  int nr() const                // 反復計算数を返す
  { return v_nr; }
  double xmin() const           // xの最小値を返す
  { return xrange.first; }
  double xmax() const           // xの最大値を返す
  { return xrange.second; }
  double ymin() const           // yの最小値を返す
  { return yrange.first; }
  double ymax() const           // yの最大値を返す
  { return yrange.second; }
  int xyn() const               // xとyの例数を返す
  { return statcxy.n(); }
  double xymx() const           // xの平均値を返す
  { return statcxy.mx(); }
  double xymy() const           // yの平均値を返す
  { return statcxy.my(); }
  double xyvxx() const          // xの不偏分散を返す
  { return statcxy.vxx(); }
  double xyvyy() const          // yの不偏分散を返す
  { return statcxy.vyy(); }
  int ydyn() const              // yとdyの例数を返す
  { return statcydy.n(); }
  int ydydft() const            // yとdyの全体自由度を返す
  { return statcydy.n(); }
  double ydymy() const          // yの平均値を返す
  { return statcydy.mx(); }
  double ydymdy() const         // dyの平均値を返す
  { return statcydy.my(); }
  double ydyvyy() const         // yの不偏分散を返す
  { return statcydy.vxx(); }
  double ydyvdydy() const       // dyの不偏分散を返す
  { return statcydy.vyy(); }
  double ydysyy() const         // yの平方和を返す
  { return statcydy.sxx(); }
  double ydysdydy() const       // dyの平方和を返す
  { return statcydy.syy(); }
  double ydysr() const          // yとdyの微分回帰分析の残差を返す
  { return smat(np(), np()); }
  int ydydfr() const            // yとdyの残差自由度を返す
  { return statcydy.n() - np(); }
  double ydyvr() const          // yとdyの微分回帰分析の残差分散を返す
  { return (ydydfr() > 0) ? ydysr() / (double)ydydfr() : 0.0; }
  double ydyb(int idx) const // yとdyの偏回帰係数を返す
  { return (idx < np()) ? smat(idx, np()) : 0.0; }
  double ydys(int iidx, int jidx) const // yとdyの積和・平方和を返す
  { return (iidx < np() && jidx < np()) ? smat(iidx, jidx) : 0.0; }
  int xlyn() const              // xとlyの例数を返す
  { return statcxly.n(); }
  double xlymx() const          // xの平均値を返す
  { return statcxly.mx(); }
  double xlymly() const         // lyの平均値を返す
  { return statcxly.my(); }
  double xlyvxx() const         // xの不偏分散を返す
  { return statcxly.vxx(); }
  double xlyvlyly() const       // lyの不偏分散を返す
  { return statcxly.vyy(); }
  double xlysxx() const         // xの平方和を返す
  { return statcxly.sxx(); }
  double xlyslyly() const       // lyの平方和を返す
  { return statcxly.syy(); }
  double xlysxly() const        // xとlyの積和を返す
  { return statcxly.sxy(); }
  double y00() const            // y0初期値を返す
  { return v_y00; }
  double c0() const             // c初期値を返す
  { return v_c0; }
  double m0() const             // M初期値を返す
  { return v_m0; }
  double y0() const             // y0を返す
  { return v_y0; }
  double c() const              // cを返す
  { return v_c; }
  double m() const              // Mを返す
  { return v_m; }
  restype resregc() const       // 微分回帰分析の結果フラッグを返す
  { return v_resregc; }
  restype resvar() const       // 処理結果のフラッグを返す
  { return v_resvar; }
};

void varclass::cal(ldata& datalst, option& opt)
// 微分ロジット分析を実行する
{
  v_y00 = yrange.first * 0.5;  v_m0 = yrange.second;
  double b = (statcxy.sxx() > 0.0) ? statcxy.sxy() / statcxy.sxx() : 0.0;
  v_c0 = (sgn(v_m0) == 0) ? 0.0 : 8.0 * b / v_m0;
  v_y0 = v_y00;  v_c = v_c0;  v_m = v_m0; // 回帰直線によって仮初期値を設定する
  int n0 = cal_smat(datalst, opt);
  if (n0 < np())           // 微分回帰分析によって初期値を設定できない
  { v_resregc = reserr; }
  else
  {
    size_t dyidx = smat.r() - 1;
    double b1 = smat(0, dyidx), b2 = smat(1, dyidx);
    if (sgn(b1 * b2) >= 0)      // 不適解
    { v_resregc = reserr; }
    else
    {
      double c = b1, m = -b1 / b2, y0 = 0.0;
      v_resregc = cal_y0(y0, c, m, datalst, opt);
      if (v_resregc == resok)
      { v_y0 = v_y00 = y0;  v_c = v_c0 = c;  v_m = v_m0 = m; }
    }
  }
  v_resvar = cal_reglog(datalst, opt);
  if (v_resvar == reserr) return; // 微分ロジット分析できず
  if (v_noutd > 0 && opt.outd() == outdcomp) v_resvar = cal_replog(datalst, opt);
}

int varclass::cal_smat(ldata& datalst, option& opt)
// 微分回帰分析用に積和行列を計算する
{
  datalst.sort(data_less(0));   // datalstをx順にソートする
  size_t dim = 1 + np();        // 回帰式は定数項なし
  size_t dyidx = dim - 1;
  smat.resize(dim, dim);
  vecdbl fdatavec(dim);
  double x1 = 0.0, x2 = 0.0, y1 = 0.0, y2 = 0.0;
  int nxy = 0;
  ldata_itr ditrpre = datalst.begin();
  for (ldata_itr ditr = datalst.begin(); ditr != datalst.end(); ++ditr)
  {
    double x = ditr->fdata(0), y = ditr->fdata(1);
    ditr->resize(3);            // dyを追加する
    ditr->assign(2, BLANKf);
    if (x != BLANKf && y != BLANKf)
    {
      statcxy.sigma(x, y);
      if (x < xrange.first) xrange.first = x;
      if (x > xrange.second) xrange.second = x;
      if (y < yrange.first) yrange.first = y;
      if (y > yrange.second) yrange.second = y;
      if (nxy < 1)
      { x1 = x;  y1 = y;  nxy = 1; }
      else if (nxy == 1)
      {
        x2 = x;  y2 = y;  nxy = 2;
        ditrpre = ditr;
      }
      else
      {
        double dx1 = x2 - x1, dx2 = x - x2;
        double w1 = (fabs(dx1) < FDeps) ? 0.0 : 1.0 / dx1;
        double w2 = (fabs(dx2) < FDeps) ? 0.0 : 1.0 / dx2;
        double b1 = (y2 - y1) * w1, b2 = (y - y2) * w2;
        w1 = fabs(w1);  w2 = fabs(w2);
        double dy = (fabs(w1 + w2) < FDeps) ? 0.0 : (b1 * w1 + b2 * w2) / (w1 + w2);
        double fd = y2;
        statcydy.sigma(y2, dy);
        for (int i = 0; i < np(); ++i)
        {
          fdatavec[i] = fd;
          fd *= y2;
       }
        fdatavec[dyidx] = dy;
        for (size_t i = 0; i < dim; ++i)
        {
          fd = fdatavec[i];
          for (size_t j = i; j < dim; ++j) smat.add(i, j, fd * fdatavec[j]);
        }
        ditrpre->assign(2, dy); // dyを1つ前のケースに保存する
        x1 = x2;  y1 = y2;
        x2 = x;  y2 = y;  ++nxy;
        ditrpre = ditr;
      }
    }
  }
  for (size_t i = 0; i < dyidx; ++i)
  {
    for (size_t j = i + 1; j <= dyidx; ++j) smat.assign(j, i, smat(i, j));
  }
  for (size_t i = 0; i < dyidx; ++i) smat.sweep(i);
  return statcydy.n();          // 有効なデータ数を返す
}

restype varclass::cal_y0(double& y0, double c, double m, ldata& datalst, option& opt)
// y0をロジット分析によって求める
{
  for (auto& cdatai : datalst)
  {
    double x = cdatai.fdata(0), y = cdatai.fdata(1);
    if (y != BLANKf && x != BLANKf)
    {
      double ly = 0.0;
      if (0.0 < y && y < m)     // 範囲内データ
      {
        ly = log(y / (m - y));
        statcxly.sigma(x, ly);
      }
      else if (opt.outd() != outdexc) // 範囲外データ除外以外
      {
        if (y <= 0.0) ly = log(EPSly / (1.0 - EPSly)); // y0は近似法で求める
        else ly = log((1.0 - EPSly) / EPSly);
        statcxly.sigma(x, ly);
      }
    }
  }
  if (statcxly.n() < 1) return reserr; // ロジットがない時
  double ly0 = statcxly.my() - c * statcxly.mx();
  y0 = m / (1.0 + exp(-ly0));
  return resok;                 // y0計算完了
}

restype varclass::cal_reglog(ldata& datalst, option& opt)
// 微分ロジット分析によってパラメータを推定する(Mは更新しない)
{
  statcxly.clear();             // statcxlyを再使用するので初期化する
  double ly0 = log(v_y0 / (v_m - v_y0)); // ロジット回帰式の切片
  v_noutd = 0;                         // 範囲外データ数
  for (auto& cdatai : datalst)
  {
    double x = cdatai.fdata(0), y = cdatai.fdata(1);
    if (y != BLANKf && x != BLANKf)
    {
      double ly = ly0 + v_c * x;  // ロジット推定値(範囲外データ補完用)
      if (0.0 < y && y < v_m)     // 範囲内データ
      {
        ly = log(y / (v_m - y));
        statcxly.sigma(x, ly);
      }
      else                      // 範囲外データ
      {
        ++v_noutd;
        if (opt.outd() == outdapp) // 近似
        {
          if (y <= 0.0) ly = log(EPSly / (1.0 - EPSly));
          else ly = log((1.0 - EPSly) / EPSly);
          statcxly.sigma(x, ly);
        }
        else if (opt.outd() == outdcomp) // 繰り返し補完
        {
          if (isnan(ly) == false && -FDmax < ly && ly < FDmax) statcxly.sigma(x, ly);
        }
      }
    }
  }
  if (statcxly.n() < 2) return reserr; // ロジット分析不可能
  double mx = statcxly.mx(), my = statcxly.my();
  double sxx = statcxly.sxx(), sxy = statcxly.sxy();
  v_c = (sgn(sxx) > 0) ? sxy / sxx : 0.0;
  ly0 = my - v_c * mx;
  v_y0 = v_m / (1.0 + exp(-ly0));
  return resok;                 // ロジット分析完了
}

restype varclass::cal_replog(ldata& datalst, option& opt)
// 反復微分ロジット分析によってパラメータを更新する
{
  double syypre = statcxly.syy();
  const double Seps = 1.0e-5 * fmax(5.0, syypre);
  restype resreg = resnon;
  v_nr = 0;
  do
  {
    ++v_nr;
    resreg = cal_reglog(datalst, opt);
    if (resreg == resok)
    {
      double syy = statcxly.syy();
      if (fabs(syy - syypre) < Seps) resreg = resok; // 収束
      else if (v_nr >= RMAXdef) resreg = resrmax; // 反復計算数最大値
      else { syypre = syy;  resreg = resnon; }
    }
  } while (resreg == resnon);
  return resreg;                // 反復ロジット分析の結果を返す
}

restype varclass::cal_mci(double& mlow, double& mup, double cc0) const
// Mの信頼区間を計算する(Newton法)
{
  double t = ptot(1.0 - cc0, ydydfr());
  double v0 = 0.0, y = m()  - cal_dycr(m(), t, false), dey = 1.0, slope = 1.0;
  double vpre = cal_dyest(m() * 0.5);
  const double Rerror = FDeps * fabs(vpre);
  const double Ydelta = 10.0 * FDeps * fabs(m());
  constexpr int Rmax = RMAXdef / 2;
  int nr = 0;
  int sgnb2 = sgn(ydyb(1));     // 通常は-1
  mlow = mup = m();
  restype resmci = resnon;
  do
  {
    vpre = cal_dyest(y) + sgnb2 * cal_dycr(y, t, false);
    if (fabs(v0 - vpre) < Rerror)
    { resmci = resok; }
    else if (nr >= Rmax)
    { resmci = resrmax; }
    else
    {
      slope = ydyb(0) + 2.0 * ydyb(1) * y;
      if (fabs(slope) < Ydelta) resmci = reserr;
      else { y += (v0 - vpre) / slope;  ++nr; }
    }
  } while (resmci == resnon);
  if (resmci == reserr || m() < y) return reserr; // 下限サーチできず
  mlow = y;
  y = m() + (m() - mlow);  nr = 0;
  resmci = resnon;
  do
  {
    vpre = cal_dyest(y) - sgnb2 * cal_dycr(y, t, false);
    if (fabs(v0 - vpre) < Rerror)
    { resmci = resok; }
    else if (nr >= Rmax)
    { resmci = resrmax; }
    else
    {
      slope = ydyb(0) + 2.0 * ydyb(1) * y;
      if (fabs(slope) < Ydelta) resmci = reserr;
      else { y += (v0 - vpre) / slope;  ++nr; }
    }
  } while (resmci == resnon);
  if (resmci == reserr || y < mlow || FDmax < y) return reserr; // 上限サーチできず
  mup = y;
  return resmci;                // サーチ完了
}

double varclass::cal_dyest(double y) const // dy推定値を返す
{
  double yp = y;
  double dyest = 0.0;
  for (int i = 0; i < np(); ++i)
  {
    dyest += ydyb(i) * yp;
    yp *= y;
  }
  return dyest;                 // dy推定値を返す
}

double varclass::cal_dycr(double y, double x, bool plim) const
// dyの限界幅を返す
{
  vec<double> yv(np() + 1);
  double yp = y;
  for (int i = 0; i < np(); ++i)
  {
    yv.assign(i, yp);
    yp *= y;
  }
  yv.assign(np(), 0.0);
  vec<double> tv(yv.trn() * smat);
  double w = inp(tv, yv);
  if (plim) w += 1.0;           // 予測限界幅
  return x * sqrt(ydyvr() * w); // vの限界幅を返す
}

class dlaclass                  // dlaメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout);
  // 分析結果のタイトルとデータ識別情報を出力する
  void put_stat(ostream& fout, const varclass& var) const;
  // 基礎統計量を出力する
  void put_regc(ostream& fout, option& opt, const varclass& var) const;
  // yとdyの2次回帰分析の結果を出力する
  void put_coef2(ostream& fout, option& opt, const varclass& var) const;
  // 2次回帰式の偏回帰係数とパラメータを出力する
  void put_reglog(ostream& fout, option& opt, const varclass& var);
  // ロジット回帰分析の結果を出力する
  void put_coef1(ostream& fout, option& opt, double ly0, double vb0, double c, double vb1, int dfr, double m) const;
  // 1次回帰式の回帰係数を出力する
  void put_dns(option& opt, const varclass& var);
  // xとyとy推測値とlyとly推測値をデータファイル形式で出力する
  void put_dnsid(ofstream& fout);
  // DANSデータのタイトルとデータ識別情報を出力する
  void plot_ngp(option& opt, const varclass& var); // Ngraph用データ出力
  void plot_ngpid(ofstream& fout);
  // Ngraph用データのタイトルとデータ識別情報を出力する
  void plot_ngpxy(ostream& fout); // Ngraph用にxとyのデータを出力する
  void plot_ngplog(ostream& fout, option& opt, const varclass& var) const;
  // Ngraph用にxとyのロジスティック曲線を出力する
  int set_xrange(double& xmin, double& xmax, double& dx, option& opt, const varclass& var) const;
  void plot_ngploglim(ostream& fout, option& opt, const varclass& var) const;
  // Ngraph用にxとyのロジスティック曲線の信頼＆予測限界を出力する
  void plot_loglim(ostream& fout, const string& title, option& opt, const varclass& var, double ccpr, double tv, double vc, double vb, double mx) const;
  // Ngraph用にロジスティック曲線の信頼／予測限界を出力する
  void plot_ngpydy(ostream& fout);
  // Ngraph用にyとdyのデータを出力する
  void plot_ngpregc(ostream& fout, option& opt, const varclass& var) const;
  // Ngraph用にyとdyの2次回帰曲線を出力する
  void plot_ngpregclim(ostream& fout, option& opt, const varclass& var) const;
  // Ngraph用にyとdyの2次回帰曲線の信頼＆予測限界を出力する
  void plot_regclim(ostream& fout, option& opt, const varclass& var, double ccpr, double tv, bool plim) const;
  // Ngraph用にyとdyの2次回帰曲線の信頼／予測限界を出力する
  void plot_ngpxdy(ostream& fout);
  // Ngraph用にxとdyのデータを出力する
  void plot_ngpregd(ostream& fout, option& opt, const varclass& var) const;
  // Ngraph用にxとdyの微分ロジスティック曲線を出力する
  void plot_ngpxly(ostream& fout, const varclass& var);
  // Ngraph用にxとlyのデータを出力する
  void plot_ngpregl(ostream& fout, option& opt, const varclass& var) const;
  // Ngraph用にxとlyのロジット直線を出力する
  void plot_ngpregllim(ostream& fout, option& opt, const varclass& var) const;
  // Ngraph用にxとlyのロジット直線の信頼＆予測限界を出力する
  void plot_regllim(ostream& fout, const string& title, option& opt, const varclass& var, double ccpr, double tv, double vc, double vb, double mx) const;
  // Ngraph用にロジット直線の信頼／予測限界を出力する
public:
  dlaclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~dlaclass()                   // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み累積する
  void put(option& opt);                     // 集計結果を出力する
};

void dlaclass::put_resid(ostream& fout)
// 分析結果のタイトルとデータ識別情報を出力する
{
  puttitle(fout, Title);
  fout << "\n";
  id.put_title(fout);
  id.put_stra(fout);
  fout << "説明変数x:" << id.item(0) << "\n";
  fout << "目的変数y:" << id.item(1) << endl;
}

void dlaclass::put_stat(ostream& fout, const varclass& var) const
// 基礎統計量を出力する
{
  constexpr int LENline = 70; // 区切り線の長さ
  fout << "\n・変数の基礎統計量\n";
  putline(fout, LENline);
  fout << "x :";  putmsd(fout, var.xyn(), var.xymx(), var.xyvxx());
  fout << "y :";  putmsd(fout, var.xyn(), var.xymy(), var.xyvyy());
  fout << "dy:";  putmsd(fout, var.ydyn(), var.ydymdy(), var.ydyvdydy());
  putline(fout, LENline);
}

void dlaclass::put_regc(ostream& fout, option& opt, const varclass& var) const
// yとdyの2次回帰分析の結果を出力する
{
  fout << "\n・yとdyの2次回帰分析  2次回帰モデル:dy=b1・y^1+b2・y^2\n";
  size_t dyidx = var.np();      // 次数(2)
  int n = var.ydyn();
  double st = var.ydysdydy() + (double)n * var.ydymdy() * var.ydymdy();
  double sr = var.ydysr();
  int dfr = var.ydydfr(), dft = var.ydydft();
  if (dfr < 1 || sr < FDeps) sr = 0.0;
  int dfb = var.np();
  double sb = st - sr;
  if (sb < FDeps) sb = 0.0;
  put_coef2(fout, opt, var);
  if (var.resregc() == reserr)
  { fout << "※計算エラー発生！" << endl; }
  else
  {
    double ly0 = (sgn(var.y00()) > 0 && sgndiff(var.m0(), var.y00()) > 0) ? log(var.y00() / (var.m0() - var.y00())) : 0.0;
    double r2 = (st > 0.0) ? sb / st : 0.0;
    fout << "ロジットモデルの切片:y0=" << var.y00() << "  ly0=" << ly0 << "\n";
    fout << "重寄与率(決定係数):R2乗=" << r2 << "  重相関係数:R=" << sqrt(r2) << "\n\n";
    putatab(fout, "回帰", sb, dfb, sr, dfr, st, dft);
  }
}

void dlaclass::put_coef2(ostream& fout, option& opt, const varclass& var) const
// 2次回帰式の偏回帰係数とパラメータを出力する
{
  constexpr int Llen = 64;      // 区切り線の長さ
  fout << "記号    偏回帰係数      標準誤差";
  fout << setw(3) << opt.cc() * 100.0 << "%信頼区間：下限          上限\n";
  putline(fout, Llen);
  int idx = 0, offset = 1;
  vecdbl blow(var.np());
  vecdbl bup(var.np());
  double vbm = 0.0, ci = 0.0;
  for (int i = 0; i < var.np(); ++i)
  {
    double bi = var.ydyb(i), sii = var.ydys(i, i);
    double vbi = var.ydyvr() * sii;
    fout << "b" << left << setw(3) << i + 1;
    fout << right << setw(14) << bi << setw(14) << sqrt(vbi);
    ci = calmcci(mcfis, opt.cc(), 1, 1, var.ydydfr(), 1, vbi);
    blow[i] = bi - ci;  bup[i] = bi + ci;
    fout << setw(18) << blow[i] << setw(14) << bup[i] << "\n";
    vbm += vbi;
  }
  putline(fout, Llen);
  if (var.resregc() == resok)   // 微分回帰分析によって初期値を設定できた時
  {
    double clow = blow[0], cup = bup[0], mlow = 0.0, mup = 0.0;
    fout << "c   " << right << setw(14) << var.c0() << setw(32) << clow << setw(14) << cup << "\n";
    if (var.cal_mci(mlow, mup, opt.cc()) == reserr)
    {
      ci = calmcci(mcfis, opt.cc(), 1, 1, var.ydydfr(), 1, vbm);
      mlow = var.m0() - ci;  mup = var.m0() + ci;
    }
    fout << "M   " << right << setw(14) << var.m0() << setw(32) << mlow << setw(14) << mup << "\n";
  }
  else
  {
    fout << "c   " << right << setw(14) << var.c0() << " (1次回帰式によって推定)\n";
    fout << "M   " << right << setw(14) << var.m0() << " (1次回帰式によって推定)\n";
  }
  putline(fout, Llen);
}

void dlaclass::put_reglog(ostream& fout, option& opt, const varclass& var)
// ロジット回帰分析結果を出力する
{
  constexpr int LENline = 70; // 区切り線の長さ
  fout << "\n・ロジット分析用変数の基礎統計量\n  ";
  if (opt.outd() == outdexc)
  { fout << "範囲外データ除外:範囲外データ数=" << var.noutd(); }
  else if (opt.outd() == outdapp)
  { fout << "範囲外データ近似:範囲外データ数=" << var.noutd(); }
  else
  {
    fout << "範囲外データ反復補完:範囲外データ数=" << var.noutd() << "  反復回数=" << var.nr();
    if (var.resvar() == resrmax) fout << "(最大回数)";
  }
  fout << "\n";
  putline(fout, LENline);
  fout << "x :";  putmsd(fout, var.xlyn(), var.xlymx(), var.xlyvxx());
  fout << "ly:";  putmsd(fout, var.xlyn(), var.xlymly(), var.xlyvlyly());
  putline(fout, LENline);
  fout << "\n・ロジット分析  ロジット回帰モデル:ly=ly0+c・x\n";
  fout << "  ロジット:ly=ln(y/{M-y})  ly0=ln(y0/{M-y0})  M=" << var.m() << "(固定)\n";
  int dft = var.xlyn() - 1, dfb = 1;
//  if (opt.outd() == outdcomp) dft -= var.noutd(); // 補完値は自由度を減らす
  int dfr = dft - dfb;
  double mx = var.xlymx(), my = var.xlymly();
  double sxx = var.xlysxx(), st = var.xlyslyly(), sxy = var.xlysxly();
  double b1 = (sxx > 0.0) ? sxy / sxx : 0.0; // =c
  double ly0 = my - b1 * mx;
  double sb = sxy * b1;
  double sr = st - sb;
  const double Seps = DBL_EPSILON * fmax(5.0, st);
  if (dfr < 1 || sr < Seps) sr = 0.0;
  double vr = (dfr > 0) ? sr / (double)dfr : 0.0;
  double vb0 = vr * (1.0 / (double)(dft + 1) + mx * mx / sxx);
  double vb1 = (sxx > 0.0) ? vr / sxx : 0.0;
  put_coef1(fout, opt, ly0, vb0, var.c(), vb1, dfr, var.m());
  if (var.resvar() == reserr)
  { fout << "※計算エラー発生！" << endl; }
  else
  {
    double r2 = (st > 0.0) ? sb / st : 0.0;
    if (opt.res("y"))
    {
      double y = opt.y();
      if (opt.y() == BLANKf) y = var.m() * 0.5;
      else if (opt.y() == DBL_MAX) y = var.xymy();
      else if (opt.y() >= var.m()) y = var.m() * 0.5;
      double ly = log(y / (var.m() - y));
      double x = (ly - ly0) / var.c();
      fout << "y=" << y << "の時:ly=" << ly << "  x=" << x << "\n";
    }
    fout << "重寄与率(決定係数):R2乗=" << r2 << "  重相関係数:R=" << sqrt(r2) << "\n\n";
    putatab(fout, "回帰", sb, dfb, sr, dfr, st, dft);
  }
}

void dlaclass::put_coef1(ostream& fout, option& opt, double ly0, double vb0, double c, double vb1, int dfr, double m) const
// 1次回帰式の回帰係数を出力する
{
  constexpr int Llen = 64;      // 区切り線の長さ
  fout << "記号      回帰係数      標準誤差";
  fout << setw(3) << opt.cc() * 100.0 << "%信頼区間：下限          上限\n";
  putline(fout, Llen);
  fout << "ly0 " << right << setw(14) << ly0 << setw(14) << sqrt(vb0);
  double ci = calmcci(mcfis, opt.cc(), 1, 1, dfr, 1, vb0);
  double ly0l = ly0 - ci, ly0u = ly0 + ci;
  fout << setw(18) << ly0l << setw(14) << ly0u << "\n";
  fout << "c   " << right << setw(14) << c << setw(14) << sqrt(vb1);
  ci = calmcci(mcfis, opt.cc(), 1, 1, dfr, 1, vb1);
  fout << setw(18) << c - ci << setw(14) << c + ci << "\n";
  putline(fout, Llen);
  double y0 = m / (1.0 + exp(-ly0));
  double y0l = m / (1.0 + exp(-ly0l)), y0u = m / (1.0 + exp(-ly0u));
  fout << "y0  " << right << setw(14) << y0 << setw(14) << " ";
  fout << setw(18) << y0l << setw(14) << y0u << "\n";
  putline(fout, Llen);
}

void dlaclass::put_dns(option& opt, const varclass& var)
// xとyとy推測値とlyとly推測値をデータファイル形式で出力する
{
  ofstream fout;
  if (ofopen(fout, (opt.dnsname()).c_str()) == false) return; // ファイルオープン失敗
  put_dnsid(fout);
  datac tdata(5);
  double ly0 = log(var.y0() / (var.m() - var.y0()));
  for (auto& cdatai : datalst)
  {
    double x = cdatai.fdata(0), y = cdatai.fdata(1);
    tdata.assign(0, x);  tdata.assign(1, y);
    if (y == BLANKf || x == BLANKf)
    {
      tdata.assign(2, BLANKs);  tdata.assign(3, BLANKs);
      tdata.assign(4, BLANKs);
    }
    else
    {
      double lye = ly0 + var.c() * x; // ロジット推測値
      double ye = var.m() / (1.0 + exp(-lye));
      tdata.assign(2, ye);  // y推測値
      double ly = lye;
      if (0.0 < y && y < var.m()) ly = log(y / (var.m() - y)); // 範囲内
      else ly = BLANKf;
      tdata.assign(3, ly);
      tdata.assign(4, lye);
    }
    tdata.put(fout);
  }
}

void dlaclass::put_dnsid(ofstream& fout)
// DANSデータのタイトルとデータ識別情報を出力する
{
  idd tid(id.title(), 5);
  ostringstream ostr;
  id.put_title(fout, formtag);
  id.put_stra(fout, formtag);
  tid.assign_item(0, id.item(0));
  tid.assign_item(1, id.item(1));
  ostr << id.item(1) << "の推測値";
  tid.assign_item(2, ostr.str());
  ostr.str("");
  ostr << id.item(1) << "のロジット";
  tid.assign_item(3, ostr.str());
  tid.assign_item(4, "ロジットの推測値");
  tid.put_item(fout, formtag);
}

void dlaclass::plot_ngp(option& opt, const varclass& var)
// Ngraph用データ出力
{
  ofstream fout;
  if (ofopen(fout, (opt.ngpname()).c_str()) == false) return; // ファイルオープン失敗
  bool delml = false;
  if (opt.ngp("logc"))          // xとyとロジスティック曲線を出力する
  {
    fout << "-s1 -dline\n";
    plot_ngpid(fout);
    plot_ngpxy(fout);
    fout << NGPdelm << "\n";
    plot_ngplog(fout, opt, var);
    fout << NGPdelm << "\n";
    plot_ngploglim(fout, opt, var);
    delml = true;
  }
  if (opt.ngp("regc"))          // yとdyと2次回帰曲線を出力する
  {
    if (delml)                  // 区切り行を出力する
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline\n";
      delml = true;
    }
    plot_ngpid(fout);
    plot_ngpydy(fout);
    fout << NGPdelm << "\n";
    plot_ngpregc(fout, opt, var);
    fout << NGPdelm << "\n";
    plot_ngpregclim(fout, opt, var);
  }
  if (opt.ngp("regd"))          // xとdyと微分ロジスティック曲線を出力する
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline\n";
      delml = true;
    }
    plot_ngpid(fout);
    plot_ngpxdy(fout);
    fout << NGPdelm << "\n";
    plot_ngpregd(fout, opt, var);
  }
  if (opt.ngp("regl"))          // xとlyとロジット直線を出力する
  {
    if (delml) fout << NGPdelm << "\n";
    else fout << "-s1 -dline\n";
    plot_ngpid(fout);
    plot_ngpxly(fout, var);
    fout << NGPdelm << "\n";
    plot_ngpregl(fout, opt, var);
    fout << NGPdelm << "\n";
    plot_ngpregllim(fout, opt, var);
  }
}

void dlaclass::plot_ngpid(ofstream& fout)
// Ngraph用データのタイトルとデータ識別情報を出力する
{
  fout << "%";  puttitle(fout, Title);
  fout << "%\n";
  id.put_title(fout, formngp);
  id.put_stra(fout, formngp);
  fout << "%説明変数x:" << id.item(0) << "\n";
  fout << "%目的変数y:" << id.item(1) << endl;
}

void dlaclass::plot_ngpxy(ostream& fout)
// Ngraph用にxとyのデータを出力する
{
  fout << "%説明変数xと目的変数yのデータ\n";
  for (auto& cdatai : datalst)
  {
    double x = cdatai.fdata(0), y = cdatai.fdata(1);
    if (y != BLANKf && x != BLANKf) fout << x << " " << y << endl;
  }
}

void dlaclass::plot_ngplog(ostream& fout, option& opt, const varclass& var) const
// Ngraph用にxとyのロジスティック曲線を出力する
{
  fout << "%説明変数xと目的変数yのロジスティック曲線\n";
  double xmin = 0.0, xmax = 1.0, dx = 1.0;
  int rd = set_xrange(xmin, xmax, dx, opt, var);
  double x = xmin, ly0 = log(var.y0() / (var.m() - var.y0()));
  for (int i = 0; i <= rd; ++i, x += dx)
  { fout << x << " " << var.m() / (1.0 + exp(-ly0 - var.c() * x)) << endl; }
}

int dlaclass::set_xrange(double& xmin, double& xmax, double& dx, option& opt, const varclass& var) const
// xの最小値と最大値とdxを返す
{
  if (opt.xmin() != BLANKf && opt.xmax() != BLANKf) // X軸の範囲指定あり
  { xmin = opt.xmin();  xmax = opt.xmax(); }
  else
  {
    constexpr double EPSYdef = 0.001; // yの微小値省略時値
    double ly0 = log(var.y0() / (var.m() - var.y0()));
    double y = var.m() * EPSYdef; // y最小値
    xmin = (log(y / (var.m() - y)) - ly0) / var.c();
    if (xmin > var.xmin()) xmin = var.xmin();
    if (sgn(xmin) < 0) xmin = 0.0;
    y = var.m() * (1.0 - EPSYdef); // y最大値
    xmax = (log(y / (var.m() - y)) - ly0) / var.c();
    if (xmax < var.xmax()) xmax = var.xmax();
  }
  dx = 1.0;
  int rd = (int)((xmax - xmin) / dx);
  if (rd < RDdef)               // 分割数が省略時値未満の時
  {
    int r = (int)((double)RDdef / (double)rd + 0.9);
    rd *= r;  dx /= (double)r;
  }
  else if (rd > RDdef * 10)     // 分割数が省略時値の10倍以上の時
  {
    int r = (int)((double)rd / (double)(RDdef * 10) + 0.9);
    rd /= r;  dx *= (double)r;
  }
  return rd;                    // 分割数を返す
}

void dlaclass::plot_ngploglim(ostream& fout, option& opt, const varclass& var) const
// Ngraph用にxとyのロジスティック曲線の信頼＆予測限界を出力する
{
  int dft = var.xlyn()- 1, dfb = 1;
  if (opt.outd() == outdcomp) dft -= var.noutd(); // 補完値は自由度を減らす
  int dfr = dft - dfb;
  double mx = var.xlymx(), my = var.xlymly();
  double sxx = var.xlysxx(), st = var.xlyslyly(), sxy = var.xlysxly();
  double b1 = (sxx > 0.0) ? sxy / sxx : 0.0;
  double sb = sxy * b1;
  double sr = st - sb;
  const double Seps = DBL_EPSILON * fmax(5.0, st);
  if (dfr < 1 || sr < Seps) sr = 0.0;
  double vr = (dfr > 0) ? sr / (double)dfr : 0.0;
  double cf = 1.0 / (double)(dft + 1), rf = 1.0 + 1.0 / (double)(dft + 1);
  double vb = (sxx > 0.0) ? vr / sxx : 0.0;
  double tv = ptot(1.0 - opt.cc(), dfr);
  plot_loglim(fout, "信頼限界", opt, var, opt.cc() * 100.0, tv, vr * cf, vb, mx);
  fout << NGPdelm << "\n";
  plot_loglim(fout, "予測限界", opt, var, opt.cc() * 100.0, tv, vr * rf, vb, mx);
}

void dlaclass::plot_loglim(ostream& fout, const string& title, option& opt, const varclass& var, double ccpr, double tv, double vc, double vb, double mx) const
// Ngraph用にロジスティック曲線の信頼／予測限界を出力する
{
  fout << "%説明変数xと目的変数yのロジスティック曲線の" << ccpr << "%" << title << "上限\n";
  double xmin = 0.0, xmax = 1.0, dx = 1.0;
  int rd = set_xrange(xmin, xmax, dx, opt, var);
  double ly0 = log(var.y0() / (var.m() - var.y0()));
  double x = xmin, dm = 0.0, ly = 0.0, y = 0.0;
  for (int i = 0; i <= rd; ++i, x += dx)
  {
    dm = x - mx;
    ly = ly0 + var.c() * x + tv * sqrt(vc + vb * dm * dm);
    fout << x << " " << var.m() / (1.0 + exp(-ly)) << "\n";
  }
  fout << NGPdelm2 << "\n"; 
  fout << "%説明変数xと目的変数yのロジスティック曲線の" << ccpr << "%" << title << "下限\n";
  x = xmin;
  for (int i = 0; i <= rd; ++i, x += dx)
  {
    dm = x - mx;
    ly = ly0 + var.c() * x - tv * sqrt(vc + vb * dm * dm);
    fout << x << " " << var.m() / (1.0 + exp(-ly)) << endl;
  }
}

void dlaclass::plot_ngpydy(ostream& fout)
// Ngraph用にyとdyのデータを出力する
{
  fout << "%目的変数yとその変化量dyのデータ\n";
  for (auto& cdatai : datalst)
  {
    double x = cdatai.fdata(0), y = cdatai.fdata(1), dy = cdatai.fdata(2);
    if (x != BLANKf && y != BLANKf && dy != BLANKf) fout << y << " " << dy << endl;
  }
}

void dlaclass::plot_ngpregc(ostream& fout, option& opt, const varclass& var) const
// Ngraph用にyとdyの2次回帰曲線を出力する
{
  fout << "%目的変数yとその変化量dyの2次回帰曲線\n";
  double yrange = var.m() - var.y0(), ymg = 0.1;
  double ymin = var.y0() - yrange * ymg, ymax = var.m() + yrange * ymg;
  double dy = (ymax - ymin) / (double)RDdef;
  double y = ymin;
  for (int i = 0; i <= RDdef; ++i, y += dy)
  { fout << y << " " << var.cal_dyest(y) << endl; }
}

void dlaclass::plot_ngpregclim(ostream& fout, option& opt, const varclass& var) const
// Ngraph用にyとdyの2次回帰曲線の信頼＆予測限界を出力する
{
  double tv = ptot(1.0 - opt.cc(), var.ydydfr());
  plot_regclim(fout, opt, var, opt.cc() * 100.0, tv, false);
  fout << NGPdelm << "\n";
  plot_regclim(fout, opt, var, opt.cc() * 100.0, tv, true);
}

void dlaclass::plot_regclim(ostream& fout, option& opt, const varclass& var, double ccpr, double tv, bool plim) const
// Ngraph用にyとdyの2次回帰曲線の信頼／予測限界を出力する
{
  fout << "%目的変数yとその変化量dyの2次回帰曲線の" << ccpr << "%";
  if (plim) fout << "予測限界";
  else fout << "信頼限界";
  fout << "上限\n";
  double yrange = var.m() - var.y0(), ymg = 0.1;
  double ymin = var.y0() - yrange * ymg, ymax = var.m() + yrange * ymg;
  double dy = (ymax - ymin) / (double)RDdef;
  double y = ymin, v = 0.0;
  for (int i = 0; i <= RDdef; ++i, y += dy)
  {
    v = var.cal_dyest(y) + var.cal_dycr(y, tv, plim);
    fout << y << " " << v << "\n";
  }
  fout << NGPdelm2 << "\n"; 
  fout << "%目的変数yとその変化量dyの2次回帰曲線の" << ccpr << "%";
  if (plim) fout << "予測限界";
  else fout << "信頼限界";
  fout << "下限\n";
  y = ymin;
  for (int i = 0; i <= RDdef; ++i, y += dy)
  {
    v = var.cal_dyest(y) - var.cal_dycr(y, tv, plim);
    fout << y << " " << v << endl;
  }
}

void dlaclass::plot_ngpxdy(ostream& fout)
// Ngraph用にxとdyのデータを出力する
{
  fout << "%説明変数xと目的変数変化量dyのデータ\n";
  for (auto& cdatai : datalst)
  {
    double x = cdatai.fdata(0), y = cdatai.fdata(1), dy = cdatai.fdata(2);
    if (x != BLANKf && y != BLANKf && dy != BLANKf) fout << x << " " << dy << endl;
  }
}

void dlaclass::plot_ngpregd(ostream& fout, option& opt, const varclass& var) const
// Ngraph用にxとdyの微分ロジスティック曲線を出力する
{
  fout << "%説明変数xと目的変数変化量dyの微分ロジスティック曲線\n";
  double xmin = 0.0, xmax = 1.0, dx = 1.0;
  int rd = set_xrange(xmin, xmax, dx, opt, var);
  double x = xmin;
  for (int i = 0; i <= rd; ++i, x += dx)
  { fout << x << " " << logden(x, var.m(), var.c(), var.y0()) << endl; }
}

void dlaclass::plot_ngpxly(ostream& fout, const varclass& var)
// Ngraph用にxとlyのデータを出力する
{
  fout << "%説明変数xと目的変数ロジットlyのデータ\n";
  for (auto& cdatai : datalst)
  {
    double x = cdatai.fdata(0), y = cdatai.fdata(1);
    if (y != BLANKf && x != BLANKf)
    {
      if (0.0 < y && y < var.m()) fout << x << " " << log(y / (var.m() - y)) << endl;
    }
  }
}

void dlaclass::plot_ngpregl(ostream& fout, option& opt, const varclass& var) const
// Ngraph用にxとlyのロジット直線を出力する
{
  fout << "%説明変数xと目的変数ロジットlyのロジット直線\n";
  double xmin = 0.0, xmax = 0.0, dx = 0.0;
  int rd = set_xrange(xmin, xmax, dx, opt, var);
  double ly0 = log(var.y0() / (var.m() - var.y0()));
  fout << xmin << " " << ly0 + var.c() * xmin << "\n";
  fout << xmax << " " << ly0 + var.c() * xmax << endl;
}

void dlaclass::plot_ngpregllim(ostream& fout, option& opt, const varclass& var) const
// Ngraph用にxとlyのロジット直線の信頼＆予測限界を出力する
{
  int dft = var.xlyn() - 1, dfb = 1;
  if (opt.outd() == outdcomp) dft -= var.noutd(); // 補完値は自由度を減らす
  int dfr = dft - dfb;
  double mx = var.xlymx(), my = var.xlymly();
  double sxx = var.xlysxx(), st = var.xlyslyly(), sxy = var.xlysxly();
  double b1 = (sxx > 0.0) ? sxy / sxx : 0.0;
  double sb = sxy * b1;
  double sr = st - sb;
  if (dfr < 1 || sr < FDeps) sr = 0.0;
  double vr = (dfr > 0) ? sr / (double)dfr : 0.0;
  double cf = 1.0 / (double)(dft + 1), rf = 1.0 + 1.0 / (double)(dft + 1);
  double vb = (sxx > 0.0) ? vr / sxx : 0.0;
  double tv = ptot(1.0 - opt.cc(), dfr);
  plot_regllim(fout, "信頼限界", opt, var, opt.cc() * 100.0, tv, vr * cf, vb, mx);
  fout << NGPdelm << "\n";
  plot_regllim(fout, "予測限界", opt, var, opt.cc() * 100.0, tv, vr * rf, vb, mx);
}

void dlaclass::plot_regllim(ostream& fout, const string& title, option& opt, const varclass& var, double ccpr, double tv, double vc, double vb, double mx) const
// Ngraph用にロジット直線の信頼／予測限界を出力する
{
  fout << "%説明変数xと目的変数ロジットlyのロジット直線の" << ccpr << "%" << title << "上限\n";
  double xmin = 0.0, xmax = 1.0, dx = 1.0;
  int rd = set_xrange(xmin, xmax, dx, opt, var);
  double ly0 = log(var.y0() / (var.m() - var.y0()));
  double x = xmin, dm = 0.0, ly = 0.0;
  for (int i = 0; i <= rd; ++i, x += dx)
  {
    dm = x - mx;
    ly = ly0 + var.c() * x + tv * sqrt(vc + vb * dm * dm);
    fout << x << " " << ly << "\n";
  }
  fout << NGPdelm2 << "\n"; 
  fout << "%説明変数xと目的変数ロジットlyのロジット直線の" << ccpr << "%" << title << "下限\n";
  x = xmin;
  for (int i = 0; i <= rd; ++i, x += dx)
  {
    dm = x - mx;
    ly = ly0 + var.c() * x - tv * sqrt(vc + vb * dm * dm);
    fout << x << " " << ly << endl;
  }
}

int dlaclass::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 dlaclass::put(option& opt) // 集計結果を出力する
{
  varclass var;                 // 計算用変数クラス
  var.cal(datalst, opt);
  if (opt.res("res"))           // 分析結果出力
  {
    put_resid(cout);
    if (opt.res("stat")) put_stat(cout, var);
    if (var.xyn() < 1)
    {
      cout << "※有効なデータなし！\n\n" << endl;
      return;
    }
    if (opt.res("regc")) put_regc(cout, opt, var);
    if (var.resvar() != reserr && opt.res("regl")) put_reglog(cout, opt, var);
    cout << "\n" << endl;
  }
  if (var.resvar() == reserr) return;
  if (opt.dnsname() != "") put_dns(opt, var); // データファイル出力
  if (opt.ngpname() != "") plot_ngp(opt, var); // Ngraphデータ出力
}

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