// roc.C Version 7.44  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-29 11:09:03 tac>
//
// 群項目と集計項目からROC分析を行う。
// カットオフポイントはROC曲線最適点、つまりROC曲線が感度1・特異度1の点
// (左上隅)に最も近い点か、感度と特異度の平均が最大となる点、
// 正診率が最大となる点または指定された点のどれかで計算。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -W roc roc.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 =            // プログラム使用法表示用
{
  "roc  Version 7.44            Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:roc  [<オプション> …]  [<データファイル> …]",
  "機能:最初の項目を群項目、次の項目を集計項目としてROC分析を行う",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定する。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -co[{low|up}]={a|c|r|y|<境界値>}:境界値(カットオフポイント)を指定する。",
  "    co:上限／下限を自動判定  colow:下限境界値  coup:上限境界値",
  "  -dnsname[=<出力ファイル名>]:DANSデータを出力するファイル名を指定する。",
  "    オプション省略時:DANSデータ出力せず  -dnsname:-dnsname=roc_dns.d",
  "    出力項目:診断率",
  "  -f=<ファイル名>:指定されたファイルからオプションを読み込む。",
  "  -g{1|2}=<群コード>:群1(疾患群)または群2(正常群)のコードを指定する。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示する。(-h:-h=man)",
  "  -ngp={roc|roci}[,…]:Ngraph用出力データを指定する。",
  "    roc:ROC曲線  roci:ROC曲線の信頼区間  オプション省略時:全て指定",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データを出力するファイル名を指定する。",
  "    オプション省略時:Ngraph用データ出力せず  -ngpname:-ngpname=roc_ngp.dat",
  "  -res=[-][{def|all|auc|or|roc|stat|uda}][,…]:解析結果の出力項目を指定する。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    auc:ROC曲線のAUC  or:オッズ比  roc:ROC分析  stat:基礎統計量",
  "    uda:単変量線形判別分析  オプション省略時:自動判定",
  "  -pi=<群1の事前確率>:群1の事前確率を指定する。(0＜事前確率＜1)"
};

const string Title = "=== ROC分析(Receiver Operating Characteristic analysis) ==="; // 表題
const string DNSNAMEdef = "roc_dns.d"; // 省略時DANSファイル名
const string NGPNAMEdef = "roc_ngp.dat"; // 省略時Ngraphファイル名
constexpr int RDdef = 50;       // 座標軸の分割数
constexpr double PIdef = 0.5;   // 省略時事前確率

enum cotype                     // カットオフポイントの種類
{
  conon,                        // 指定なし
  coac,                         // 正診率最大点
  cocm,                         // 最大χ2乗法
  coroc,                        // ROC曲線最適点
  coyi,                         // Youden index:(感度+特異度-1)最大点
  colim,                        // 境界値(上下限は自動判定)
  lowlim,                       // 下限境界値
  uplim,                        // 上限境界値
  coetc                         // その他
};

class option                    // オプション集合体
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  cotype v_co;                  // カットオフポイントの種類
  double v_lim;                 // 境界値
  string v_dnsname;             // DANSファイル名
  vecstr gvec;                  // 群コードベクトル
  string v_ngpname;             // Ngraphファイル名
  lstr ngplst;                  // Ngraph出力リスト
  mapsb ngpmap;                 // Ngraph出力の有無マップ
  double v_pi;                  // 群1の事前確率
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  double cc() const             // 信頼係数を返す
  { return v_cc; }
  cotype co() const             // カットオフポイントの種類を返す
  { return v_co; }
  double lim() const            // 境界値を返す
  { return v_lim; }
  const string& dnsname() const // DANSデータファイル名を返す
  { return v_dnsname; }
  const string& g(int idx) const // 群コードを返す
  { return gvec[idx]; }
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  bool ngp(const string& code); // Ngraph出力フラッグを返す
  double pi() const             // 群1の事前確率を返す
  { return v_pi; }
  bool res(const string& code); // 分析結果出力フラッグを返す
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int search_g(const string& gcode) const; // 群コードを検索する
  int check(size_t ni);         // オプションのチェックとセット
  int set_ngp();                // -ngpオプションのセット
  int set_res();                // -resオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_cc(CCdef), v_co(coroc), v_lim(BLANKf), gvec(2), v_dnsname(), v_ngpname(), ngplst(), ngpmap(), v_pi(BLANKf), reslst(), resmap()
// コンストラクタ
{
  if (c_pname == 0) return;     // 初期設定せず
  gvec[0] = "1";                // 群1のコードだけ初期設定
  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][0] == 'c' && opt[0][1] == 'o') // カットオフポイントの種類の指定
  {
    if (opt[0] == "co")
    {
      if (opt[1].empty() || opt[1] == "r") v_co = coroc;
      else if (opt[1] == "a") v_co = coac;
      else if (opt[1] == "c") v_co = cocm;
      else if (opt[1] == "y") v_co = coyi;
      else v_co = colim;
    }
    else if (opt[0] == "colow") v_co = lowlim;
    else if (opt[0] == "coup") v_co = uplim;
    else err = badopt(pname(), optstr, USAGE);
    if (v_co == colim || v_co == lowlim || v_co == uplim) // 指定境界値の時
    {
      v_lim = strtodf(opt[1]);
      if (v_lim == BLANKf) 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] == "f")       // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "g1")      // 群1のコード指定
  {
    if (opt[1].empty()) gvec[0].erase();
    else gvec[0] = opt[1];
  }
  else if (opt[0] == "g2")      // 群2のコード指定
  {
    if (opt[1].empty()) gvec[1].erase();
    else gvec[1] = opt[1];
  }
  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] == "pi")      // 群1の事前確率指定
  {
    v_pi = opt[1].empty() ? PIdef : strtodf(opt[1]);
    if (v_pi <= 0.0 || 1.0 <= v_pi) err = badopt(pname(), optstr, USAGE, 0.0, 1.0);
  }
  else if (opt[0] == "res")     // 分析結果の出力指定
  {
    if (strsplit(reslst, opt[1], 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::search_g(const string& gcode) const // 群コードを検索する
{
  if (gcode == BLANKs) return -1; // 欠損データはヒットさせない
  int idx = -1;
  if (gvec[0].empty())          // 群1のコード指定がない場合
  {
    idx = 0;
    if (gcode == gvec[1]) idx = 1;
  }
  else if (gcode == gvec[0])
  { idx = 0; }
  else if (gvec[1].empty())     // 群2のコード指定がない場合
  { idx = 1; }
  else if (gcode == gvec[1])
  { idx = 1; }
  return idx;                   // -1:ヒットせず,0:群1,1:群2
}

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["roc"] = ngpmap["roci"] = true;
  if (ngplst.size() > 0)
  {
    if (v_ngpname == "") return disperr(pname(), "-ngpnameオプションの指定が必要！");
    ngpmap["roc"] = ngpmap["roci"] = 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["auc"] = resmap["or"] = resmap["roc"] = resmap["stat"] = true;
  resmap["uda"] = false;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["auc"] = resmap["or"] = resmap["roc"] = resmap["stat"] = 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") // 省略時出力
      {
        resmap["def"] = resmap["auc"] = resmap["or"] = resmap["roc"] = resmap["stat"] = true;
      }
      else if (resstr == "all") // 全て出力
      {
        resmap["def"] = resmap["auc"] = resmap["or"] = resmap["roc"] = resmap["stat"] = resmap["uda"] = 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 diag                      // 診断率とROCクラス
{
private:
  cotype v_co;                  // 境界値の種類
  cotype v_colu;                // 境界値の上限／下限フラッグ
  double v_pi;                  // 群1の事前確率
  double v_lim;                 // 境界値
  int v_nt;                     // 合計例数
  int v_tp;                     // True-Positive例数
  int v_tn;                     // True-Negative例数
  int v_fp;                     // False-Positive例数
  int v_fn;                     // False-Negative例数
  double v_auc;                 // ROC曲線のAUC(曲線下面積)
  double v_aucv;                // ROC曲線のAUCの分散
  double v_u;                   // Mann-WhitneyのU値
  double v_uv;                  // Mann-WhitneyのU値の分散
  double v_nlim;                // 理論的境界値
  double v_nsn;                 // 理論的境界値による感度
  double v_nsp;                 // 理論的境界値による特異度
  double v_nac;                 // 理論的境界値による正診率
  double v_nppv;                // 理論的境界値によるPPV
  double v_nnpv;                // 理論的境界値によるNPV
  double v_nauc;                // 理論的AUC
  double v_d2;                  // マハラノビスの汎距離
  double v_ppr;                 // 誤判別確率
  vector<statv<double>> statvec; // 基礎統計量ベクトル
  void cal_roc(ldata& datalst); // ROC分析を行う
  void cal_dis();               // 単変量判別分析を行う
public:
  diag() : v_co(conon), v_colu(lowlim), v_pi(BLANKf), v_lim(BLANKf), v_nt(0), v_tp(0), v_tn(0), v_fp(0), v_fn(0), v_auc(BLANKf), v_aucv(BLANKf), v_u(BLANKf), v_uv(BLANKf), v_nlim(BLANKf), v_nsn(BLANKf), v_nsp(BLANKf), v_nac(BLANKf), v_nppv(BLANKf), v_nnpv(BLANKf), v_nauc(BLANKf), v_d2(BLANKf), v_ppr(BLANKf), statvec(2)
  // デフォールトコンストラクタ
  {}
  explicit diag(cotype co0, double pi0 = BLANKf, double lim0 = BLANKf) : v_co(co0), v_colu(lowlim), v_pi(pi0), v_lim(lim0), v_nt(0), v_tp(0), v_tn(0), v_fp(0), v_fn(0), v_auc(BLANKf), v_aucv(BLANKf), v_u(BLANKf), v_uv(BLANKf), v_nlim(BLANKf), v_nsn(BLANKf), v_nsp(BLANKf), v_nac(BLANKf), v_nppv(BLANKf), v_nnpv(BLANKf), v_nauc(BLANKf), v_d2(BLANKf), v_ppr(BLANKf), statvec(2)
  // コンストラクタ
  {}
  virtual ~diag()               // デストラクタ
  {}
  diag(const diag& diag0) : v_co(diag0.v_co), v_colu(diag0.v_colu), v_pi(diag0.v_pi), v_lim(diag0.v_lim), v_tp(diag0.v_tp), v_tn(diag0.v_tn), v_fp(diag0.v_fp), v_fn(diag0.v_fn), v_auc(diag0.v_auc), v_aucv(diag0.v_aucv), v_u(diag0.v_u), v_uv(diag0.v_uv), v_nlim(diag0.v_nlim), v_nsn(diag0.v_nsn), v_nsp(diag0.v_nsp), v_nac(diag0.v_nac), v_nppv(diag0.v_nppv), v_nnpv(diag0.v_nnpv), v_nauc(diag0.v_nauc), v_d2(diag0.v_d2), v_ppr(diag0.v_ppr), statvec(diag0.statvec)
  // コピーコンストラクタ
  {}
  const diag& operator=(const diag& diag0); // コピー代入演算子
  cotype co() const             // 境界値の種類を返す
  { return v_co; }
  cotype co(cotype co0)         // 境界値の種類を設定する
  { return (v_co = co0); }
  cotype colu() const           // 境界値の上限／下限フラッグを返す
  { return v_colu; }
  cotype colu(cotype colu0)    // 境界値の上限／下限フラッグを設定する
  { return (v_colu = colu0); }
  double pi()                   // 群1の事前確率を返す
  { return v_pi; }
  double pi(double pi0)         // 群1の事前確率を設定する
  { return (v_pi = pi0); }
  double lim() const            // 境界値を返す
  { return v_lim; }
  double lim(double lim0)       // 境界値を設定する
  { return (v_lim = lim0); }
  int nt() const                // 合計例数を返す
  { return v_nt; }
  int n1() const                // 群1の例数を返す
  { return statvec[0].n(); }
  int n2() const                // 群2の例数を返す
  { return statvec[1].n(); }
  int tp() const                // True-Positive例数を返す
  { return v_tp; }
  int tn() const                // True-Negative例数を返す
  { return v_tn; }
  int fp() const                // False-Positive例数を返す
  { return v_fp; }
  int fn() const                // False-Negative例数を返す
  { return v_fn; }
  double sn() const             // 感度を返す
  { return ((tp() + fn()) < 1) ? 0.0 : (double)(tp()) / (double)(tp() + fn()); }
  double sp() const             // 特異度を返す
  { return ((fp() + tn()) < 1) ? 0.0 : (double)(tn()) / (double)(fp() + tn()); }
  double ac() const             // 正診率を返す
  { return ((tp() + fn() + fp() + tn()) < 1) ? 0.0 : (double)(tp() + tn()) / (double)(tp() + fn() + fp() + tn()); }
  double plr() const;           // 陽性尤度比を返す
  double nlr() const;           // 陰性尤度比を返す
  double ppv();                 // 陽性予測値を返す
  double npv();                 // 陰性予測値を返す
  double auc() const            // ROC曲線のAUC(曲線下面積)を返す
  { return v_auc; }
  double aucv() const           // ROC曲線AUCの分散を返す
  { return v_aucv; }
  double u() const              // Mann-WhitneyのU値を返す
  { return v_u; }
  double uv() const             // Mann-WhitneyのU値の分散を返す
  { return v_uv; }
  double nlim() const           // 理論的境界値を返す
  { return v_nlim; }
  double nsn() const            // 理論的境界値による感度を返す
  { return v_nsn; }
  double nsp() const            // 理論的境界値による特異度を返す
  { return v_nsp; }
  double nac() const            // 理論的境界値による正診率を返す 
  { return v_nac; }
  double nppv() const           // 理論的境界値によるPPVを返す 
  { return v_nppv; }
  double nnpv() const           // 理論境界値によるNPVを返す 
  { return v_nnpv; }
  double nauc() const           // 理論的AUCを返す
  { return v_nauc; }
  double d2() const             // マハラノビスの汎距離を返す
  { return v_d2; }
  double ppr() const            // 誤判別確率を返す
  { return v_ppr; }
  statv<double>& dstat(size_t idx) // 診断項目の基礎統計量を返す
  { return statvec[idx]; }
  void cal_diag(ldata& datalst); // ROC分析と単変量判別分析で診断率を計算する
  int cal_wtest(double& dev, double& vd, double& t, int& df);
  // Welchの検定結果を返す
};

const diag& diag::operator=(const diag& diag0) // コピー代入演算子
{
  if (&diag0 == this) return *this; // 自分自身への代入を防ぐ
  v_co = diag0.v_co;  v_colu = diag0.v_colu;
  v_pi = diag0.v_pi;  v_lim = diag0.v_lim;
  v_nt = diag0.v_nt;  v_tp = diag0.v_tp;  v_tn = diag0.v_tn;
  v_fp = diag0.v_fp;  v_fn = diag0.v_fn;
  v_auc = diag0.v_auc;  v_aucv = diag0.v_aucv;
  v_u = diag0.v_u;  v_uv = diag0.v_uv;
  v_nlim = diag0.v_nlim;  v_nsn = diag0.v_nsn;
  v_nsp = diag0.v_nsp;  v_nac = diag0.v_nac;
  v_nppv = diag0.v_nppv;  v_nnpv = diag0.v_nnpv;
  v_nac = diag0.v_nac;  v_nauc = diag0.v_nauc;
  v_d2 = diag0.v_d2;  v_ppr = diag0.v_ppr;
  statvec.resize(diag0.statvec.size());
  statvec[0] = diag0.statvec[0];  statvec[1] = diag0.statvec[1];
  return *this;
}

double diag::plr() const        // 陽性尤度比を返す
{
  if (sgn(1.0 - sp()) <= 0) return (sn() > 0.0) ? BLANKf : 0.0; // BLANKf:∞
  return sn() / (1.0 - sp());
}

double diag::nlr() const        // 陰性尤度比を返す
{
  if (sgn(sp()) <= 0) return (sn() < 1.0) ? BLANKf : 0.0;
  return (1.0 - sn()) / sp();
}

double diag::ppv()              // 陽性予測値を返す
{
  if (sgn(pi()) <= 0) return BLANKf;
  if (sgn(sn()) <= 0) return (sp() < 1.0) ? 0.0 : 1.0;
  return 1.0 / (1.0 + (1.0 - pi()) / pi() * (1.0 - sp()) / sn());
}

double diag::npv()              // 陰性予測値を返す
{
  if (sgn(1.0 - pi()) <= 0) return BLANKf;
  if (sgn(sp()) <= 0) return (sn() < 1.0) ? 0.0 : 1.0;
  return 1.0 / (1.0 + pi() / (1.0 - pi()) * (1.0 - sn()) / sp());
}

void diag::cal_roc(ldata& datalst) // ROC分析を行う
{
  int n1 = statvec[0].n(), n2 = statvec[1].n();
  int n = n1 + n2, tp = n1, dfn = 0, tn = 0, dtn = 0;
  int ac_max = 0, tp_max = 0, tn_max = n2;
  double csqr_max = 0.0, yi_max = 0.0;
  double dist_min = 2.0, m1 = statvec[0].mean(), m2 = statvec[1].mean();
  cotype co = v_co;
  v_colu = lowlim;
  if (v_co == uplim)                     // 上限境界値指定の場合
  { datalst.reverse();  v_colu = uplim; } // データリストを逆順にする
  else if (v_co != lowlim)                // 上下限自動判定の場合
  {
    co = v_colu = lowlim;
    if (m1 < m2) { co = v_colu = uplim;  datalst.reverse(); }
  }
  double w = 0.0, q1 = 0.0, q2 = 0.0, k = 0.0;
  double fpre = (datalst.begin())->fdata(1), lim = v_lim;
  bool firstf = true;
  for (ldata_itr ditr = datalst.begin(); ditr != datalst.end(); ++ditr)
  {
    double f = ditr->fdata(1);
    if (f != fpre || firstf)
    {
      tp -= dfn;  tn += dtn;
      w += (double)(tp * dtn) + 0.5 * (double)(dfn * dtn);
      q1 += (double)dtn * ((double)(tp * (tp + dfn)) + (double)(dfn * dfn) / 3.0);
      q2 += (double)dfn * ((double)(tn * (tn - dtn)) + (double)(dtn * dtn) / 3.0);
      int g = dfn + dtn;
      k += (double)g * (double)(g * g - 1);
      int ac = tp + tn;
      int fn = n1 - tp, fp = n2 - tn;
      double csqr = 0.0;
      double yi = 0.0, x = 0.0, y = 0.0, dist = 0.0;
      fpre = f;  dfn = dtn = 0;
      switch (v_co)
      {
      case lowlim :             // 下限境界値指定
        if (tp_max == 0 && lim <= f) { tp_max = tp;  tn_max = tn; }
        break;
      case uplim :              // 上限境界値指定
        if (tp_max == 0 && f <= lim) { tp_max = tp;  tn_max = tn; }
        break;
      case colim :              // 境界値指定(上下限は自動判定)
        if (co == lowlim)
        {
          if (tp_max == 0 && lim <= f) { tp_max = tp;  tn_max = tn; }
        }
        else
        {
          if (tp_max == 0 && f <= lim) { tp_max = tp;  tn_max = tn; }
        }
        break;
      case coac :             // 正診率基準
        if (ac > ac_max)
        { ac_max = ac;  tp_max = tp;  tn_max = tn;  lim = f; }
        break;
      case cocm : csqr = csqr22(tn, fp, fn, tp);
        if (csqr > csqr_max)    // 最大χ2乗法
        { csqr_max = csqr;  tp_max = tp;  tn_max = tn;  lim = f; }
        break;
      case coyi : yi = (double)tp / (double)n1 - (double)fp / (double)n2;
        // Youden index
        if (yi > yi_max)
        { yi_max = yi;  ac_max = ac;  tp_max = tp;  tn_max = tn;  lim = f; }
        break;
      case coroc : x = (double)fn / (double)n1; // ROC曲線最適点基準
        y = (double)fp / (double)n2;
        dist = x * x + y * y;
        if (dist < dist_min)
        {
          dist_min = dist;
          ac_max = ac;  tp_max = tp;  tn_max = tn;  lim = f;
        }
        break;
      }
      firstf = false;
    }
    if ((int)(ditr->fdata(0)) == 0) ++dfn;
    else ++dtn;
  }
  v_co = co;  v_lim = lim;
  v_tp = tp_max;  v_fn = n1 - tp_max;  v_fp = n2 - tn_max;  v_tn = tn_max;
  w += 0.5 * (double)(dfn * dtn);
  q1 += (double)dtn * ((double)(dfn * dfn) / 3.0);
  q2 += (double)dfn * ((double)(n2 * tn) + (double)(dtn * dtn) / 3.0);
  k = 1.0 - k / (double)(n - 1) / (double)n / (double)(n + 1);
  v_u = w;  v_uv = k * (double)n1 * (double)n2 * (double)(n + 1) / 12.0;
  w = w / (double)n1 / (double)n2;
  q1 = q1 / (double)n1 / (double)n1 / (double)n2;
  q2 = q2 / (double)n1 / (double)n2 / (double)n2;
  v_auc = w;
  v_aucv = (w * (1.0 - w) + (n1 - 1) * (q1 - w * w) + (n2 - 1) * (q2 - w * w)) / (double)n1 / (double)n2;
  const double Veps = FDeps * 5.0;
  if (v_aucv < Veps) v_aucv = 0.0;
}

void diag::cal_dis()            // 単変量判別分析を行う
{
  int n1 = statvec[0].n(), n2 = statvec[1].n();
  double m1 = statvec[0].mean(), m2 = statvec[1].mean();
  double sd1 = statvec[0].sd(), sd2 = statvec[1].sd();
  double v1 = statvec[0].vr(), v2 = statvec[1].vr();
  v_d2 = (v1 + v2 > 0.0) ? 2.0 * (m1 - m2) * (m1 - m2) / (v1 + v2) : BLANKf;
  v_nlim = (m1 + m2) * 0.5;     // 判別分析の境界値
  if (fabs(v_nlim) < FDeps) v_nlim = 0.0;
  double dev = m1 - v_nlim;
  if (v_colu == uplim) dev = -dev;
  v_nsn = (sd1 > 0.0) ? sndis(dev / sd1) : 1.0;
  dev = v_nlim - m2;
  if (v_colu == uplim) dev = -dev;
  v_nsp = (sd2 > 0.0) ? sndis(dev / sd2) : 1.0;
  v_nac = (n1 + n2 > 0) ? ((double)n1 * v_nsn + (double)n2 * v_nsp) / (double)(n1 + n2) : 1.0;
  double orpi = (1.0 - pi()) / pi();
  v_nppv = (v_nsn > 0.0) ? 1.0 / (1.0 + orpi * (1.0 - v_nsp) / v_nsn) : 0.0;
  v_nnpv = (v_nsp > 0.0) ? 1.0 / (1.0 + (1.0 - v_nsn) / v_nsp / orpi) : 0.0;
  v_nauc = 1.0;  v_ppr = 0.0;
  if (v_d2 == BLANKf) return;
  double z = sqrt(v_d2 * 0.5);
  v_nauc = sndis(z);            // 理論AUC
  z = sqrt(v_d2) * 0.5;
  v_ppr = 1.0 - (double)(sndis(z));
}

void diag::cal_diag(ldata& datalst)
// ROC分析と単変量判別分析で診断率を計算する
{
  for (auto& cdatai : datalst) statvec[(size_t)(cdatai.fdata(0))].sigma(cdatai.fdata(1));
  v_nt = statvec[0].n() + statvec[1].n();
  if (v_nt < 1 || statvec[0].n() < 1 || statvec[1].n() < 1) return; // 計算できず
  if (v_pi == BLANKf) v_pi = (double)(statvec[0].n()) / (double)v_nt; // -pi指定なしの時
  cal_roc(datalst);             // ROC分析
  cal_dis();                    // 単変量判別分析
}

int diag::cal_wtest(double& dev, double& vd, double& t, int& df)
// Welchの検定結果を返す
{
  dev = fabs(statvec[0].mean() - statvec[1].mean());
  int n1 = statvec[0].n(), n2 = statvec[1].n();
  double v1 = statvec[0].vr() / (double)n1;
  double v2 = statvec[1].vr() / (double)n2;
  vd = v1 + v2;
  df = (vd > 0.0) ? (int)(vd * vd / (v1 * v1 / (n1 - 1) + v2 * v2 / (n2 -1))) : n1 + n2  - 2;
  t = (vd > 0.0) ? dev / sqrt(vd) : 0.0;
  return df;                    // 自由度を返す
}

class rocclass                  // rocメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout); // 分析結果のタイトルとデータ識別情報を出力する
  void put_result(ostream& fout, option& opt, diag& dg);
  // 診断率の計算結果を出力する
  void put_stat(ostream& fout, option& opt, diag& dg) const;
  // 基礎統計量を出力する
  void put_roca(ostream& fout, option& opt, diag& dg); // ROC分析の結果を出力する
  void put_pnpv(ostream& fout, option& opt, diag& dg);
  // 陽性予測値と陰性予測値を出力する
  void put_binci(ostream& fout, int n, int nt, double cc) const;
  // F分布を利用した出現率の信頼区間を出力する
  void put_lrci(ostream& fout, double& lrl, double& lru, int n1, int nt1, int n2, int nt2, double cc) const;
  // 尤度比の信頼区間を出力し、値を返す
  void put_pvci(ostream& fout, double& pvl, double& pvu, bool ppv, double p0, double lrl, double lru, double cc) const;
  // 陽性/陰性予測値の信頼区間を出力し、値を返す
  void put_auc(ostream& fout, option& opt, diag& dg);
  // ROC曲線のAUC(曲線下面積)とMann-WhitneyのU値と検定結果を出力する
  void put_ozr(ostream& fout, option& opt, diag& dg);
  // オッズ比(True群/False群)とその検定結果を出力する
  void put_uda(ostream& fout, option& opt, diag& dg);
  // 単変量判別分析結果と理論的AUCを出力する
  void put_dns(option& opt, diag& dg);
  // データファイル形式の診断率を出力する
  void put_dnsid(ofstream& fout);
  // DANSデータのタイトルとデータ識別情報を出力する
  void plot_ngraph(option& opt, diag& dg); // Ngraph用データを出力する
  void plot_ngpid(ofstream& fout);
  // Ngraph用データのタイトルとデータ識別情報を出力する
  void plot_roc(ofstream& fout, diag& dg); // Ngraph用ROC曲線を出力する
  void plot_roci(ofstream& fout, option& opt, diag& dg);
  // Ngraph用ROC曲線信頼区間を出力する
  void plot_nroc(ofstream& fout, diag& dg) const;
  // Ngraph用理論ROC曲線を出力する
  void plot_nroci(ofstream& fout, option& opt, diag& dg) const;
  // Ngraph用理論ROC曲線信頼区間を出力する
public:
  rocclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~rocclass()                   // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み累積する
  void put(option& opt);                     // 集計結果を出力する
};

void rocclass::put_resid(ostream& fout)
// 分析結果のタイトルとデータ識別情報を出力する
{
  puttitle(fout, Title);
  fout << "\n";
  id.put_title(fout);
  id.put_stra(fout);
  fout << "群項目  :" << id.item(0) << "\n";
  fout << "集計項目:" << id.item(1) << "\n";
}

void rocclass::put_result(ostream& fout, option& opt, diag& dg)
// 診断率の計算結果を出力する
{
  constexpr int LENline = 74;   // 区切り線の長さ
  put_resid(fout);
  bool putres = false;
  if (opt.res("stat"))
  {
    fout << "\n";
    putline(fout, LENline);
    put_stat(fout, opt, dg);
    putres = true;
  }
  if (dg.n1() * dg.n2() < 1)
  {
    putline(fout, LENline);
    fout << "※ROC分析不可能！\n\n" << endl;
    return;
  }
  if (opt.res("roc"))
  {
    if (putres == false) fout << "\n";
    putline(fout, LENline);
    put_roca(fout, opt, dg);
    putres = true;
  }
  if (opt.res("auc"))
  {
    if (putres == false) fout << "\n";
    putline(fout, LENline);
    put_auc(fout, opt, dg);
    putres = true;
  }
  if (opt.res("or"))
  {
    if (putres == false) fout << "\n";
    putline(fout, LENline);
    put_ozr(fout, opt, dg);
    putres = true;
  }
  if (opt.res("uda"))
  {
    if (putres == false) fout << "\n";
    putline(fout, LENline);
    put_uda(fout, opt, dg);
    putres = true;
  }
  if (putres) putline(fout, LENline);
  fout << "\n" << endl;
}

void rocclass::put_stat(ostream& fout, option& opt, diag& dg) const
// 基礎統計量を出力する
{
  fout << "群1(疾患群):コード=";
  if ((opt.g(0)).empty()) fout << opt.g(1) << " 以外  ";
  else fout << opt.g(0) << "       ";
  fout << "True-Positive=" << dg.tp() << "例 False-Negative=" << dg.fn() << "例\n";
  fout << "  ";
  putmsd(fout, (dg.dstat(0)).n(), (dg.dstat(0)).mean(), (dg.dstat(0)).vr());
  fout << "群2(正常群):コード=";
  if ((opt.g(1)).empty()) fout << opt.g(0) << " 以外  ";
  else fout << opt.g(1) << "       ";
  fout << "False-Positive=" << dg.fp() << "例 True-Negative=" << dg.tn() << "例\n";
  fout << "  ";
  putmsd(fout, (dg.dstat(1)).n(), (dg.dstat(1)).mean(), (dg.dstat(1)).vr());
}

void rocclass::put_roca(ostream& fout, option& opt, diag& dg)
// ROC分析の結果を出力する
{
  string tstr("以上");
  if (dg.colu() == uplim) tstr = "以下";
  fout << "ROC分析:境界値=" << dg.lim() << " " << tstr;
  if (opt.co() == coac) tstr = "(正診率最大点)";
  else if (opt.co() == cocm) tstr = "(最大χ2乗法)";
  else if (opt.co() == coyi) tstr = "(Youden index)";
  else if (opt.co() == coroc) tstr = "(ROC曲線最適点)";
  else if (opt.co() == colim) tstr = "(指定値)";
  else tstr.erase();
  fout << "を群1と診断した場合" << tstr << "\n";
  int nt0 = dg.fn() + dg.tp();
  fout << "  感度(sensitivity)=" << dg.sn() << "(" << dg.tp() <<"/" << nt0 << ")";
  put_binci(fout, dg.tp(), nt0, opt.cc());
  fout << "\n";
  int nt1 = dg.fp() + dg.tn();
  fout << "  特異度(specificity)=" << dg.sp() << "(" << dg.tn() << "/" << nt1 << ")";
  put_binci(fout, dg.tn(), nt1, opt.cc());
  fout << "\n";
  int nt = nt0 + nt1;
  fout << "  正診率(accuracy)=" << dg.ac() << "(" << dg.tp() + dg.tn() <<"/" << nt << ")";
  put_binci(fout, dg.tp() + dg.tn(), nt, opt.cc());
  fout << "\n";
  put_pnpv(fout, opt, dg);
}

void rocclass::put_pnpv(ostream& fout, option& opt, diag& dg)
// 陽性予測値と陰性予測値を出力する
{
  int nt0 = dg.fn() + dg.tp();
  int nt1 = dg.fp() + dg.tn();
  double lrpl = 1.0, lrpu = 1.0, lrnl = 1.0, lrnu = 1.0;
  fout << "  陽性尤度比(LR+)=";
  if (dg.plr() == BLANKf)
  { fout << "∞\n"; }
  else
  {
    fout << dg.plr();
    if (dg.plr() > 0.0) put_lrci(fout, lrpl, lrpu, dg.tp(), nt0, dg.fp(), nt1, opt.cc());
  }
  fout << "  陰性尤度比(LR-)=";
  if (dg.nlr() == BLANKf)
  { fout << "∞\n"; }
  else
  {
    fout << dg.nlr();
    if (dg.nlr() > 0.0) put_lrci(fout, lrnl, lrnu, dg.fn(), nt0, dg.tn(), nt1, opt.cc());
  }
  fout << "群1の事前確率=" << dg.pi() << "  群2の事前確率=" << 1.0 - dg.pi() << "の時\n";
  double ppvl = 1.0, ppvu = 1.0, npvl = 1.0, npvu = 1.0;
  fout << "  陽性予測値(PPV)=" << dg.ppv();
  if (0.0 < dg.ppv() && dg.ppv() < 1.0) put_pvci(fout, ppvl, ppvu, true, dg.pi(), lrpl, lrpu, opt.cc());
  fout << "  陰性予測値(NPV)=" << dg.npv();
  if (0.0 < dg.npv() && dg.npv() < 1.0) put_pvci(fout, npvl, npvu, false, dg.pi(), lrnl, lrnu, opt.cc());
  fout << "  リスク比(RR)=";
  if (dg.npv() < 1.0)
  {
    double rr = dg.ppv() / (1.0 - dg.npv());
    fout << rr;
    double rrl = (npvl < 1.0) ? ppvl / (1.0 - npvl) : 1.0;
    double rru = (npvu < 1.0) ? ppvu / (1.0 - npvu) : 1.0;
    fout << " (" << opt.cc() * 100.0 << "%CI:" << rrl << "-" << rru << ")";
  }
  else
  { fout << "∞"; }
  fout << endl;
}

void rocclass::put_binci(ostream& fout, int n, int nt, double cc) const
// F分布を利用した出現率の信頼区間を出力する
{
  double plow = 0.0, pup = 1.0;
  calbexci(plow, pup, n, nt, cc);
  fout << " (" << cc * 100.0 << "%CI:" << plow << "-" << pup << ")";
}

void rocclass::put_lrci(ostream& fout, double& lrl, double& lru, int n1, int nt1, int n2, int nt2, double cc) const
// 尤度比の信頼区間を出力し、値を返す
{
  double lr = 1.0, b = 0.0, v = 0.0, ci = 0.0, bl = 0.0, bu = 0.0;
  lrl = 1.0;  lru = 1.0;
  if (n1 > 0 && nt1 > 0 && n2 > 0 && nt2 > 0)
  {
    lr = (double)n1 / (double)nt1 * (double)nt2 / (double)n2;
    b = log(lr);
    v = 1.0 / (double)n1 + 1.0 / (double)n2 - 1.0 / (double)nt1 - 1.0 / (double)nt2;
    if (v < FDeps) v = 0.0;
    ci = ptoz(1.0 - cc) * sqrt(v);
    bl = b - ci;  bu = b + ci;
    lrl = exp(bl);  lru = exp(bu);
  }
  fout << " (" << cc * 100.0 << "%CI:" << lrl << "-" << lru << ")\n";
}

void rocclass::put_pvci(ostream& fout, double& pvl, double& pvu, bool ppv, double p0, double lrl, double lru, double cc) const
// 陽性/陰性予測値の信頼区間を出力し、値を返す
{
  if (p0 <= 0.0 || 1.0 <= p0) return;
  pvl = pvu = 1.0;
  if (lrl > 0.0 && lru > 0.0)
  {
    if (ppv)                    // 陽性予測値
    {
      pvl = 1.0 / (1.0 + (1.0 - p0) / p0 / lrl);
      pvu = 1.0 / (1.0 + (1.0 - p0) / p0 / lru);
    }
    else                        // 陰性予測値
    {
      pvl = 1.0 / (1.0 + p0 / (1.0 - p0) * lru);
      pvu = 1.0 / (1.0 + p0 / (1.0 - p0) * lrl);
    }
  }
  if (pvl < 0.0) pvl = 0.0;
  if (pvu > 1.0) pvu = 1.0;
  fout << " (" << cc * 100.0 << "%CI:" << pvl << "-" << pvu << ")\n";
}

void rocclass::put_auc(ostream& fout, option& opt, diag& dg)
// ROC曲線のAUC(曲線下面積)とMann-WhitneyのU値と検定結果を出力する
{
  double se = sqrt(dg.aucv()), z = 0.0, ci = 0.0;
  fout << "ROC曲線のAUC=" << dg.auc();
  z = ptoz(1.0 - opt.cc());
  ci = z * se;
  double aucl = dg.auc() - ci, aucu = dg.auc() + ci;
  if (aucl < 0.0) aucl = 0.0;
  if (aucu > 1.0) aucu = 1.0;
  fout << " (" << opt.cc() * 100.0 << "%CI:" << aucl << "-" << aucu << ")";
  fout << "  標準誤差=" << se << "\n";
  int n1 = (dg.dstat(0)).n(), n2 = (dg.dstat(1)).n();
  double um = (double)n1 * (double)n2 * 0.5;
  se = sqrt(dg.uv());
  fout << "Mann-WhitneyのU値=" << dg.u() << "=AUC×" << n1 << "×" << n2;
  ci = z * se;
  double ul = dg.u() - ci, uu = dg.u() + ci;
  if (ul < 0.0) ul = 0.0;
  if (uu > n1 * n2) uu = double(n1) * double(n2);
  fout << " (" << opt.cc() * 100.0 << "%CI:" << ul << "-" << uu << ")";
  z = (sgn(se) > 0) ? (dg.u() - um) / se : 0.0;
  ldbl pr = ztop(z);
  fout << "\n  正規分布:z=" << z << "  有意確率:p=" << pr << pind(pr) << endl;
}

void rocclass::put_ozr(ostream& fout, option& opt, diag& dg)
// オッズ比(群1/群2)とその検定結果を出力する
{
  double n11 = dg.tp(), n12 = dg.fn(), n21 = dg.fp(), n22 = dg.tn();
  if (dg.tp() == 0 || dg.fn() == 0 || dg.fp() == 0 || dg.tn() == 0)
  { n11 += 0.5;  n12 += 0.5;  n21 += 0.5;  n22 += 0.5; }
  double ozr = n11 / n12 * n22 / n21;
  double b = log(ozr), v = 1.0 / n11 + 1.0 / n12 + 1.0 / n21 + 1.0 / n22;
  double csqr = b * b / v, se = sqrt(v), ci = 0.0;
  fout << "オッズ比(群1/群2):OR=" << ozr;
  ci = ptoz(1.0 - opt.cc()) * se;
  fout << " (" << opt.cc() * 100.0 << "%CI:" << exp(b - ci) << "-" << exp(b + ci) << ")\n";
  fout << "  ln(OR)=" << b;
  fout << " (" << opt.cc() * 100.0 << "%CI:" << b - ci << "-" << b + ci << ")";
  ldbl pr = ctop(csqr, 1);
  fout << "\n  χ2乗=" << csqr << "  自由度=1  有意確率:p=" << pr << pind(pr) << endl;
}

void rocclass::put_uda(ostream& fout, option& opt, diag& dg)
// 単変量判別分析結果と理論的AUCを出力する
{
  string tstr("以上");
  if (dg.colu() == uplim) tstr = "以下";
  fout << "単変量線型判別分析:理論的境界値=" << dg.nlim() << " " << tstr;
  fout << "を群1と診断した場合\n";
  fout << "  理論的感度=" << dg.nsn() << "  理論的特異度=" << dg.nsp() << "  理論的正診率=" << dg.nac() << "\n";
  if (0.0 < opt.pi()) fout << "  理論的PPV=" << dg.nppv() << "  理論的NPV=" << dg.nnpv() << "\n";
  if (dg.d2() == BLANKf) return;
  double n1 = (double)((dg.dstat(0)).n()), n2 = (double)((dg.dstat(1)).n());
  double f = dg.d2() * n1 * n2 / (n1 + n2);
  double r2 = f / (f + n1 + n2 - 2.0);
  fout << "判別効率(Mahalanobisの平方距離):D2乗=" << dg.d2() << "  相関比:η=" << sqrt(r2) << "\n";
  fout << "  判別確率:Φ(D/2)=Φ(" << sqrt(dg.d2()) * 0.5 << ")=" << 1.0 - dg.ppr() << "(誤判別確率:1-p=" << dg.ppr() << ")\n";
  fout << "理論的ROC曲線のAUC=Φ(D/√2)=Φ(" << sqrt(dg.d2() * 0.5) << ")=" << dg.nauc() << "\n";
  double dev = 0.0, vd = 0.0, t = 0.0;
  int df = 0;
  dg.cal_wtest(dev, vd, t, df);
  double pr = ttop(t, df);
  fout << "Welchの検定:t=" << t << "  自由度=" << df << "  有意確率:p=" << pr << pind(pr) << "\n";
  double ci = calmcci(mcfis, opt.cc(), 2, 1, df, 1, vd);
  fout << "  平均値の差の" << opt.cc() * 100.0 << "%CI=";
  fout << dev << "±" << ci << "(" << dev - ci << "-" << dev + ci << ")\n";
}

void rocclass::put_dns(option& opt, diag& dg) // データファイル形式の診断率を出力する
{
  ofstream fout;
  if (ofopen(fout, (opt.dnsname()).c_str()) == false) return; // ファイルオープン失敗
  put_dnsid(fout);
  datac cdata(3);
  int n1 = (dg.dstat(0)).n(), n2 = (dg.dstat(1)).n();
  int tp = n1, dfn = 0, tn = 0, dtn = 0;
  double fpre = (datalst.begin())->fdata(1);
  for (auto& cdatai : datalst)
  {
    double f = cdatai.fdata(1);
    if (f != fpre)
    {
      tp -= dfn;  tn += dtn;
      fpre = f;  dfn = dtn = 0;
      double sn = (double)tp / (double)n1, sp = (double)tn / (double)n2;
      cdata.assign(0, f);  cdata.assign(1, sn);  cdata.assign(2, sp);
      cdata.put(fout);
    }
    if ((int)(cdatai.fdata(0)) == 0) ++dfn;
    else ++dtn;
  }
}

void rocclass::put_dnsid(ofstream& fout)
// DANSデータのタイトルとデータ識別情報を出力する
{
  id.put_title(fout, formtag);
  id.put_stra(fout, formtag);
  idd tid(id.title(), 3);
  tid.assign_item(0, "境界値");
  tid.assign_item(1, "感度");
  tid.assign_item(2, "特異度");
  tid.put_item(fout, formtag);
}

void rocclass::plot_ngraph(option& opt, diag& dg) // Ngraph用データを出力する
{
  ofstream fout;
  if (ofopen(fout, (opt.ngpname()).c_str()) == false) return; // ファイルオープン失敗
  bool delml = false;
  if (opt.ngp("roc"))
  {
    fout << "-s1 -dline -w80 -zx1,0,0.2 -zy0,1,0.2\n";
    plot_ngpid(fout);
    plot_roc(fout, dg);
    fout << NGPdelm << "\n";
    plot_nroc(fout, dg);
    delml = true;
  }
  if (opt.ngp("roci"))
  {
    if (delml) fout << NGPdelm << "\n";
    else fout << "-s1 -dline -w80 -zx1,0,0.2 -zy0,1,0.2\n";
    plot_ngpid(fout);
    plot_roci(fout, opt, dg);
    fout << NGPdelm << "\n";
    plot_nroci(fout, opt, dg);
  }
}

void rocclass::plot_ngpid(ofstream& fout)
// Ngraph用データのタイトルとデータ識別情報を出力する
{
  fout << "%";  puttitle(fout, Title);
  fout << "%\n";
  id.put_title(fout, formngp);
  id.put_stra(fout, formngp);
  fout << "%群項目  :" << id.item(0) << "\n";
  fout << "%集計項目:" << id.item(1) << endl;
}

void rocclass::plot_roc(ofstream& fout, diag& dg) // Ngraph用ROC曲線を出力する
{
  int n1 = (dg.dstat(0)).n(), n2 = (dg.dstat(1)).n();
  int tp = n1, dfn = 0, tn = 0, dtn = 0;
  fout << "%ROC曲線\n";
  fout << "0 1\n";
  double fpre = (datalst.begin())->fdata(1);
  double f = 0.0, sn = 0.0, sp = 0.0;
  for (auto& cdatai : datalst)
  {
    f = cdatai.fdata(1);
    if (f != fpre)
    {
      tp -= dfn;  tn += dtn;
      fpre = f;  dfn = dtn = 0;
      sn = (double)tp / (double)n1;  sp = (double)tn / (double)n2;
      fout << sp << " " << sn << "\n";
    }
    if ((int)(cdatai.fdata(0)) == 0) ++dfn;
    else ++dtn;
  }
  fout << "1 0\n";
}

void rocclass::plot_roci(ofstream& fout, option& opt, diag& dg)
// Ngraph用ROC曲線の信頼区間を出力する
{
  int n1 = (dg.dstat(0)).n(), n2 = (dg.dstat(1)).n();
  int tp = n1, dfn = 0, tn = 0, dtn = 0;
  double fpre = (datalst.begin())->fdata(1);
  double f = 0.0, sn = 0.0, sp = 0.0, plow = 0.0, pup = 1.0;
  fout << "%ROC曲線の" << opt.cc() * 100.0 << "%CI下限\n"; // 信頼区間下限
  fout << "0 1\n";
  for (auto& cdatai : datalst)
  {
    f = cdatai.fdata(1);
    if (f != fpre)
    {
      tp -= dfn;  tn += dtn;
      fpre = f;  dfn = dtn = 0;
      sp = (double)tn / (double)n2;
      calbexci(plow, pup, tp, n1, opt.cc());
      fout << sp << " " << plow << "\n";
    }
    if ((int)(cdatai.fdata(0)) == 0) ++dfn;
    else ++dtn;
  }
  fout << "1 0\n";
  tp = n1;  dfn = 0;  tn = 0;  dtn = 0;
  fpre = (datalst.begin())->fdata(1);
  fout << NGPdelm2 << "\n";
  fout << "%ROC曲線の" << opt.cc() * 100.0 << "%CI上限\n"; // 信頼区間上限
  fout << "0 1\n";
  for (auto& cdatai : datalst)
  {
    f = cdatai.fdata(1);
    if (f != fpre)
    {
      tp -= dfn;  tn += dtn;
      fpre = f;  dfn = dtn = 0;
      sp = (double)tn / (double)n2;
      calbexci(plow, pup, tp, n1, opt.cc());
      fout << sp << " " << pup << "\n";
    }
    if (int(cdatai.fdata(0)) == 0) ++dfn;
    else ++dtn;
  }
  fout << "1 0\n";
}

void rocclass::plot_nroc(ofstream& fout, diag& dg) const
// Ngraph用理論的ROC曲線を出力する
{
  double m1 = (dg.dstat(0)).mean(), sd1 = (dg.dstat(0)).sd();
  double m2 = (dg.dstat(1)).mean(), sd2 = (dg.dstat(1)).sd();
  fout << "%理論的ROC曲線\n";
  if (sgndiff(m1, m2) == 0 && sgn(sd1 + sd2) == 0) // データが定数の時
  {
    fout << "0 1\n";
    fout << "1 0\n";
    return;
  }
  double xmin = m2 - sd2 * 4.0, xmax = m1 + sd1 * 4.0;
  if (m1 < m2)                  // 群1の平均値の方が小さい時
  { xmin = m1 - sd1 * 4.0;  xmax = m2 + sd2 * 4.0; }
  double dx = (xmax - xmin) / (double)RDdef;
  double sp = 0.0, sn = 0.0;
  if (dg.colu() == lowlim) fout << "0 1\n";
  else fout << "1 0\n";
  double x = xmin;
  for (int i = 0; i <= RDdef; ++i, x += dx)
  {
    if (dg.colu() == lowlim)
    { sp = sndis((x - m2) / sd2);  sn = 1.0 - sndis((x - m1) / sd1); }
    else
    { sp = 1.0 - sndis((x - m2) / sd2);  sn = sndis((x - m1) / sd1); }
    fout << sp << " " << sn << "\n";
  }
  if (dg.colu() == lowlim) fout << "1 0\n";
  else fout << "0 1\n";
}

void rocclass::plot_nroci(ofstream& fout, option& opt, diag& dg) const
// Ngraph用理論的ROC曲線の信頼区間を出力する
{
  int n1 = (dg.dstat(0)).n(), n2 = (dg.dstat(1)).n();
  double m1 = (dg.dstat(0)).mean(), sd1 = (dg.dstat(0)).sd();
  double m2 = (dg.dstat(1)).mean(), sd2 = (dg.dstat(1)).sd();
  if (sgndiff(m1, m2) == 0 && sgn(sd1 + sd2) == 0) // データが定数の時
  {
    fout << "%理論的ROC曲線の" << opt.cc() * 100.0 << "%CI下限\n";
    fout << "0 1\n";
    fout << "1 0\n";
    fout << "%理論的ROC曲線の" << opt.cc() * 100.0 << "%CI上限\n";
    fout << "0 1\n";
    fout << "1 0\n";
    return;
  }
  double xmin = m2 - sd2 * 4.0, xmax = m1 + sd1 * 4.0;
  if (m1 < m2)                  // 群1の平均値の方が小さい時
  { xmin = m1 - sd1 * 4.0;  xmax = m2 + sd2 * 4.0; }
  double dx = (xmax - xmin) / (double)RDdef;
  double sp = 0.0, sn = 0.0, plow = 0.0, pup = 1.0;
  fout << "%理論的ROC曲線の" << opt.cc() * 100.0 << "%CI下限\n"; // 信頼区間下限
  if (dg.colu() == lowlim) fout << "0 1\n";
  else fout << "1 0\n";
  double x = xmin;
  for (int i = 0; i <= RDdef; ++i, x += dx)
  {
    if (dg.colu() == lowlim)
    { sp = sndis((x - m2) / sd2);  sn = 1.0 - sndis((x - m1) / sd1); }
    else
    { sp = 1.0 - sndis((x - m2) / sd2);  sn = sndis((x - m1) / sd1); }
    calbnmci(plow, pup, sn, n1, opt.cc());
    fout << sp << " " << plow << "\n";
  }
  if (dg.colu() == lowlim) fout << "1 0\n";
  else fout << "0 1\n";
  fout << NGPdelm2 << "\n";
  fout << "%理論的ROC曲線の" << opt.cc() * 100.0 << "%CI上限\n"; // 信頼区間上限
  if (dg.colu() == lowlim) fout << "0 1\n";
  else fout << "1 0\n";
  x = xmin;
  for (int i = 0; i <= RDdef; ++i, x += dx)
  {
    if (dg.colu() == lowlim)
    { sp = sndis((x - m2) / sd2);  sn = 1.0 - sndis((x - m1) / sd1); }
    else
    { sp = 1.0 - sndis((x - m2) / sd2);  sn = sndis((x - m1) / sd1); }
    calbnmci(plow, pup, sn, n1, opt.cc());
    fout << sp << " " << pup << "\n";
  }
  if (dg.colu() == lowlim) fout << "1 0\n";
  else fout << "0 1\n";
}

int rocclass::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("実行時オプションエラー！");
  nc = (int)(id.nc());          // ケース数変更用
  ldata_itr ditr = datalst.begin();
  while (ditr != datalst.end())
  {
    int idx = opt.search_g(ditr->sdata(0));
    double fdata = ditr->fdata(1);
    if (idx < 0 || fdata == BLANKf) // 欠損データはリストから削除する
    { ditr = datalst.erase(ditr);  --nc; }
    else
    {
      ditr->assign(0, (double)idx); // 群コードを群インデックスに変える
      ++ditr;
    }
  }
  if (nc < (int)(id.nc())) id.nc(nc); // ケース数を減らす
  return nc;
}

void rocclass::put(option& opt) // 集計結果を出力する
{
  datalst.sort(data_less(1));    // 項目データ順にソートする
  diag dg(opt.co(), opt.pi(), opt.lim());
  dg.cal_diag(datalst);
  if (opt.res("res")) put_result(cout, opt, dg);
  if (dg.n1() * dg.n2() < 1) return;
  if (opt.dnsname() != "") put_dns(opt, dg);
  if (opt.ngpname() != "") plot_ngraph(opt, dg);
}

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