// plot.C Version 7.42  Copyleft (c) 杉本典夫
// Time-stamp: <2025-11-03 14:20:21 tac>
//
// Ngraph用プロットデータの出力。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -W plot plot.cpp -I/usr/local/include
//              -L/usr/local/lib -ldans -lfun

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

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

const vecstr USAGE =            // プログラム使用法
{
  "plot  Version 7.42           Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:plot [<オプション> …]",
  "機能:Ngraph用プロットデータを出力します。",
  "<オプション>",
  "  -ct={-|+}:連続修正の有無を指定します。",
  "    -:無(オプション省略時)  +:有",
  "  -dtype={as|b|beta|c|d|dun|e|etc|exp|f|fis|ln|log|n|ncc|ncf|nct|npv|per",
  "          |plog|poi|poly|ppv|pro|sn|t|tuk|w1|w2|wb}",
  "    分布の種類を指定します。",
  "  -f=<ファイル名>:指定されたファイルからオプションを読み込みます。",
  "  -fpara=[<実数1>,…]:実数パラメータを指定します。",
  "  -ftype={den|dis|haz|pinv|pval|sur}:関数の種類を指定します。",
  "  -h[=<ページャ名>]:オンラインマニュアルを表示します。(-h:-h=man)",
  "  -ipara=[<整数1>,…]:整数パラメータを指定します。",
  "  -range=<最小値>-<最小値>:座標軸の最小値と最大値を指定します。"
};

constexpr ldbl Pmin = (ldbl)FDmin; // pr最小値
constexpr ldbl Peps = 1.0e-10L;    // pr収束誤差
constexpr ldbl Delta = 1.0e-5L;    // きざみ幅
constexpr uint RDdef = 50;         // 省略時座標軸分割数

enum distype                    // 分布の種類
{
  disnon,                       // 指定なし
  disas,                        // 逆正弦曲線
  disb,                         // 二項分布
  disbeta,                      // ベータ分布
  disc,                         // χ2乗分布
  disd,                         // Mahalanobisの汎距離分布
  disdun,                       // Dunnettのd分布
  dise,                         // 指数分布
  disexp,                       // 指数関数
  disf,                         // F分布
  disfis,                       // 超幾何分布(Fisherの正確検定)
  disln,                        // 対数正規分布
  dislog,                       // ロジット曲線
  disn,                         // 正規分布
  disncc,                       // 非心χ2乗分布
  disncf,                       // 非心F分布
  disnct,                       // 非心t分布
  disnpv,                       // NPV曲線
  disper,                       // 周期曲線
  displog,                      // ロジスティック関数の陽性率曲線
  dispoi,                       // Poisson分布
  dispoly,                      // 多項関数
  disppv,                       // PPV曲線
  dispro,                       // プロビット曲線
  dissn,                        // 標準正規分布
  dist,                         // t分布
  distuk,                       // Tukeyのq分布(Student化された範囲のq)
  disw1,                        // Wilcoxonの符号付順位分布
  disw2,                        // Wilcoxonの順位和分布
  diswb,                        // Weibull分布
  disetc                        // その他の分布
};

enum funtype                    // 関数の種類
{
  funnon,                       // 指定なし
  funden,                       // 分布密度関数
  fundis,                       // 分布関数
  funhaz,                       // ハザード関数
  funpval,                      // 確率
  funpinv,                      // 確率逆関数
  funsur,                       // 生存関数
  funetc                        // その他の関数
};

ldbl ddis(ldbl d, int g, int df, ldbl r) // Dunnettのd分布関数
{
  if (d < 0.0L) return dtop(-d, g, df, r) * 0.5L;
  ldbl pr = 1.0L - dtop(d, g, df, r) * 0.5L;
  return (sgn(pr) == 0) ? 0.0L : pr;
}

ldbl dpval(ldbl d, int g, int df, ldbl r) // Dunnettのd分布両側確率関数
{
  return dtop(d, g, df, r);
}

ldbl qdis(ldbl q, int g, int df) // Tukeyのq分布関数
{
  if (q < 0.0L) return qtop(-q, g, df) * 0.5L;
  ldbl pr = 1.0L - qtop(q, g, df) * 0.5L;
  return (sgn(pr) == 0) ? 0.0L : pr;
}

ldbl qpval(ldbl q, int g, int df) // Tukeyのq分布両側確率関数
{
  ldbl pr = qtop(q, g, df);
  return (sgn(pr) == 0) ? 0.0L : pr;
}

ldbl eden(ldbl t, ldbl lamda)   // 指数分布密度関数(死亡率関数)
{
  ldbl pr = lamda * exp(-lamda * t);
  return (sgn(pr) == 0) ? 0.0L : pr;
}

ldbl edis(ldbl t, ldbl lamda)   // 指数分布関数(死亡関数)
{
  ldbl pr = 1.0L - exp(-lamda * t);
  return (sgn(pr) == 0) ? 0.0L : pr;
}

ldbl esur(ldbl t, ldbl lamda)   // 指数分布生存関数(片側確率関数)
{
  ldbl pr = exp(-lamda * t);
  return (sgn(pr) == 0) ? 0.0L : pr;
}

ldbl ehaz(ldbl t, ldbl lamda)   // 指数分布ハザード関数
{
  return lamda;
}

ldbl wbden(ldbl t, ldbl lamda, int a) // Weibull分布密度関数(死亡率関数)
{
  ldbl pr = pow(lamda, a) * (ldbl)a * pow(t, a - 1) * exp(-pow(lamda * t, a));
  return (sgn(pr) == 0) ? 0.0L : pr;
}

ldbl wbdis(ldbl t, ldbl lamda, int a) // Weibull分布関数(死亡関数)
{
  ldbl pr = 1.0L - exp(-pow(lamda * t, a));
  return (sgn(pr) == 0) ? 0.0L : pr;
}

ldbl wbsur(ldbl t, ldbl lamda, int a) // Weibull分布生存関数(片側確率関数)
{
  ldbl pr = exp(-pow(lamda * t, a));
  return (sgn(pr) == 0) ? 0.0L : pr;
}

ldbl wbhaz(ldbl t, ldbl lamda, int a) // Weibull分布ハザード関数
{
  return pow(lamda, a) * (ldbl)a * pow(t, a - 1);
}

ldbl asdis(ldbl h)              // 逆正弦分布関数
{
  return sin(h) * sin(h);
}

ldbl asden(ldbl h)              // 逆正弦分布密度関数
{
  constexpr ldbl SD = sqrt(M_PI * M_PI / 16.0L - 0.5);
  constexpr ldbl Dh = Delta * SD;
  return (asdis(h + 0.5L * Dh) - asdis(h - 0.5L * Dh)) / Delta;
}

ldbl snpinv(ldbl pr)       // 標準正規分布両側確率逆変換関数(Newton法)
{
  if (pr <= 0.0L || pr >= 1.0L) return 0.0L;
  ldbl z = guessz(pr);
  ldbl ppre = ztop(z);
  if (ppre < Pmin)              // 初期値が非常に小さい時
  {
    z = 1.0L;
    ppre = ztop(z);
  }
  while (sgndiff(ppre, pr) != 0 && fabs(ppre - pr) >= Peps)
  {
    ldbl slope = -2.0L * snden(z) / ppre;
    if (sgn(slope) == 0) break; // 収束
    z += log(pr / ppre) / slope;
    if (sgn(z) <= 0) z = Delta * z;
    ppre = ztop(z);
    if (sgn(ppre) <= 0)         // 収束失敗
    {
      z = HUGE_VALL;
      break;
    }
  }
  return z;
}

ldbl snsur(ldbl z)              // 標準正規分布生存関数(片側確率関数)
{
  ldbl pr = 1.0L - sndis(z);
  return (sgn(pr) == 0) ? 0.0L : pr;
}

ldbl snhaz(ldbl z)              // 標準正規分布ハザード関数
{
  return snden(z) / snsur(z);
}

ldbl cpval(ldbl c, int df)      // χ2乗分布片側確率関数
{
  return ctop(c, df);
}

ldbl csur(ldbl c, int df)       // χ2乗分布生存関数(片側確率関数)
{
  return ctop(c, df);
}

ldbl chaz(ldbl c, int df)       // χ2乗分布ハザード関数
{
  return cden(c, df) / csur(c, df);
}

ldbl tpval(ldbl t, int df)      // t分布両側確率関数
{
  if (df < 1) return ztop(t); // df=∞の時は標準正規分布
  return ttop(t, df);
}

ldbl fpval(ldbl f, int df1, int df2) // F分布片側確率関数
{
  if (df2 < 1) return ctop(df1 * f, df1); // df2=∞の時はχ2乗分布利用
  ldbl pr = ftop(f, df1, df2);
  return (sgn(pr) == 0) ? 0.0L : pr;
}

ldbl betaden(ldbl t, ldbl df1, ldbl df2) // ベータ分布密度関数
{
  if (t <= 0.0L || 1.0L <= t || df1 <= 0.0L || df2 <= 0.0L) return 0.0L;
  ldbl f = df2 / df1 * t / (1.0L - t);
  ldbl df = df2 / df1 / (1.0L - t) / (1.0 - t);
  ldbl pr = df * fden(f, 2 * (int)df1, 2 * (int)df2);
  return (sgn(pr) == 0) ? 0.0L : pr;
}

ldbl betadis(ldbl t, ldbl df1, ldbl df2) // ベータ分布関数
{
  if (t <= 0.0L || df1 <= 0.0L || df2 <= 0.0L) return 0.0L;
  if (1.0L <= t) return 1.0L;
  ldbl f = df2 * t / df1 / (1.0L - t);
  ldbl pr = fdis(f, 2 * (int)df1, 2 * (int)df2);
  return (sgn(pr) == 0) ? 0.0L : pr;
}

ldbl poiden(int x, ldbl lamda)  // Poisson分布密度関数
{
  ldbl pr = exp(-lamda);
  for (int i = 1; i <= x; ++i) pr = pr * lamda / (ldbl)i;
  return (sgn(pr) == 0) ? 0.0L : pr;
}

ldbl poidis(int x, ldbl lamda)  // Poisson分布関数
{
  ldbl pri = exp(-lamda);
  ldbl pr = pri;
  for (int i = 1; i <= x; ++i)
  {
    pri = pri * lamda / (ldbl)i;
    pr += pri;
  }
  return (sgn(pr) == 0) ? 0.0L : pr;
}

ldbl fw1den(int n, int t)       // Wilcoxonの符号付順位分布度数関数
{
  int tn = n * (n + 1) / 2;
  if (n <= 0 || t < 0 || t > tn) return 0.0L;
  if (t <= 1 || t >= tn - 1) return 1.0L;
  return fw1den(n - 1, t) + fw1den(n - 1, t - n); // 度数を返す
}

ldbl fw1dis(int n, int t)      // Wilcoxonの符号付順位分布累積度数関数
{
  if (n < 0 || t < 0) return 0.0L;
  if (n == 0 || t == 0) return 1.0L;
  if (t == 1) return 2.0L;
  return fw1dis(n - 1, t) + fw1dis(n - 1, t - n); // 左側累積度数を返す
}

ldbl fuden(int m, int n, int u) // Mann-WhitneyのU分布度数関数
{
  int umax = m * n;
  if (m <= 0 || n <= 0 || u < 0 || u > umax) return 0.0L;
  if (u == 0 || u == 1 || u == umax || u == umax - 1) return 1.0L;
  return fuden(m, n - 1, u) + fuden(m - 1, n, u - n); // 度数を返す
}

ldbl fudis(int m, int n, int u) // Mann-WhitneyのU分布累積度数布関数
{
  if (m < 0 || n < 0 || u < 0) return 0.0L;
  if (m == 0 || n == 0 || u == 0) return 1.0L;
  if (u == 1) return 2.0L;
  return fudis(m, n - 1, u) + fudis(m - 1, n, u - n); // 左側累積度数を返す
}

ldbl nccden(ldbl c, int df, ldbl lamda) // 非心χ2乗分布密度関数
{
  constexpr ldbl DChalf = Delta * 0.5L;
  return (ncctolp(c + DChalf, df, lamda) - ncctolp(c - DChalf, df, lamda)) / Delta;
}

ldbl ncfden(ldbl f, int df1, int df2, ldbl lamda) // 非心F分布密度関数
{
  constexpr ldbl DChalf = Delta * 0.5L;
  return (ncftolp(f + DChalf, df1, df2, lamda) - ncftolp(f - DChalf, df1, df2, lamda)) / Delta;
}

ldbl nctden(ldbl t, int df, ldbl lamda) // 非心t分布密度関数
{
  constexpr ldbl DChalf = Delta * 0.5L;
  return (ncttolp(t + DChalf, df, lamda) - ncttolp(t - DChalf, df, lamda)) / Delta;
}

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  bool v_ct;                    // 連続修正フラッグ
  distype v_dtype;              // 分布の種類
  lstr fparalst;                // 実数パラメータリスト
  vecdbl fparavec;              // 実数パラメータベクトル配列
  funtype v_ftype;              // 関数の種類
  lstr iparalst;                // 整数パラメータリスト
  vecint iparavec;              // 整数パラメータベクトル配列
  pairdbl xrange;               // 描画範囲の最小値と最大値
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  bool ct() const               // 連続修正フラッグを返す
  { return v_ct; }
  distype dtype() const         // 分布の種類を返す
  { return v_dtype; }
  double fpara(size_t idx) const // 実数パラメーターを返す
  { return (idx < fparavec.size()) ? fparavec[idx] : 0.0; }
  funtype ftype() const         // 関数の種類を返す
  { return v_ftype; }
  int ipara(size_t idx) const   // 整数パラメーターを返す
  { return (idx < iparavec.size()) ? iparavec[idx] : 0.0; }
  size_t nfpara() const         // 実数パラメーター数を返す
  { return fparavec.size(); }
  size_t nipara() const         // 整数パラメーター数を返す
  { return iparavec.size(); }
  double rmin() const           // 描画範囲の最小値を返す
  { return xrange.first; }
  double rmax() const           // 描画範囲の最大値を返す
  { return xrange.second; }
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check();                 // オプションのチェックとセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_ct(false), v_dtype(disn), fparalst(), fparavec(), v_ftype(funden), iparalst(), iparavec(), 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); // プログラム名削除
}

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] == "ct")           // 連続修正フラッグの指定
  {
    if (opt[1].empty() || opt[1] == "+") v_ct = true;
    else if (opt[1] == "-") v_ct = false;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "dtype")   // 分布の種類の指定
  {
    if (opt[1].empty() || opt[1] == "n") v_dtype = disn;
    else if (opt[1] == "as") v_dtype = disas;
    else if (opt[1] == "b") v_dtype = disb;
    else if (opt[1] == "beta") v_dtype = disbeta;
    else if (opt[1] == "c") v_dtype = disc;
    else if (opt[1] == "d") v_dtype = disd;
    else if (opt[1] == "dun") v_dtype = disdun;
    else if (opt[1] == "e") v_dtype = dise;
    else if (opt[1] == "etc") v_dtype = disetc;
    else if (opt[1] == "exp") v_dtype = disexp;
    else if (opt[1] == "f") v_dtype = disf;
    else if (opt[1] == "fis") v_dtype = disfis;
    else if (opt[1] == "ln") v_dtype = disln;
    else if (opt[1] == "log") v_dtype = dislog;
    else if (opt[1] == "ncc") v_dtype = disncc;
    else if (opt[1] == "ncf") v_dtype = disncf;
    else if (opt[1] == "npv") v_dtype = disnpv;
    else if (opt[1] == "nct") v_dtype = disnct;
    else if (opt[1] == "per") v_dtype = disper;
    else if (opt[1] == "plog") v_dtype = displog;
    else if (opt[1] == "poi") v_dtype = dispoi;
    else if (opt[1] == "poly") v_dtype = dispoly;
    else if (opt[1] == "ppv") v_dtype = disppv;
    else if (opt[1] == "pro") v_dtype = dispro;
    else if (opt[1] == "sn") v_dtype = dissn;
    else if (opt[1] == "t") v_dtype = dist;
    else if (opt[1] == "tuk") v_dtype = distuk;
    else if (opt[1] == "w1") v_dtype = disw1;
    else if (opt[1] == "w2") v_dtype = disw2;
    else if (opt[1] == "wb") v_dtype = diswb;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "f")       // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "fpara")   // 実数パラメーターの指定
  {
    if (strsplit(fparalst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "ftype")   // 関数の種類の指定
  {
    if (opt[1].empty() || opt[1] == "den") v_ftype = funden;
    else if (opt[1] == "dis") v_ftype = fundis;
    else if (opt[1] == "haz") v_ftype = funhaz;
    else if (opt[1] == "pval") v_ftype = funpval;
    else if (opt[1] == "pinv") v_ftype = funpinv;
    else if (opt[1] == "sur") v_ftype = funsur;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else if (opt[0] == "ipara")   // 整数パラメーターの指定
  {
    if (strsplit(iparalst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "range")  // 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()            // オプションのチェックとセット
{
  if (fparalst.size() > 0)      // 実数パラメータが指定されている時
  {
    fparavec.resize(fparalst.size(), 0.0);
    size_t i = 0;
    for (string fpstr : fparalst)
    {
      double fpara = strtodf(fpstr);
      if (fpara == BLANKf) return disperr(pname(), "実数パラメータの指定が不適切！");
      fparavec[i] = fpara;
      ++i;
    }
  }
  if (iparalst.size() > 0)      // 整数パラメータが指定されている時
  {
    iparavec.resize(iparalst.size(), 0.0);
    size_t i = 0;
    for (string ipstr : iparalst)
    {
      int ipara = (int)(strtodf(ipstr));
      if (ipara == INT_MAX) return disperr(pname(), "整数パラメータの指定が不適切！");
      iparavec[i] = ipara;
      ++i;
    }
  }
  switch (v_dtype)
  {
  case disas :
    if (v_ftype != funden && v_ftype != fundis) return disperr(pname(), "逆正弦分布は分布密度関数と分布関数だけサポート");
    break;
  case disb :
    if (v_ftype != funden && v_ftype != fundis && v_ftype != funpval) return disperr(pname(), "二項分布は分布密度関数と分布関数と確率関数だけサポート");
    break;
  case disbeta :
    if (v_ftype != funden && v_ftype != fundis) return disperr(pname(), "ベータ分布は分布密度関数と分布関数だけサポート");
    break;
  case disc :
    if (v_ftype == funpinv) return disperr(pname(), "χ2乗分布の確率逆関数は非サポート");
    break;
  case disdun :
    if (v_ftype != fundis && v_ftype != funpval) return disperr(pname(), "Dunnettのd分布は分布関数と確率関数だけサポート");
    break;
  case dise :
    if (v_ftype == funpval || v_ftype == funpinv) return disperr(pname(), "指数分布の確率計算と確率逆関数は非サポート");
    break;
  case disf :
    if (v_ftype != funden && v_ftype != fundis && v_ftype != funpval) return disperr(pname(), "F分布は分布密度関数と分布関数と確率関数だけサポート");
    break;
  case disfis :
    if (v_ftype != funden && v_ftype != fundis && v_ftype != funpval) return disperr(pname(), "超幾何分布は分布密度関数と分布関数と確率関数だけサポート");
    break;
  case disln :
    if (v_ftype == funpinv) return disperr(pname(), "対数正規分布の確率逆関数は非サポート");
    break;
  case dislog :
    if (v_ftype != funden && v_ftype != fundis) return disperr(pname(), "ロジスティック分布は分布密度関数と分布関数だけサポート");
    break;
  case disn :
    if (v_ftype != funden && v_ftype != fundis && v_ftype != funpval && v_ftype != funpinv) return disperr(pname(), "正規分布は分布密度関数と分布関数と確率関数と逆確率関数だけサポート");
    break;
  case disncc :
    if (v_ftype != funden && v_ftype != fundis) return disperr(pname(), "非心χ2乗分布は分布密度関数と分布関数だけサポート");
    break;
  case disncf :
    if (v_ftype != funden && v_ftype != fundis) return disperr(pname(), "非心F分布は分布密度関数と分布関数だけサポート");
    break;
  case disnct :
    if (v_ftype != funden && v_ftype != fundis) return disperr(pname(), "非心t分布は分布密度関数と分布関数だけサポート");
    break;
  case dispoi :
    if (v_ftype != funden && v_ftype != fundis) return disperr(pname(), "Poisson分布は分布密度関数と分布関数だけサポート");
    break;
  case dispro :
    if (v_ftype != funden && v_ftype != fundis) return disperr(pname(), "プロビット分布は分布密度関数と分布関数だけサポート");
    break;
  case dissn :
    if (v_ftype != funden && v_ftype != fundis && v_ftype != funpval && v_ftype != funpinv) return disperr(pname(), "標準正規分布は分布密度関数と分布関数と確率関数と逆確率関数だけサポート");
    break;
  case dist :
    if (v_ftype != funden && v_ftype != fundis && v_ftype != funpval) return disperr(pname(), "t分布は分布密度関数と分布関数と確率関数だけサポート");
    break;
  case distuk :
    if (v_ftype != fundis && v_ftype != funpval) return disperr(pname(), "Tukeyのq分布は分布関数と確率関数だけサポート");
    break;
  case disw1 :
    if (v_ftype != funden && v_ftype != fundis && v_ftype != funpval) return disperr(pname(), "Wilcoxonの符号付き順位分布は分布密度関数と分布関数と確率関数だけサポート");
    break;
  case disw2 :
    if (v_ftype != funden && v_ftype != fundis && v_ftype != funpval) return disperr(pname(), "Wilcoxonの順位和分布は分布密度関数と分布関数と確率関数だけサポート");
    break;
  case diswb :
    if (v_ftype == funpval || v_ftype == funpinv) return disperr(pname(), "Weibull分布の確率計算と確率逆関数は非サポート");
  }
  return 0;                     // エラーなし
}

class plotclass                 // メインクラス
{
private:
  void plot_asden(ostream& fout, const option& opt) const;
  // 逆正弦分布密度関数出力
  void plot_asdis(ostream& fout, const option& opt) const;
  // 逆正弦分布関数(逆正弦曲線)出力
  void plot_bden(ostream& fout, const option& opt) const;
  // 二項分布密度関数出力
  void plot_bdis(ostream& fout, const option& opt) const; // 二項分布関数出力
  void plot_bpval(ostream& fout, const option& opt) const;
  // 二項分布片側確率出力
  void plot_betaden(ostream& fout, const option& opt) const;
  // ベータ分布密度関数出力
  void plot_betadis(ostream& fout, const option& opt) const;
  // ベータ分布関数出力
  void plot_cden(ostream& fout, const option& opt) const; // χ2乗分布密度関数出力
  void plot_cdis(ostream& fout, const option& opt) const;
  // χ2乗分布関数出力
  void plot_chaz(ostream& fout, const option& opt) const;
  // χ2乗分布ハザード関数出力
  void plot_cpval(ostream& fout, const option& opt) const;
  // χ2乗分布片側確率出力
  void plot_csur(ostream& fout, const option& opt) const;
  // χ2乗分布生存関数出力
  void plot_delp(ostream& fout, const option& opt) const;
  // Mahalanobisの汎距離楕円出力
  void plot_ddis(ostream& fout, const option& opt) const;
  // Dunnettのd分布関数出力
  void plot_dpval(ostream& fout, const option& opt) const;
  // Dunnettのd分布両側確率出力
  void plot_eden(ostream& fout, const option& opt) const;
  // 指数分布密度関数出力
  void plot_edis(ostream& fout, const option& opt) const; // 指数分布関数出力
  void plot_ehaz(ostream& fout, const option& opt) const;
  // 指数分布ハザード関数出力
  void plot_esur(ostream& fout, const option& opt) const;
  // 指数分布生存関数出力
  void plot_exp(ostream& fout, const option& opt) const; // 指数関数出力
  void plot_fden(ostream& fout, const option& opt) const; // F分布密度関数出力
  void plot_fdis(ostream& fout, const option& opt) const; // F分布関数出力
  void plot_fpval(ostream& fout, const option& opt) const; // F分布片側確率出力
  void plot_fisden(ostream& fout, const option& opt) const;
  // 超幾何分布密度関数出力
  void plot_fisdis(ostream& fout, const option& opt) const; // 超幾何分布関数出力
  void plot_fispval(ostream& fout, const option& opt) const;
  // 超幾何分布片側確率出力(Fisherの正確計算)
  void plot_lnden(ostream& fout, const option& opt) const;
  // 対数正規分布密度関数出力
  void plot_lndis(ostream& fout, const option& opt) const;
  // 対数正規分布関数出力
  void plot_lnhaz(ostream& fout, const option& opt) const;
  // 対数正規分布ハザード関数出力
  void plot_lnpval(ostream& fout, const option& opt) const;
  // 対数正規分布両側確率出力
  void plot_lnsur(ostream& fout, const option& opt) const;
  // 対数正規分布生存関数出力
  void plot_logden(ostream& fout, const option& opt) const;
  // ロジスティック分布密度関数出力
  void plot_logdis(ostream& fout, const option& opt) const;
  // ロジスティック分布関数出力(ロジスティック曲線)
  void plot_nden(ostream& fout, const option& opt) const; // 正規分布密度関数出力
  void plot_ndis(ostream& fout, const option& opt) const; // 正規分布関数出力
  void plot_npval(ostream& fout, const option& opt) const;
  // 正規分布両側確率出力
  void plot_npinv(ostream& fout, const option& opt) const;
  // 正規分布両側確率逆関数出力
  void plot_nccden(ostream& fout, const option& opt) const;
  // 非心χ2乗分布密度関数出力
  void plot_nccdis(ostream& fout, const option& opt) const;
  // 非心χ2乗分布関数出力
  void plot_ncfden(ostream& fout, const option& opt) const;
  // 非心F分布密度関数出力
  void plot_ncfdis(ostream& fout, const option& opt) const; // 非心F分布関数出力
  void plot_nctden(ostream& fout, const option& opt) const;
  // 非心t分布密度関数出力
  void plot_nctdis(ostream& fout, const option& opt) const; // 非心t分布関数出力
  void plot_perdis(ostream& fout, const option& opt) const; // 周期曲線出力
  void plot_plogdis(ostream& fout, const option& opt) const;
  // ロジスティック関数の陽性率曲線出力
  void plot_poiden(ostream& fout, const option& opt) const;
  // Poisson分布密度関数出力
  void plot_poidis(ostream& fout, const option& opt) const;
  // Poisson分布関数出力
  void plot_poly(ostream& fout, const option& opt) const; // 多項関数出力
  double cal_poly(const option& opt, int n, double x) const;
  // 多項関数値を計算して返す
  void plot_npv(ostream& fout, const option& opt) const; // NPV曲線出力
  void plot_ppv(ostream& fout, const option& opt) const; // PPV曲線出力
  void plot_proden(ostream& fout, const option& opt) const;
  // プロビット分布密度関数出力(正規分布)
  void plot_prodis(ostream& fout, const option& opt) const;
  // プロビット分布関数出力(プロビット曲線)
  void plot_snden(ostream& fout, const option& opt) const;
  // 標準正規分布密度関数出力
  void plot_sndis(ostream& fout, const option& opt) const;
  // 標準正規分布関数出力
  void plot_snpval(ostream& fout, const option& opt) const;
  // 標準正規分布両側確率出力
  void plot_snpinv(ostream& fout, const option& opt) const;
  // 標準正規分布両側確率逆関数出力
  void plot_tden(ostream& fout, const option& opt) const; // t分布密度関数出力
  void plot_tdis(ostream& fout, const option& opt) const; // t分布関数出力
  void plot_tpval(ostream& fout, const option& opt) const; // t分布両側確率出力
  void plot_qdis(ostream& fout, const option& opt) const;
  // Tukeyのq分布関数出力
  void plot_qpval(ostream& fout, const option& opt) const;
  // Tukeyのq分布両側確率出力
  void plot_w1den(ostream& fout, const option& opt) const;
  // Wilcoxonの符号付順位分布密度関数出力
  void plot_w1dis(ostream& fout, const option& opt) const;
  // Wilcoxonの符号付順位分布関数出力
  void plot_w1pval(ostream& fout, const option& opt) const;
  // Wilcoxonの符号付順位分布両側確率出力
  void plot_w2den(ostream& fout, const option& opt) const;
  // Wilcoxonの順位和分布密度関数出力
  void plot_w2dis(ostream& fout, const option& opt) const;
  // Wilcoxonの順位和分布関数出力
  void plot_w2pval(ostream& fout, const option& opt) const;
  // Wilcoxonの順位和分布両側確率出力
  void plot_wbden(ostream& fout, const option& opt) const;
  // Weibull分布密度関数出力
  void plot_wbdis(ostream& fout, const option& opt) const;
  // Weibull分布関数出力
  void plot_wbhaz(ostream& fout, const option& opt) const;
  // Weibull分布ハザード関数出力
  void plot_wbsur(ostream& fout, const option& opt) const;
  // Weibull分布生存関数出力
  void plot_etc(ostream& fout, const option& opt) const; // その他の関数出力
public:
  plotclass()                   // コンストラクタ
  {}
  ~plotclass()                  // デストラクタ
  {}
  void put(option& opt);        // 各種分布を出力する
};

void plotclass::plot_asden(ostream& fout, const option& opt) const
// 逆正弦分布密度関数出力
{
  constexpr ldbl SD = sqrt(M_PI * M_PI / 16.0L - 0.5);
  constexpr ldbl HMAXdef = M_PI * 0.5L; // θ最大値
  constexpr ldbl Hmean = M_PI * 0.25L;  // θ平均値
  ldbl hmin = 0.0L, hmax = HMAXdef;
  if (opt.rmin() != BLANKf && opt.rmin() >= 0.0) hmin = opt.rmin();
  if (opt.rmax() != BLANKf && opt.rmax() > 0.0) hmax = opt.rmax();
  ldbl dh = (hmax - hmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%逆正弦分布密度関数:1SD=" << SD << "(RAD)\n";
  ldbl h = hmin;
  for (uint i = 0; i <= RDdef; ++i, h += dh)
  {
    ldbl x = 5.0L + (h - Hmean) / SD;
    fout << x << " " << asden(h) << "\n";
  }
}

void plotclass::plot_asdis(ostream& fout, const option& opt) const
// 逆正弦分布関数(逆正弦曲線)出力
{
  constexpr ldbl SD = sqrt(M_PI * M_PI / 16.0L - 0.5L);
  constexpr ldbl HMAXdef = M_PI * 0.5L; // θ最大値
  constexpr ldbl Hmean = M_PI * 0.25L; // θ平均値
  ldbl hmin = 0.0L, hmax = HMAXdef;
  if (opt.rmin() != BLANKf && opt.rmin() >= 0.0) hmin = opt.rmin();
  if (opt.rmax() != BLANKf && opt.rmax() > 0.0) hmax = opt.rmax();
  ldbl dh = (hmax - hmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%逆正弦曲線:1SD=" << SD << "(RAD)\n";
  ldbl h = hmin;
  for (uint i = 0; i <= RDdef; ++i, h += dh)
  {
    ldbl x = 5.0L + (h - Hmean) / SD;
    fout << x << " " << asdis(h) << "\n";
  }
}

void plotclass::plot_bden(ostream& fout, const option& opt) const
// 二項分布密度関数出力
{
  constexpr int Ndef = 10;      // 省略時例数
  constexpr ldbl P0def = 0.5L;  // 省略時理論確率
  int n = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : Ndef; // 全例数
  ldbl p0 = (opt.nfpara() > 0 && 0.0 < opt.fpara(0) && opt.fpara(0) < 1.0) ? (ldbl)(opt.fpara(0)) : P0def; // 理論確率
  int imin = 0, imax = n;
  if (opt.rmin() != BLANKf && opt.rmin() >= 0.0) imin = (int)(opt.rmin());
  if (opt.rmax() != BLANKf && opt.rmax() <= (double)n) imax = (int)(opt.rmax());
  fout << "-s1 -drectangle_solid_fill -x1 -y3 -w80\n";
  fout << "%二項分布密度関数:理論確率=" << p0 << " n=" << n << "\n";
  for (int i = imin; i <= imax; ++i)
  {
    double x = (double)i;
    fout << x - 0.1 << " 0 " << x + 0.1 << " " << bnmden(i, n, p0) << "\n";
  }
}

void plotclass::plot_bdis(ostream& fout, const option& opt) const
// 二項分布関数出力
{
  constexpr int Ndef = 10;      // 省略時例数
  constexpr ldbl P0def = 0.5L;  // 省略時理論確率
  int n = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : Ndef; // 全例数
  ldbl p0 = (opt.nfpara() > 0 && 0.0 < opt.fpara(0) && opt.fpara(0) < 1.0) ? (ldbl)(opt.fpara(0)) : P0def; // 理論確率
  int imin = 0, imax = n;
  if (opt.rmin() != BLANKf && opt.rmin() >= 0.0) imin = (int)(opt.rmin());
  if (opt.rmax() != BLANKf && opt.rmax() <= (double)n) imax = (int)(opt.rmax());
  fout << "-s1 -dline -w80\n";
  fout << "%二項分布関数:理論確率=" << p0 << " n=" << n << "\n";
  ldbl ppre = 0.0L;
  for (int i = imin; i <= imax; ++i)
  {
    ldbl pr = bnmdis(i, n, p0);
    fout << i << " " << ppre << "\n";
    fout << i << " " << pr << "\n";
    ppre = pr;
  }
}

void plotclass::plot_bpval(ostream& fout, const option& opt) const
// 二項分布片側確率出力
{
  constexpr int Ndef = 10;      // 省略時例数
  constexpr ldbl P0def = 0.5L;  // 省略時理論確率
  int n = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : Ndef; // 全例数
  ldbl p0 = (opt.nfpara() > 0 && 0.0 < opt.fpara(0) && opt.fpara(0) < 1.0) ? (ldbl)(opt.fpara(0)) : P0def; // 理論確率
  int imin = 0, imax = n;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin() && opt.rmin() < (double)n) imin = int(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax() && opt.rmax() <= (double)n) imax = (int)(opt.rmax());
  fout << "-s1 -dline -w80\n";
  fout << "%二項分布片側確率:理論確率=" << p0 << " n=" << n << "\n";
  if (imin == 0)
  {
    fout << 0 << " " << 1.0 << "\n";
    imin = 1;
  }
  for (int i = imin; i <= imax; ++i)
  { fout << i << " " << 1.0 - bnmdis(i - 1, n, p0) << "\n"; }
}

void plotclass::plot_betaden(ostream& fout, const option& opt) const
// ベータ分布密度関数出力
{
  constexpr ldbl DF1def = 5.0L; // 省略時第1自由度
  constexpr ldbl DF2def = 5.0L; // 省略時第2自由度
  ldbl df1 = (opt.nfpara() > 0 && 0.0 < opt.fpara(0)) ? (ldbl)(opt.fpara(0)) : DF1def; // 第1自由度
  ldbl df2 = (opt.nfpara() > 1 && 0.0 <= opt.fpara(1)) ? (ldbl)(opt.fpara(1)) : DF2def; // 第2自由度
  ldbl tmin = 0.0, tmax = 1.0;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin() && opt.rmin() <= 1.0) tmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && tmin < opt.rmax() && opt.rmax() <= 1.0) tmax = (ldbl)(opt.rmax());
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%ベータ分布密度関数:df1=" << df1 << " df2=" << df2 << "\n";
  ldbl dt = (tmax - tmin) / (ldbl)RDdef, t = tmin;
  for (uint i = 0; i <= RDdef; ++i, t += dt)
  { fout << t << " " << betaden(t, df1, df2) << "\n"; }
}

void plotclass::plot_betadis(ostream& fout, const option& opt) const
// ベータ分布関数出力
{
  constexpr ldbl DF1def = 5.0L; // 省略時第1自由度
  constexpr ldbl DF2def = 5.0L; // 省略時第2自由度
  ldbl df1 = (opt.nfpara() > 0 && 0.0 < opt.fpara(0)) ? (ldbl)(opt.fpara(0)) : DF1def; // 第1自由度
  ldbl df2 = (opt.nfpara() > 1 && 0.0 <= opt.fpara(1)) ? (ldbl)(opt.fpara(1)) : DF2def; // 第2自由度
  ldbl tmin = 0.0, tmax = 1.0;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin() && opt.rmin() <= 1.0) tmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && tmin < opt.rmax() && opt.rmax() <= 1.0) tmax = (ldbl)(opt.rmax());
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%ベータ分布関数:df1=" << df1 << " df2=" << df2 << "\n";
  ldbl dt = (tmax - tmin) / (ldbl)RDdef, t = tmin;
  for (uint i = 0; i <= RDdef; ++i, t += dt)
  { fout << t << " " << betadis(t, df1, df2) << "\n"; }
}

void plotclass::plot_cden(ostream& fout, const option& opt) const
// χ2乗分布密度関数出力
{
  constexpr int DFdef = 1;                       // 省略時自由度
  constexpr ldbl CMAXdef = 20.0L;                // 省略時χ2乗最大値
  constexpr ldbl LCMINdef = log(CMAXdef / 1000.0L); // 対数χ2乗最小値
  int df = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : DFdef;
  ldbl cmin = 0.0, cmax = CMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) cmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) cmax = (ldbl)(opt.rmax());
  ldbl lcmin = (cmin > 0.0L) ? log(cmin) : LCMINdef, lcmax = log(cmax);
  ldbl ldc = (lcmax - lcmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%χ2乗分布密度関数:df=" << df << "\n";
  if (cmin == 0.0) fout << "0 " << cden(cmin, df) << "\n";
  ldbl lc = lcmin;
  for (uint i = 0; i <= RDdef; ++i, lc += ldc)
  {
    ldbl c = exp(lc);
    fout << c << " " << cden(c, df) << "\n";
  }
}

void plotclass::plot_cdis(ostream& fout, const option& opt) const
// χ2乗分布関数出力
{
  constexpr int DFdef = 1;        // 省略時自由度
  constexpr ldbl CMAXdef = 20.0L; // 省略時χ2乗最大値
  constexpr ldbl LCMINdef = log(CMAXdef / 1000.0L); // 対数χ2乗最小値
  int df = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : DFdef;
  ldbl cmin = 0.0, cmax = CMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) cmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) cmax = (ldbl)(opt.rmax());
  ldbl lcmin = (cmin > 0.0L) ? log(cmin) : LCMINdef, lcmax = log(cmax);
  ldbl ldc = (lcmax - lcmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%χ2乗分布関数:df=" << df << "\n";
  if (cmin == 0.0) fout << "0 0\n";
  ldbl lc = lcmin;
  for (uint i = 0; i <= RDdef; ++i, lc += ldc)
  {
    ldbl c = exp(lc);
    fout << c << " " << cdis(c, df) << "\n";
  }
}

void plotclass::plot_chaz(ostream& fout, const option& opt) const
// χ2乗分布ハザード関数出力
{
  constexpr int DFdef = 1;        // 省略時自由度
  constexpr ldbl CMAXdef = 20.0L; // 省略時χ2乗最大値
  int df = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : DFdef;
  ldbl cmin = 0.0, cmax = CMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) cmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) cmax = (ldbl)(opt.rmax());
  ldbl dc = (cmax - cmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%χ2乗分布ハザード関数:df=" << df << "\n";
  ldbl c = cmin;
  for (uint i = 0; i <= RDdef; ++i, c += dc)
  { fout << c << " " << chaz(c, df) << "\n"; }
}

void plotclass::plot_cpval(ostream& fout, const option& opt) const
// χ2乗分布片側確率出力
{
  constexpr int DFdef = 1;        // 省略時自由度
  constexpr ldbl CMAXdef = 20.0L; // 省略時χ2乗最大値
  constexpr ldbl LCMINdef = log(CMAXdef / 1000.0L); // 対数χ2乗最小値
  int df = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : DFdef;
  ldbl cmin = 0.0, cmax = CMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) cmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) cmax = (ldbl)(opt.rmax());
  ldbl lcmin = (cmin > 0.0) ? log(cmin) : LCMINdef, lcmax = log(cmax);
  ldbl ldc = (lcmax - lcmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%χ2乗分布片側確率:df=" << df << "\n";
  if (cmin == 0.0) fout << "0 1\n";
  ldbl lc = lcmin;
  for (uint i = 0; i <=RDdef; ++i, lc += ldc)
  {
    ldbl c = exp(lc);
    fout << c << " " << cpval(c, df) << "\n";
  }
}

void plotclass::plot_csur(ostream& fout, const option& opt) const
// χ2乗分布生存関数出力
{
  constexpr int DFdef = 1;        // 省略時自由度
  constexpr ldbl CMAXdef = 20.0L; // 省略時χ2乗最大値
  constexpr ldbl LCMINdef = log(CMAXdef / 1000.0L); // 対数χ2乗最小値
  int df = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : DFdef;
  ldbl cmin = 0.0, cmax = CMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) cmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) cmax = (ldbl)(opt.rmax());
  ldbl lcmin = (cmin > 0.0L) ? log(cmin) : LCMINdef, lcmax = log(cmax);
  ldbl ldc = (lcmax - lcmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%χ2乗分布生存関数:df=" << df << "\n";
  if (cmin == 0.0) fout << "0 1\n";
  ldbl lc = lcmin;
  for (uint i = 0; i <= RDdef; ++i, lc += ldc)
  {
    ldbl c = exp(lc);
    fout << c << " " << csur(c, df) << "\n";
  }
}

void plotclass::plot_delp(ostream& fout, const option& opt) const
// Mahalanobisの汎距離楕円出力
{
  constexpr ldbl Ddef = 1.0L;            // 省略時D
  constexpr ldbl Rdef = 0.0L;            // 省略時相関係数
  constexpr ldbl TMAXdef = 2.0L * M_PI;  // 省略時t最大値
  ldbl d = (opt.nfpara() > 0 && 0.0 < opt.fpara(0)) ? opt.fpara(0) : Ddef;
  ldbl r = (opt.nfpara() > 1 && 0.0 <= opt.fpara(1) && opt.fpara(1) <= 1.0) ? opt.fpara(1) : Rdef;
  int df = 10000;               // 本来は無限大
  double a = 0.0, b = 0.0, h = 0.0;
  double cc = 1.0 - (double)(ftop((ldbl)(0.5 * d * d), 2, df));
  calelp(a, b, h, cc, 2, df, 1.0, 1.0, r);
  double tmin = 0, tmax = TMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin() && opt.rmin() < TMAXdef) tmin = opt.rmin();
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax() && opt.rmax() <= TMAXdef) tmax = opt.rmax();
  double dt = (tmax - tmin) / (double)RDdef;
  fout << "-s1 -dcurve,spline_close -w80\n";
  fout << "%Mahalanobisの汎距離楕円:D=" << d << "(信頼係数=" << cc << ") r=" << r << "\n";
  double t = tmin;
  for (uint i = 0; i <= RDdef; ++i, t += dt)
  {
    double x0 = a * cos(t), y0 = b * sin(t);
    double x = x0 * cos(h) - y0 * sin(h), y = x0 * sin(h) + y0 * cos(h);
    fout << x << " " << y << "\n";
  }
}

void plotclass::plot_ddis(ostream& fout, const option& opt) const
// Dunnettのd分布関数出力
{
  constexpr int Gdef = 2;       // 省略時群数
  constexpr int DFdef = 0;      // 省略時自由度
  constexpr ldbl Rdef = 0.5L;   // 省略時相関係数
  constexpr ldbl DMINdef = -5.0L; // 省略時dmin
  constexpr ldbl DMAXdef = 5.0L;  // 省略時dmax
  int g = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : Gdef;
  int df = (opt.nipara() > 1 && opt.ipara(1) >= 0) ? opt.ipara(1) : DFdef;
  ldbl r = (opt.nfpara() > 0 && 0.0 <= opt.fpara(0) && opt.fpara(0) <= 1.0) ? (ldbl)(opt.fpara(0)) : Rdef;
  ldbl dmin = DMINdef, dmax = DMAXdef;
  if (opt.rmin() != BLANKf) dmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) dmax = (ldbl)(opt.rmax());
  ldbl dd = (dmax - dmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%Dunnettのd分布関数:g=" << g << " df=" << df << " r=" << r << "\n";
  ldbl d = dmin;
  const ldbl Deps = dd * 0.01L; // d計算誤差
  for (uint i = 0; i <= RDdef; ++i, d += dd)
  {
    if (fabs(d) < Deps) d = 0.0L;
    fout << d << " " << ddis(d, g, df, r) << "\n";
  }
}

void plotclass::plot_dpval(ostream& fout, const option& opt) const
// Dunnettのd分布両側確率出力
{
  constexpr int Gdef = 2;       // 省略時群数
  constexpr int DFdef = 0;      // 省略時自由度
  constexpr ldbl Rdef = 0.5L;   // 省略時相関係数
  constexpr ldbl DMAXdef = 5.0L; // 省略時dmax
  int g = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : Gdef;
  int df = (opt.nipara() > 1 && opt.ipara(1) >= 0) ? opt.ipara(1) : DFdef;
  ldbl r = (opt.nfpara() > 0 && 0.0 <= opt.fpara(0) && opt.fpara(0) <= 1.0) ? (ldbl)(opt.fpara(0)) : Rdef;
  ldbl dmin = 0.0, dmax = DMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) dmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) dmax = (ldbl)(opt.rmax());
  ldbl dd = (dmax - dmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%Dunnettのd分布両側確率:g=" << g << " df=" << df << " r=" << r << "\n";
  ldbl d = dmin;
  const ldbl Deps = dd * 0.01L; // d計算誤差
  for (uint i = 0; i <= RDdef; ++i, d += dd)
  {
    if (fabs(d) < Deps) d = 0.0L;
    fout << d << " " << dpval(d, g, df, r) << "\n";
  }
}

void plotclass::plot_eden(ostream& fout, const option& opt) const
// 指数分布密度関数出力
{
  constexpr ldbl Ldef = 1.0L;      // 省略時λ
  constexpr ldbl TMAXdef = 20.0L;  // 省略時t最大値
  constexpr ldbl EPSt = 0.001L;    // t微小値
  ldbl lamda = (opt.nfpara() > 0 && opt.fpara(0) > 0.0) ? opt.fpara(0) : Ldef;
  ldbl tmin = 0.0, tmax = TMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) tmin = opt.rmin();
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) tmax = opt.rmax();
  ldbl ltmin = log((tmin > 0.0L) ? tmin : tmax * EPSt), ltmax = log(tmax);
  ldbl ldt = (ltmax - ltmin) / ldbl(RDdef);
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%指数分布密度関数:λ=" << lamda << "\n";
  if (tmin == 0.0) fout << "0 " << eden(tmin, lamda) << "\n";
  ldbl lt = ltmin;
  for (uint i = 0; i <= RDdef; ++i, lt += ldt)
  {
    ldbl t = exp(lt);
    fout << t << " " << eden(t, lamda) << "\n";
  }
}

void plotclass::plot_edis(ostream& fout, const option& opt) const
// 指数分布関数出力
{
  constexpr ldbl Ldef = 1.0;     // 省略時λ
  constexpr ldbl TMAXdef = 20.0; // 省略時t最大値
  constexpr ldbl EPSt = 0.001;   // tの微小値
  ldbl lamda = (opt.nfpara() > 0 && opt.fpara(0) > 0.0) ? opt.fpara(0) : Ldef;
  ldbl tmin = 0.0, tmax = TMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) tmin = opt.rmin();
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) tmax = opt.rmax();
  ldbl ltmin = log((tmin > 0.0) ? tmin : tmax * EPSt), ltmax = log(tmax);
  ldbl ldt = (ltmax - ltmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%指数分布関数:λ=" << lamda << "\n";
  if (tmin == 0.0) fout << "0 " << edis(tmin, lamda) << "\n";
  ldbl lt = ltmin;
  for (uint i = 0; i <= RDdef; ++i, lt += ldt)
  {
    ldbl t = exp(lt);
    fout << t << " " << edis(t, lamda) << "\n";
  }
}

void plotclass::plot_ehaz(ostream& fout, const option& opt) const
// 指数分布ハザード関数出力
{
  constexpr ldbl Ldef = 1.0L;     // 省略時λ
  constexpr ldbl TMAXdef = 20.0L; // 省略時tmax
  ldbl lamda = (opt.nfpara() > 0 && opt.fpara(0) > 0.0) ? opt.fpara(0) : Ldef;
  ldbl tmin = 0.0, tmax = TMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) tmin = opt.rmin();
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) tmax = opt.rmax();
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%指数分布ハザード関数:λ=" << lamda << "\n";
  fout << tmin << " " << ehaz(tmin, lamda) << "\n";
  fout << tmax << " " << ehaz(tmax, lamda) << endl;
}

void plotclass::plot_esur(ostream& fout, const option& opt) const
// 指数分布生存関数出力
{
  constexpr ldbl Ldef = 1.0L;     // 省略時λ
  constexpr ldbl TMAXdef = 20.0L; // 省略時t最大値
  constexpr ldbl EPSt = 0.001L;   // tの微小値
  ldbl lamda = (opt.nfpara() > 0 && opt.fpara(0) > 0.0) ? opt.fpara(0) : Ldef;
  ldbl tmin = 0.0, tmax = TMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) tmin = opt.rmin();
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) tmax = opt.rmax();
  ldbl ltmin = log((tmin > 0.0) ? tmin : tmax * EPSt), ltmax = log(tmax);
  ldbl ldt = (ltmax - ltmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%指数分布生存関数:λ=" << lamda << "\n";
  if (tmin == 0.0) fout << "0 " << esur(tmin, lamda) << "\n";
  ldbl lt = ltmin;
  for (uint i = 0; i <= RDdef; ++i, lt += ldt)
  {
    ldbl t = exp(lt);
    fout << t << " " << esur(t, lamda) << "\n";
  }
}

void plotclass::plot_exp(ostream& fout, const option& opt) const
// 指数関数出力
{
  constexpr ldbl Y0def = 1.0L;    // 省略時初期値
  constexpr ldbl Kdef = 1.0L;     // 省略時k
  constexpr ldbl TMAXdef = 20.0L; // 省略時t最大値
  constexpr ldbl EPSt = 0.001L;   // t微小値
  ldbl y0 = (opt.nfpara() > 0 && opt.fpara(0) > 0.0) ? opt.fpara(0) : Y0def;
  ldbl k = (opt.nfpara() > 1 && opt.fpara(1) != 0.0) ? opt.fpara(1) : Kdef;
  ldbl tmin = 0.0, tmax = TMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) tmin = opt.rmin();
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) tmax = opt.rmax();
  ldbl ltmin = log((tmin > 0.0) ? tmin : tmax * EPSt), ltmax = log(tmax);
  ldbl ldt = (ltmax - ltmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%指数関数:y0=" << y0 << " k=" << k << "\n";
  if (tmin == 0.0) fout << "0 " << y0 << "\n";
  ldbl lt = ltmin;
  for (uint i = 0; i <= RDdef; ++i, lt += ldt)
  {
    ldbl t = exp(lt);
    fout << t << " " << y0 * exp(k * t) << "\n";
  }
}

void plotclass::plot_fden(ostream& fout, const option& opt) const
// F分布密度関数出力
{
  constexpr int DF1def = 1;     // 省略時DF1
  constexpr int DF2def = 0;     // 省略時DF2
  constexpr ldbl FMAXdef = 10.0L; // 省略時F最大値
  constexpr ldbl EPSf = 0.001L;   // Fの微小値
  int df1 = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : DF1def;
  int df2 = (opt.nipara() > 1 && opt.ipara(1) >= 0) ? opt.ipara(1) : DF2def;
  ldbl fmin = 0.0, fmax = FMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) fmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) fmax = (ldbl)(opt.rmax());
  ldbl lfmin = log((fmin > 0.0) ? fmin : fmax * EPSf), lfmax = log(fmax);
  ldbl ldf = (lfmax - lfmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%F分布密度関数:df1=" << df1 << " df2=" << df2 << "\n";
  if (fmin == 0.0) fout << "0 " << fden(fmin, df1, df2) << "\n";
  ldbl lf = lfmin;
  for (uint i = 0; i <= RDdef; ++i, lf += ldf)
  {
    ldbl f = exp(lf);
    fout << f << " " << fden(f, df1, df2) << "\n";
  }
}

void plotclass::plot_fdis(ostream& fout, const option& opt) const
// f分布関数出力
{
  constexpr int DF1def = 1;     // 省略時DF1
  constexpr int DF2def = 0;     // 省略時DF2
  constexpr ldbl FMAXdef = 10.0L; // 省略時F最大値
  constexpr ldbl EPSf = 0.001L;   // Fの微小値
  int df1 = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : DF1def;
  int df2 = (opt.nipara() > 1 && opt.ipara(1) >= 0) ? opt.ipara(1) : DF2def;
  ldbl fmin = 0.0, fmax = FMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) fmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) fmax = (ldbl)(opt.rmax());
  ldbl lfmin = log((fmin > 0.0) ? fmin : fmax * EPSf), lfmax = log(fmax);
  ldbl ldf = (lfmax - lfmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%F分布関数:df1=" << df1 << " df2=" << df2 << "\n";
  if (fmin == 0.0) fout << "0 " << fdis(fmin, df1, df1) << "\n";
  ldbl lf = lfmin;
  for (uint i = 0; i <= RDdef; ++i, lf += ldf)
  {
    ldbl f = exp(lf);
    fout << f << " " << fdis(f, df1, df2) << "\n";
  }
}

void plotclass::plot_fpval(ostream& fout, const option& opt) const
// F分布片側確率出力
{
  constexpr int DF1def = 1;     // 省略時DF1
  constexpr int DF2def = 0;     // 省略時DF2
  constexpr ldbl FMAXdef = 10.0L; // 省略時F最大値
  constexpr ldbl EPSf = 0.001L;   // Fの微小値
  int df1 = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : DF1def;
  int df2 = (opt.nipara() > 1 && opt.ipara(1) >= 0) ? opt.ipara(1) : DF2def;
  ldbl fmin = 0.0, fmax = FMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) fmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) fmax = (ldbl)(opt.rmax());
  ldbl lfmin = log((fmin > 0.0) ? fmin : fmax * EPSf), lfmax = log(fmax);
  ldbl ldf = (lfmax - lfmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%F分布片側確率:df1=" << df1 << " df2=" << df2 << "\n";
  if (fmin == 0.0) fout << "0 " << fpval(fmin, df1, df1) << "\n";
  ldbl lf = lfmin;
  for (uint i = 0; i <= RDdef; ++i, lf += ldf)
  {
    ldbl f = exp(lf);
    fout << f << " " << fpval(f, df1, df2) << "\n";
  }
}

void plotclass::plot_fisden(ostream& fout, const option& opt) const
// 超幾何分布密度関数出力
{
  constexpr int Fdef = 1;       // 省略時度数
  int a = (opt.nipara() > 0 && opt.ipara(0) >= 0) ? opt.ipara(0) : Fdef;
  int b = (opt.nipara() > 1 && opt.ipara(1) >= 0) ? opt.ipara(1) : Fdef;
  int c = (opt.nipara() > 2 && opt.ipara(2) >= 0) ? opt.ipara(2) : Fdef;
  int d = (opt.nipara() > 3 && opt.ipara(3) >= 0) ? opt.ipara(3) : Fdef;
  if (a > d) swap(a, d);
  if (b > c) swap(b, c);
  if (b * c > a * d) { swap(a, b);  swap(c, d); }
  int s = a + b, f = c + d, m = a + c, n = b + d;
  int imin = 0, imax = s;
  if (opt.rmin() != BLANKf && opt.rmin() >= 0.0) imin = (int)(opt.rmin());
  if (opt.rmax() != BLANKf && opt.rmax() <= s) imax = (int)(opt.rmax());
  fout << "-s1 -drectangle_solid_fill -x1 -y3 -w80\n";
  fout << "%超幾何分布密度関数:a=" << a << " b=" << b << " c=" << c << " d=" << d << "\n";
  for (int i = imin; i <= imax; ++i)
  {
    double x = (double)i;
    fout << x - 0.1 << " 0 " << x + 0.1 << " " << fisden(i, s - i, m - i, n - s + i) << "\n";
  }
}

void plotclass::plot_fisdis(ostream& fout, const option& opt) const
// 超幾何分布関数出力
{
  constexpr int Fdef = 1;       // 省略時度数
  int a = (opt.nipara() > 0 && opt.ipara(0) >= 0) ? opt.ipara(0) : Fdef;
  int b = (opt.nipara() > 1 && opt.ipara(1) >= 0) ? opt.ipara(1) : Fdef;
  int c = (opt.nipara() > 2 && opt.ipara(2) >= 0) ? opt.ipara(2) : Fdef;
  int d = (opt.nipara() > 3 && opt.ipara(3) >= 0) ? opt.ipara(3) : Fdef;
  if (a > d) swap(a, d);
  if (b > c) swap(b, c);
  if (b * c > a * d) { swap(a, b);  swap(c, d); }
  int s = a + b, f = c + d, m = a + c, n = b + d;
  int imin = 0, imax = s;
  if (opt.rmin() != BLANKf && opt.rmin() >= 0.0) imin = (int)(opt.rmin());
  if (opt.rmax() != BLANKf && opt.rmax() <= s) imax = (int)(opt.rmax());
  fout << "-s1 -dline -w80\n";
  fout << "%超幾何分布関数:a=" << a << " b=" << b << " c=" << c << " d=" << d << "\n";
  double ppre = 0.0;
  for (int i = imin; i <= imax; ++i)
  {
    double p = fisdis(i, m - i, s - i, n - s + i);
    fout << i << " " << ppre << "\n";
    fout << i << " " << p << "\n";
    ppre = p;
  }
}

void plotclass::plot_fispval(ostream& fout, const option& opt) const
// 超幾何分布片側確率出力(Fisherの正確検定)
{
  constexpr int Fdef = 1;       // 省略時度数
  int a = (opt.nipara() > 0 && opt.ipara(0) >= 0) ? opt.ipara(0) : Fdef;
  int b = (opt.nipara() > 1 && opt.ipara(1) >= 0) ? opt.ipara(1) : Fdef;
  int c = (opt.nipara() > 2 && opt.ipara(2) >= 0) ? opt.ipara(2) : Fdef;
  int d = (opt.nipara() > 3 && opt.ipara(3) >= 0) ? opt.ipara(3) : Fdef;
  if (a > d) swap(a, d);
  if (b > c) swap(b, c);
  if (b * c > a * d) { swap(a, b);  swap(c, d); }
  int s = a + b, f = c + d, m = a + c, n = b + d;
  int imin = 0, imax = s;
  if (opt.rmin() != BLANKf && opt.rmin() >= 0.0) imin = (int)(opt.rmin());
  if (opt.rmax() != BLANKf && opt.rmax() <= s) imax = (int)(opt.rmax());
  fout << "-s1 -dline -w80\n";
  fout << "%超幾何分布片側確率:a=" << a << " b=" << b << " c=" << c << " d=" << d << "\n";
  if (imin == 0) fout << "0 1\n";
  for (int i = imin; i < imax; ++i)
  { fout << i + 1 << " " << 1.0 - fisdis(i, m - i, s - i, n - s + i) << "\n"; }
}

void plotclass::plot_nccden(ostream& fout, const option& opt) const
// 非心χ2乗分布密度関数出力
{
  constexpr int DFdef = 1;      // 省略時自由度
  constexpr ldbl Ldef = 0.0L;   // 省略時λ
  constexpr ldbl CMAXdef = 30.0L; // 省略時χ2乗最大値
  constexpr ldbl EPSc = 0.001L;   // χ2乗微小値
  int df = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : DFdef;
  ldbl lamda = (opt.nfpara() > 0 && opt.fpara(0) >= 0.0) ? (ldbl)(opt.fpara(0)) : Ldef;
  ldbl cmin = 0.0, cmax = CMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) cmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) cmax = (ldbl)(opt.rmax());
  ldbl lcmin = log((cmin > 0.0) ? cmin : cmax * EPSc), lcmax = log(cmax);
  ldbl ldc = (lcmax - lcmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%非心χ2乗分布密度関数:df=" << df << " 非心度=" << lamda << "\n";
  if (cmin == 0.0L) fout << "0 0\n";
  ldbl lc = lcmin;
  for (uint i = 0; i <= RDdef; ++i, lc += ldc)
  {
    ldbl c = exp(lc);
    fout << c << " " << nccden(c, df, lamda) << "\n";

  }
}

void plotclass::plot_nccdis(ostream& fout, const option& opt) const
// 非心χ2乗分布関数出力
{
  constexpr int DFdef = 1;      // 省略時自由度
  constexpr ldbl Ldef = 0.0L;   // 省略時λ
  constexpr ldbl CMAXdef = 30.0L; // 省略時χ2乗最大値
  constexpr ldbl EPSc = 0.001L;   // χ2乗微小値
  int df = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : DFdef;
  ldbl lamda = (opt.nfpara() > 0 && opt.fpara(0) >= 0.0) ? (ldbl)(opt.fpara(0)) : Ldef;
  ldbl cmin = 0.0, cmax = CMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) cmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) cmax = (ldbl)(opt.rmax());
  ldbl lcmin = log((cmin > 0.0) ? cmin : cmax * EPSc), lcmax = log(cmax);
  ldbl ldc = (lcmax - lcmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%非心χ2乗分布関数:df=" << df << " 非心度=" << lamda << "\n";
  if (cmin == 0.0L) fout << "0 0\n";
  ldbl lc = lcmin;
  for (uint i = 0; i <= RDdef; ++i, lc += ldc)
  {
    ldbl c = exp(lc);
    fout << c << " " << ncctolp(c, df, lamda) << "\n";
  }
}

void plotclass::plot_ncfden(ostream& fout, const option& opt) const
// 非心F分布密度関数出力
{
  constexpr int DFdef = 1;      // 省略時自由度
  constexpr ldbl Ldef = 0.0L;   // 省略時λ
  constexpr ldbl FMAXdef = 10.0L; // 省略時F最大値
  constexpr ldbl EPSf = 0.001L;   // F微小値
  int df1 = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : DFdef;
  int df2 = (opt.nipara() > 1 && opt.ipara(1) >= 0) ? opt.ipara(1) : DFdef;
  ldbl lamda = (opt.nfpara() > 0 && opt.fpara(0) >= 0.0) ? (ldbl)(opt.fpara(0)) : Ldef;
  ldbl fmin = 0.0, fmax = FMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) fmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) fmax = (ldbl)(opt.rmax());
  ldbl lfmin = log((fmin > 0.0) ? fmin : fmax * EPSf), lfmax = log(fmax);
  ldbl ldf = (lfmax - lfmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%非心F分布密度関数:df1=" << df1 << " df2=" << df2 << " 非心度=" << lamda << "\n";
  if (fmin == 0.0L) fout << "0 0\n";
  ldbl lf = lfmin;
  for (uint i = 0; i <= RDdef; ++i, lf += ldf)
  {
    ldbl f = exp(lf);
    fout << f << " " << ncfden(f, df1, df2, lamda) << "\n";
  }
}

void plotclass::plot_ncfdis(ostream& fout, const option& opt) const
// 非心F分布関数出力
{
  constexpr int DFdef = 1;      // 省略時自由度
  constexpr ldbl Ldef = 0.0L;   // 省略時λ
  constexpr ldbl FMAXdef = 10.0L; // 省略時F最大値
  constexpr ldbl EPSf = 0.001L;   // F微小値
  int df1 = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : DFdef;
  int df2 = (opt.nipara() > 1 && opt.ipara(1) >= 0) ? opt.ipara(1) : DFdef;
  ldbl lamda = (opt.nfpara() > 0 && opt.fpara(0) >= 0.0) ? (ldbl)(opt.fpara(0)) : Ldef;
  ldbl fmin = 0.0, fmax = FMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) fmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) fmax = (ldbl)(opt.rmax());
  ldbl lfmin = log((fmin > 0.0) ? fmin : fmax * EPSf), lfmax = log(fmax);
  ldbl ldf = (lfmax - lfmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%非心F分布関数:df1=" << df1 << " df2=" << df2 << " 非心度=" << lamda << "\n";
  if (fmin == 0.0L) fout << "0 0\n";
  ldbl lf = lfmin;
  for (uint i = 0; i <= RDdef; ++i, lf += ldf)
  {
    ldbl f = exp(lf);
    fout << f << " " << ncftolp(f, df1, df2, lamda) << "\n";
  }
}

void plotclass::plot_nctden(ostream& fout, const option& opt) const
// 非心t分布密度関数出力
{
  constexpr int DFdef = 1;      // 省略時自由度
  constexpr ldbl Ldef = 0.0L;   // 省略時λ
  constexpr ldbl TMAXdef = 5.0L; // 省略時t最大値
  int df = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : DFdef;
  ldbl lamda = (opt.nfpara() > 0) ? (ldbl)(opt.fpara(0)) : Ldef;
  ldbl tmin = -TMAXdef, tmax = TMAXdef;
  if (opt.rmin() != BLANKf) tmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) tmax = (ldbl)(opt.rmax());
  ldbl dt = (tmax - tmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%非心t分布密度関数:df=" << df << " 非心度=" << lamda << "\n";
  ldbl t = tmin;
  const ldbl Teps = dt * 0.01L; // t計算誤差
  for (uint i = 0; i <= RDdef; ++i, t += dt)
  {
    if (fabs(t) < Teps) t = 0.0L;
    fout << t << " " << nctden(t, df, lamda) << "\n";
  }
}

void plotclass::plot_nctdis(ostream& fout, const option& opt) const
// 非心t分布関数出力
{
  constexpr int DFdef = 1;      // 省略時自由度
  constexpr ldbl Ldef = 0.0L;   // 省略時λ
  constexpr ldbl TMAXdef = 5.0L; // 省略時t最大値
  int df = (opt.nipara() > 0 && opt.ipara(0) > 0) ? opt.ipara(0) : DFdef;
  ldbl lamda = (opt.nfpara() > 0) ? (ldbl)(opt.fpara(0)) : Ldef;
  ldbl tmin = -TMAXdef, tmax = TMAXdef;
  if (opt.rmin() != BLANKf) tmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) tmax = (ldbl)(opt.rmax());
  ldbl dt = (tmax - tmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%非心t分布関数:df=" << df << " 非心度=" << lamda << "\n";
  ldbl t = tmin;
  ldbl teps = dt * 0.01L;       // t計算誤差
  for (uint i = 0; i <= RDdef; ++i, t += dt)
  {
    if (fabs(t) < teps) t = 0.0L;
    fout << t << " " << ncttolp(t, df, lamda) << "\n";
  }
}

void plotclass::plot_logden(ostream& fout, const option& opt) const
// ロジスティック分布密度関数出力
{
  constexpr ldbl YMdef = 0.0L;  // 省略時ym
  constexpr ldbl Mdef = 1.0L;   // 省略時M
  constexpr ldbl Cdef = 1.0L;   // 省略時c
  constexpr ldbl PROy0 = 0.5L;  // 省略時y0の割合
  constexpr ldbl Yeps = 0.01L;  // yの微小値
  ldbl ym = (opt.nfpara() > 0) ? (ldbl)(opt.fpara(0)) : YMdef;
  ldbl m = (opt.nfpara() > 1 && ym < (ldbl)(opt.fpara(1))) ? (ldbl)(opt.fpara(1)) : Mdef;
  ldbl c = (opt.nfpara() > 2) ? (ldbl)(opt.fpara(2)) : Cdef;
  ldbl y0 = (opt.nfpara() > 3 && ym < (ldbl)(opt.fpara(3)) && (ldbl)(opt.fpara(3)) < ym + m) ? (ldbl)(opt.fpara(3)) : ym + m * PROy0;
  ldbl t0 = (opt.nfpara() > 4) ? (ldbl)(opt.fpara(4)) : 0.0L;
  ldbl y0m = y0 - ym;
  ldbl ly0 = log((m - y0m) / y0m), y = m * Yeps;
  ldbl tmin = (ly0 - log((m - y) / y)) / c;
  y = m * (1.0L - Yeps);
  ldbl tmax = (ly0 - log((m - y) / y)) / c;
  if (opt.rmin() != BLANKf) tmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) tmax = (ldbl)(opt.rmax());
  ldbl dt = (tmax - tmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%ロジスティック分布密度関数:ym=" << ym << " M=" << m << " c=" << c << " y0=" << y0 << " t0=" << t0 << "\n";
  ldbl t = tmin;
  ldbl teps = dt * 0.01L;       // t計算誤差
  for (uint i = 0; i <= RDdef; ++i, t += dt)
  {
    if (fabs(t) < teps) t = 0.0L;
    fout << t << " " << logden(t - t0, m, c, y0m) << "\n";
  }
}

void plotclass::plot_logdis(ostream& fout, const option& opt) const
// ロジスティック分布関数出力(ロジスティック曲線)
{
  constexpr ldbl YMdef = 0.0L;  // 省略時ym
  constexpr ldbl Mdef = 1.0L;   // 省略時M
  constexpr ldbl Cdef = 1.0L;   // 省略時c
  constexpr ldbl PROy0 = 0.5L;  // 省略時y0の割合
  constexpr ldbl Yeps = 0.01L;  // yの微小値
  ldbl ym = (opt.nfpara() > 0) ? opt.fpara(0) : YMdef;
  ldbl m = (opt.nfpara() > 1 && ym < (ldbl)(opt.fpara(1))) ? (ldbl)(opt.fpara(1)) : Mdef;
  ldbl c = (opt.nfpara() > 2) ? (ldbl)(opt.fpara(2)) : Cdef;
  ldbl y0 = (opt.nfpara() > 3 && ym < (ldbl)(opt.fpara(3)) && (ldbl)(opt.fpara(3)) < ym + m) ? (ldbl)(opt.fpara(3)) : ym + m * PROy0;
  ldbl t0 = (opt.nfpara() > 4) ? (ldbl)(opt.fpara(4)) : 0.0L;
  ldbl y0m = y0 - ym;
  ldbl ly0 = log((m - y0m) / y0m), y = m * Yeps;
  ldbl tmin = (ly0 - log((m - y) / y)) / c;
  y = m * (1.0L - Yeps);
  ldbl tmax = (ly0 - log((m - y) / y)) / c;
  if (opt.rmin() != BLANKf) tmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) tmax = (ldbl)(opt.rmax());
  ldbl dt = (tmax - tmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%ロジスティック分布関数:ym=" << ym << " M=" << m << " c=" << c << " y0=" << y0 << " t0=" << t0 << "\n";
  ldbl t = tmin;
  ldbl teps = dt * 0.01L;       // t計算誤差
  for (uint i = 0; i <= RDdef; ++i, t += dt)
  {
    if (fabs(t) < teps) t = 0.0L;
    fout << t << " " << ym + logdis(t - t0, m, c, y0m) << "\n";
  }
}

void plotclass::plot_nden(ostream& fout, const option& opt) const
// 正規分布密度関数出力
{
  constexpr int Ndef = 1;       // 省略時例数
  constexpr ldbl Mdef = 0.0L;   // 省略時mean
  constexpr ldbl SDdef = 1.0L;  // 省略時SD
  int n = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : Ndef;
  ldbl mean = (opt.nfpara() > 0) ? (ldbl)(opt.fpara(0)) : Mdef;
  ldbl sd = (opt.nfpara() > 1 && opt.fpara(1) > 0.0) ? (ldbl)(opt.fpara(1)) : SDdef;
  ldbl cf = opt.ct() ? 0.5L : 0.0L;
  ldbl xmin = mean - sd * 4.0L, xmax = mean + sd * 4.0L;
  if (opt.rmin() != BLANKf) xmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) xmax = (ldbl)(opt.rmax());
  ldbl dx = (xmax - xmin) / ldbl(RDdef), hx = 1.0L;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%正規分布関数:例数=" << n << " 平均=" << mean << " 標準偏差=" << sd;
  if (opt.ct()) fout << " 連続修正有";
  fout << "\n";
  ldbl x = xmin;
  const ldbl Xeps = dx * 0.01L; // x計算誤差
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    if (fabs(x) < Xeps) x = 0.0L;
    ldbl dev = x - mean;
    ldbl devct = dev - sgn(dev) * cf;
    fout << x << " " << (ldbl)n * snden(devct / sd) << "\n";
  }
}

void plotclass::plot_ndis(ostream& fout, const option& opt) const
// 正規分布関数出力
{
  constexpr int Ndef = 1;       // 省略時例数
  constexpr ldbl Mdef = 0.0L;   // 省略時mean
  constexpr ldbl SDdef = 1.0L;  // 省略時SD
  int n = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : Ndef;
  ldbl mean = (opt.nfpara() > 0) ? (ldbl)(opt.fpara(0)) : Mdef;
  ldbl sd = (opt.nfpara() > 1 && opt.fpara(1) > 0.0) ? (ldbl)(opt.fpara(1)) : SDdef;
  ldbl cf = opt.ct() ? 0.5L : 0.0L;
  ldbl xmin = mean - sd * 4.0L, xmax = mean + sd * 4.0L;
  if (opt.rmin() != BLANKf) xmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) xmax = (ldbl)(opt.rmax());
  ldbl dx = (xmax - xmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%正規分布関数:例数=" << n << " 平均=" << mean << " 標準偏差=" << sd;
  if (opt.ct()) fout << " 連続修正有";
  fout << "\n";
  ldbl x = xmin;
  ldbl xeps = dx * 0.01L;       // x計算誤差
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    if (fabs(x) < xeps) x = 0.0L;
    ldbl dev = x - mean;
    fout << x << " " << (ldbl)n * sndis(dev / sd) << "\n";
  }
}

void plotclass::plot_npval(ostream& fout, const option& opt) const
// 正規分布両側確率出力
{
  constexpr ldbl Mdef = 0.0L;   // 省略時mean
  constexpr ldbl SDdef = 1.0L;  // 省略時SD
  ldbl mean = (opt.nfpara() > 0) ? (ldbl)(opt.fpara(0)) : Mdef;
  ldbl sd = (opt.nfpara() > 1 && opt.fpara(1) > 0.0) ? (ldbl)(opt.fpara(1)) : SDdef;
  ldbl cf = opt.ct() ? 0.5L : 0.0L;
  ldbl xmin = 0.0L, xmax = mean + sd * 4.0L;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) xmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) xmax = (ldbl)(opt.rmax());
  ldbl dx = (xmax - xmin) / ldbl(RDdef);
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%正規分布両側確率:平均=" << mean << " 標準偏差=" << sd;
  if (opt.ct()) fout << " 連続修正有";
  fout << "\n";
  ldbl x = xmin;
  const ldbl Xeps = dx * 0.01L; // x計算誤差
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    if (fabs(x) < Xeps) x = 0.0L;
    ldbl dev = x - mean;
    ldbl pr = ztop((dev - cf) / sd);
    fout << x << " " << pr << "\n";
  }
}

void plotclass::plot_npinv(ostream& fout, const option& opt) const
// 正規分布両側確率逆関数出力
{
  constexpr ldbl Mdef = 0.0L;    // 省略時mean
  constexpr ldbl SDdef = 1.0L;   // 省略時SD
  constexpr ldbl PMINdef = 1.0e-5L; // 省略時最小p値
  constexpr ldbl PMAXdef = 1.0L;    // 省略時最大p値
  constexpr ldbl Xeps = LDBL_EPSILON * 10.0L; // x計算誤差
  ldbl mean = (opt.nfpara() > 0) ? (ldbl)(opt.fpara(0)) : Mdef;
  ldbl sd = (opt.nfpara() > 1 && opt.fpara(1) > 0.0) ? (ldbl)(opt.fpara(1)) : SDdef;
  ldbl cf = opt.ct() ? 0.5L : 0.0L;
  ldbl pmin = PMINdef, pmax = PMAXdef;
  if (opt.rmin() != BLANKf && 0.0 < opt.rmin() && opt.rmin() < 1.0) pmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax() && opt.rmax() <= 1.0) pmax = (ldbl)(opt.rmax());
  ldbl lpmin = log(pmin), lpmax = log(pmax);
  ldbl ldp = (lpmax - lpmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%正規分布両側確率逆関数:平均=" << mean << " 標準偏差=" << sd;
  if (opt.ct()) fout << " 連続修正有";
  fout << "\n";
  ldbl lp = lpmax;
  for (uint i = 0; i <= RDdef; ++i, lp -= ldp)
  {
    ldbl pr = exp(lp);
    ldbl x = ptoz(pr) * sd + mean + cf;
    if (fabs(x) < Xeps) x = 0.0L;
    fout << x << " " << pr << "\n";
  }
}

void plotclass::plot_perdis(ostream& fout, const option& opt) const
// 周期曲線出力
{
  constexpr ldbl T0def = 24.0L; // 省略時基本周期
  constexpr uint RTdef = 48;    // 省略時曲線の分割数
  ldbl t0 = (opt.nfpara() > 0) ? (ldbl)(opt.fpara(0)) : T0def;
  size_t m = (opt.nipara() > 0) ? opt.nipara() : 1; // 周期成分数
  vecint pcvec(m);
  vector<pairldbl> aapvec(m + 1);
  ldbl radf = 2.0L * (ldbl)M_PI / t0;
  aapvec[0].first = (opt.nfpara() > 1) ? (ldbl)(opt.fpara(1)) : 1.0L;
  aapvec[0].second = 0.0L;
  for (size_t k = 1; k <= m; ++k)
  {
    pcvec[k - 1] = (opt.nipara() >= k) ? opt.ipara(k - 1) : k;
    size_t i = k * 2;
    aapvec[k].first = (opt.nfpara() > i) ? (ldbl)(opt.fpara(i)) : hypot(0.1L, 0.1L);
    ldbl th = (opt.nfpara() > i + 1) ? (ldbl)(opt.fpara(i + 1)) : t0 / 8.0L;
    aapvec[k].second = th * radf;
  }
  pfun<ldbl> pf(t0, pcvec, aapvec);
  ldbl tmin = 0.0, tmax = t0;
  if (opt.rmin() != BLANKf) tmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) tmax = (ldbl)(opt.rmax());
  ldbl dt = (tmax - tmin) / (ldbl)RTdef;
  ldbl degf = 180.0L / (ldbl)M_PI;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%周期曲線:T=" << pf.t(0) << " y=" << pf.amp(0);
  for (size_t k = 1; k <= pf.m(); ++k)
  {
    fout << showpos << pf.amp(k) << noshowpos << "・cos(" << pf.w(k) * degf << "(t " << showpos << -pf.acr(k) / radf << "))";
  }
  fout << noshowpos << "\n";
  ldbl t = tmin;
  for (uint i = 0; i <= RTdef; ++i, t += dt)
  { fout << t << " " << pf.fun(t) << "\n"; }
}

void plotclass::plot_plogdis(ostream& fout, const option& opt) const
// ロジスティック関数の陽性率曲線出力
{
  constexpr ldbl YMdef = 0.0L;  // 省略時ym
  constexpr ldbl Mdef = 1.0L;   // 省略時M
  constexpr ldbl Cdef = 1.0L;   // 省略時c
  constexpr ldbl PROy0 = 0.5L;  // 省略時y0の割合
  constexpr ldbl COdef = 0.5L;  // 省略時陽性の境界値
  constexpr ldbl VRdef = 0.1L;  // 省略時残差分散
  constexpr ldbl Yeps = 0.01L;  // yの微小値
  constexpr int Ndef = 100;     // 省略時例数
  ldbl ym = (opt.nfpara() > 0) ? (ldbl)(opt.fpara(0)) : YMdef;
  ldbl m = (opt.nfpara() > 1 && ym < (ldbl)(opt.fpara(1))) ? (ldbl)(opt.fpara(1)) : Mdef;
  ldbl c = (opt.nfpara() > 2) ? (ldbl)(opt.fpara(2)) : Cdef;
  ldbl y0 = (opt.nfpara() > 3 && ym < (ldbl)(opt.fpara(3)) && (ldbl)(opt.fpara(3)) < ym + m) ? (ldbl)(opt.fpara(3)) : ym + m * PROy0;
  ldbl co = (opt.nfpara() > 4) ? (ldbl)(opt.fpara(4)) : COdef;
  ldbl vr = (opt.nfpara() > 5) ? (ldbl)(opt.fpara(5)) : VRdef;
  int n = (opt.nipara() > 0) ? (ldbl)(opt.ipara(0)) : Ndef;
  ldbl yh = y0 - ym, ly0 = log((m - yh) / yh);
  yh = co - ym;
  ldbl lco = log((m - yh) / yh);
  ldbl sdl = sqrt(vr * (1.0 + 1.0 / ldbl(n)));
  int dfr = n - 2;
  ldbl y = m * Yeps;
  ldbl tmin = (ly0 - log((m - y) / y)) / c;
  y = m * (1.0 - Yeps);
  ldbl tmax = (ly0 - log((m - y) / y)) / c;
  if (opt.rmin() != BLANKf) tmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) tmax = (ldbl)(opt.rmax());
  ldbl dt = (tmax - tmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%ロジスティック関数の陽性率曲線:ym=" << ym << " M=" << m << " c=" << c << " y0=" << y0 << " CO=" << co << " Vr=" << vr << " n=" << n << "\n";
  ldbl t = tmin, ly = 0.0;
  for (uint i = 0; i <= RDdef; ++i, t += dt)
  {
    ly = ly0 - c * t;
    ldbl ldif = ly - lco;
    ldbl lned = ldif / sdl;
    ldbl pr = 0.5L * ttop(lned, dfr);
    if (ldif < 0.0L) pr = 1.0L - pr;
    fout << t << " " << pr << "\n";
  }
}

void plotclass::plot_poiden(ostream& fout, const option& opt) const
// Poisson分布密度関数出力
{
  constexpr int Rdef = 1;       // 省略時r
  int r = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : Rdef; // 例数
  double ldef = (double)r;      // 省略時λ
  double lamda = (opt.nfpara() > 0 && 0.0 < opt.fpara(0)) ? opt.fpara(0) : ldef;
  int imin = 0, imax = r;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin() && opt.rmin() < (double)r) imin = (int)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax() && opt.rmax() <= (double)r) imax = (int)(opt.rmax());
  fout << "-s1 -drectangle_solid_fill -x1 -y3 -w80\n";
  fout << "%Poisson分布密度関数:λ=" << lamda << " r=" << r << "\n";
  for (int i = imin; i <= imax; ++i)
  {
    double x = (double)i;
    fout << x - 0.1 << " 0 " << x + 0.1 << " " << poiden(i, lamda) << "\n";
  }
}

void plotclass::plot_poidis(ostream& fout, const option& opt) const
// Poisson分布関数出力
{
  constexpr int Rdef = 1;       // 省略時r
  int r = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : Rdef; // 例数
  if (r < 1) r = Rdef;
  double ldef = (double)r;   // 省略時λ
  double lamda = (opt.nfpara() > 0 && 0.0 < opt.fpara(0)) ? opt.fpara(0) : ldef;
  int imin = 0, imax = r;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin() && opt.rmin() < (double)r) imin = (int)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax() && opt.rmax() <= (double)r) imax = (int)(opt.rmax());
  fout << "-s1 -dline -w80\n";
  fout << "%Poisson分布関数:λ=" << lamda << " r=" << r << "\n";
  for (int i = imin; i <= imax; ++i)
  { fout << i << " " << poidis(i, lamda) << "\n"; }
}

void plotclass::plot_poly(ostream& fout, const option& opt) const
// 多項関数出力
{
  constexpr int Ndef = 1;       // 省略時次数
  int n = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : Ndef; // 次数
  if (n < 1) n = Ndef;
  double xmin = 0.0, xmax = 100.0;
  if (opt.rmin() != BLANKf) xmin = opt.rmin();
  if (opt.rmax() != BLANKf && xmin < opt.rmax()) xmax = opt.rmax();
  double dx = (xmax - xmin) / (double)RDdef;
  fout << "-s1 -dline -w80\n";
  fout << "%多項関数:次数=" << n << "\n";
  double bi = (opt.nfpara() > 0) ? opt.fpara(0) : 0.0;
  fout << "%b0=" << bi << "\n";
  for (int i = 1; i <= n; ++i)
  {
    bi = (opt.nfpara() > size_t(i)) ? opt.fpara(i) : 1.0;
    fout << "%b" << i << "=" << bi << "\n";
  }
  double x = xmin;
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  { fout << x << " " << cal_poly(opt, n, x) << "\n"; }
}

double plotclass::cal_poly(const option& opt, int n, double x) const
// 多項関数値を計算して返す
{
  double y = (opt.nfpara() > 0) ? opt.fpara(0) : 0.0;
  double xi = x;
  for (int i = 1; i <= n; ++i)
  {
    double bi = (opt.nfpara() > size_t(i)) ? opt.fpara(i) : 1.0;
    y += bi * xi;
    xi *= x;
  }
  return y;                     // 多項関数値を返す
}

void plotclass::plot_npv(ostream& fout, const option& opt) const // NPV曲線出力
{
  constexpr double SNdef = 0.5; // 省略時感度
  constexpr double SPdef = 0.5; // 省略時特異度
  constexpr double PIMINdef = 0.01; // 省略時π最小値
  constexpr double PIMAXdef = 0.99; // 省略時π最大値
  double sn = (opt.nfpara() > 0 && 0.0 < opt.fpara(0) && opt.fpara(0) < 1.0) ? opt.fpara(0) : SNdef;
  double sp = (opt.nfpara() > 1 && 0.0 < opt.fpara(1) && opt.fpara(1) < 1.0) ? opt.fpara(1) : SPdef;
  double pimin = PIMINdef, pimax = PIMAXdef;
  if (opt.rmin() != BLANKf && 0.0 < opt.rmin() && opt.rmin() < 1.0) pimin = opt.rmin();
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax() && opt.rmax() < 1.0) pimax = opt.rmax();
  double dpi = (pimax - pimin) / (double)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%NPV曲線:感度=" << sn << " 特異度=" << sp << "\n";
  fout << "0 1\n";
  double pi = pimin;
  for (uint i = 0; i <= RDdef; ++i, pi += dpi)
  {
    double npv = 1.0 / (1.0 + pi / (1.0 - pi) * (1.0 - sn) / sp);
    fout << pi << " " << npv << "\n";
  }
  fout << "1 0\n";
}

void plotclass::plot_ppv(ostream& fout, const option& opt) const // PPV曲線出力
{
  constexpr double SNdef = 0.5; // 省略時感度
  constexpr double SPdef = 0.5; // 省略時特異度
  constexpr double PIMINdef = 0.01; // 省略時π最小値
  constexpr double PIMAXdef = 0.99; // 省略時π最大値
  double sn = (opt.nfpara() > 0 && 0.0 < opt.fpara(0) && opt.fpara(0) < 1.0) ? opt.fpara(0) : SNdef;
  double sp = (opt.nfpara() > 1 && 0.0 < opt.fpara(1) && opt.fpara(1) < 1.0) ? opt.fpara(1) : SPdef;
  double pimin = PIMINdef, pimax = PIMAXdef;
  if (opt.rmin() != BLANKf && 0.0 < opt.rmin() && opt.rmin() < 1.0) pimin = opt.rmin();
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax() && opt.rmax() < 1.0) pimax = opt.rmax();
  double dpi = (pimax - pimin) / (double)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%PPV曲線:感度=" << sn << " 特異度=" << sp << "\n";
  fout << "0 0\n";
  double pi = pimin;
  for (uint i = 0; i <= RDdef; ++i, pi += dpi)
  {
    double ppv = 1.0 / (1.0 + (1.0 - pi) / pi * (1.0 - sp) / sn);
    fout << pi << " " << ppv << "\n";
  }
  fout << "1 1\n";
}

void plotclass::plot_proden(ostream& fout, const option& opt) const
// プロビット分布密度関数出力(正規分布)
{
  constexpr ldbl ZMINdef = -3.0L; // 省略時最小z値
  constexpr ldbl ZMAXdef = 3.0L;  // 省略時最大z値
  ldbl zmin = ZMINdef, zmax = ZMAXdef;
  if (opt.rmin() != BLANKf) zmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) zmax = (ldbl)(opt.rmax());
  ldbl dz = (zmax - zmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%プロビット分布密度関数:1SD=1\n";
  ldbl z = zmin;
  for (uint i = 0; i <= RDdef; ++i, z += dz)
  {
    ldbl x = z + 5.0L;
    fout << x << " " << snden(z) << "\n";
  }
}

void plotclass::plot_prodis(ostream& fout, const option& opt) const
// プロビット分布関数出力(プロビット曲線)
{
  constexpr ldbl ZMINdef = -3.0L; // 省略時最小z値
  constexpr ldbl ZMAXdef = 3.0L;  // 省略時最大z値
  ldbl zmin = ZMINdef, zmax = ZMAXdef;
  if (opt.rmin() != BLANKf) zmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) zmax = (ldbl)(opt.rmax());
  ldbl dz = (zmax - zmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%プロビット曲線:1SD=1\n";
  ldbl z = zmin;
  for (uint i = 0; i <= RDdef; ++i, z += dz)
  {
    ldbl x = z + 5.0L;
    fout << x << " " << sndis(z) << "\n";
  }
}

void plotclass::plot_snden(ostream& fout, const option& opt) const
// 標準正規分布密度関数出力
{
  constexpr ldbl ZMINdef = -4.0L; // 省略時最小z値
  constexpr ldbl ZMAXdef = 4.0L;  // 省略時最大z値
  ldbl zmin = ZMINdef, zmax = ZMAXdef;
  if (opt.rmin() != BLANKf) zmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) zmax = (ldbl)(opt.rmax());
  ldbl dz = (zmax - zmin) / (ldbl)RDdef;
  ldbl cf = opt.ct() ? 0.5L : 0.0L;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%標準正規分布密度関数:平均=0 標準偏差=1";
  if (opt.ct()) fout << " 連続修正有";
  fout << "\n";
  ldbl z = zmin;
  ldbl zeps = dz * 0.01L;       // z計算誤差
  for (uint i = 0; i <= RDdef; ++i, z += dz)
  {
    if (fabs(z) < zeps) z = 0.0L;
    fout << z << " " << snden(z - sgn(z) * cf) << "\n";
  }
}

void plotclass::plot_sndis(ostream& fout, const option& opt) const
// 標準正規分布関数出力
{
  constexpr ldbl ZMINdef = -4.0L; // 省略時最小z値
  constexpr ldbl ZMAXdef = 4.0L;  // 省略時最大z値
  ldbl zmin = ZMINdef, zmax = ZMAXdef;
  if (opt.rmin() != BLANKf) zmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) zmax = (ldbl)(opt.rmax());
  ldbl dz = (zmax - zmin) / (ldbl)RDdef;
  ldbl cf = opt.ct() ? 0.5L : 0.0L;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%標準正規分布関数:平均=0 標準偏差=1";
  if (opt.ct()) fout << " 連続修正有";
  fout << "\n";
  ldbl z = zmin;
  ldbl zeps = dz * 0.01L;       // z計算誤差
  for (uint i = 0; i <= RDdef; ++i, z += dz)
  {
    if (fabs(z) < zeps) z = 0.0L;
    fout << z << " " << sndis(z - sgn(z) * cf) << "\n";
  }
}

void plotclass::plot_snpval(ostream& fout, const option& opt) const
// 標準正規分布両側確率出力
{
  constexpr ldbl ZMAXdef = 4.0L; // 省略時最大z値
  ldbl zmin = 0.0, zmax = ZMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) zmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) zmax = (ldbl)(opt.rmax());
  ldbl dz = (zmax - zmin) / (ldbl)RDdef;
  ldbl cf = opt.ct() ? 0.5L : 0.0L;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%標準正規分布両側確率:平均=0 標準偏差=1";
  if (opt.ct()) fout << " 連続修正有";
  fout << "\n";
  ldbl z = zmin;
  const ldbl Zeps = dz * 0.01L; // z計算誤差
  for (uint i = 0; i <= RDdef; ++i, z += dz)
  {
    if (fabs(z) < Zeps) z = 0.0L;
    fout << z << " " << ztop(z - sgn(z) * cf) << "\n";
  }
}

void plotclass::plot_snpinv(ostream& fout, const option& opt) const
// 標準正規分布両側確率逆関数出力
{
  constexpr ldbl PMINdef = 1.0e-5L; // 省略時最小p値
  constexpr ldbl PMAXdef = 1.0L;    // 省略時最大p値
  constexpr ldbl Zeps = LDBL_EPSILON * 4.0L; // z計算誤差
  ldbl pmin = PMINdef, pmax = PMAXdef;
  if (opt.rmin() != BLANKf && 0.0 < opt.rmin() && opt.rmin() < 1.0) pmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax() && opt.rmax() <= 1.0) pmax = (ldbl)(opt.rmax());
  ldbl lpmin = log(pmin), lpmax = log(pmax);
  ldbl ldp = (lpmax - lpmin) / (ldbl)RDdef;
  ldbl cf = opt.ct() ? 0.5L : 0.0L;
  ldbl z = 0.0L;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%標準正規分布両側確率逆関数:平均=0 標準偏差=1";
  if (opt.ct()) fout << " 連続修正有";
  fout << "\n";
  ldbl lp = lpmax;
  for (uint i = 0; i <= RDdef; ++i, lp -= ldp)
  {
    ldbl pr = exp(lp);
    ldbl z = snpinv(pr) + cf;
    if (fabs(z) < Zeps) z = 0.0L;
    fout << z << " " << pr << "\n";
  }
}

void plotclass::plot_lnden(ostream& fout, const option& opt) const
// 標準対数正規密度関数出力
{
  constexpr ldbl ZMINdef = -4.0L; // 省略時最小z値
  constexpr ldbl ZMAXdef = 4.0L;  // 省略時最大z値
  ldbl zmin = ZMINdef, zmax = ZMAXdef;
  if (opt.rmin() != BLANKf) zmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) zmax = (ldbl)(opt.rmax());
  ldbl dz = (zmax - zmin) / (ldbl)RDdef;
  ldbl cf = opt.ct() ? 0.5L : 0.0L;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%標準対数正規分布密度関数:ln(x)平均=0 ln(x)標準偏差=1";
  if (opt.ct()) fout << " 連続修正有";
  fout << "\n";
  ldbl z = zmin;
  const ldbl Zeps = dz * 0.01L; // z計算誤差
  for (uint i = 0; i <= RDdef; ++i, z += dz)
  {
    if (fabs(z) < Zeps) z = 0.0L;
    fout << exp(z) << " " << snden(z - sgn(z) * cf) << "\n";
  }
}

void plotclass::plot_lndis(ostream& fout, const option& opt) const
// 標準対数正規分布関数出力
{
  constexpr ldbl ZMINdef = -4.0L; // 省略時最小z値
  constexpr ldbl ZMAXdef = 4.0L;  // 省略時最大z値
  ldbl zmin = ZMINdef, zmax = ZMAXdef;
  if (opt.rmin() != BLANKf) zmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) zmax = (ldbl)(opt.rmax());
  ldbl dz = (zmax - zmin) / (ldbl)RDdef;
  ldbl cf = opt.ct() ? 0.5L : 0.0L;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%標準対数正規分布関数:ln(x)平均=0 ln(x)標準偏差=1";
  if (opt.ct()) fout << " 連続修正有";
  fout << "\n";
  ldbl z = zmin;
  const ldbl Zeps = dz * 0.01L; // z計算誤差
  for (uint i = 0; i <= RDdef; ++i, z += dz)
  {
    if (fabs(z) < Zeps) z = 0.0L;
    fout << exp(z) << " " << sndis(z - sgn(z) * cf) << "\n";
  }
}

void plotclass::plot_lnhaz(ostream& fout, const option& opt) const
// 標準対数正規分布ハザード関数出力
{
  constexpr ldbl ZMINdef = -4.0L; // 省略時最小z値
  constexpr ldbl ZMAXdef = 4.0L;  // 省略時最大z値
  ldbl zmin = ZMINdef, zmax = ZMAXdef;
  if (opt.rmin() != BLANKf) zmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) zmax = (ldbl)(opt.rmax());
  ldbl dz = (zmax - zmin) / (ldbl)RDdef;
  ldbl cf = opt.ct() ? 0.5L : 0.0L;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%標準対数正規分布ハザード関数:ln(x)平均=0 ln(x)標準偏差=1";
  if (opt.ct()) fout << " 連続修正有";
  fout << "\n";
  ldbl z = zmin;
  const ldbl Zeps = dz * 0.01L; // z計算誤差
  for (uint i = 0; i <= RDdef; ++i, z += dz)
  {
    if (fabs(z) < Zeps) z = 0.0L;
    fout << exp(z) << " " << snhaz(z - sgn(z) * cf) << "\n";
  }
}

void plotclass::plot_lnpval(ostream& fout, const option& opt) const
// 標準対数正規分布両側確率出力
{
  constexpr ldbl ZMAXdef = 4.0l; // 省略時最大z値
  ldbl zmin = 0.0L, zmax = ZMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) zmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) zmax = (ldbl)(opt.rmax());
  ldbl dz = (zmax - zmin) / (ldbl)RDdef;
  ldbl cf = opt.ct() ? 0.5L : 0.0L;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%標準対数正規分布両側確率:ln(x)平均=0 ln(x)標準偏差=1";
  if (opt.ct()) fout << " 連続修正有";
  fout << "\n";
  ldbl z = zmin;
  const ldbl Zeps = dz * 0.01L; // z計算誤差
  for (uint i = 0; i <= RDdef; ++i, z += dz)
  {
    if (fabs(z) < Zeps) z = 0.0L;
    fout << exp(z) << " " << ztop(z - sgn(z) * cf) << "\n";
  }
}

void plotclass::plot_lnsur(ostream& fout, const option& opt) const
// 標準対数正規分布生存関数出力
{
  constexpr ldbl ZMINdef = -4.0L; // 省略時z値
  constexpr ldbl ZMAXdef = 4.0L;  // 省略時z値
  ldbl zmin = ZMINdef, zmax = ZMAXdef;
  if (opt.rmin() != BLANKf) zmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) zmax = (ldbl)(opt.rmax());
  ldbl dz = (zmax - zmin) / (ldbl)RDdef;
  ldbl cf = opt.ct() ? 0.5L : 0.0L;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%標準対数正規分布生存関数:ln(x)平均=0 ln(x)標準偏差=1";
  if (opt.ct()) fout << " 連続修正有";
  fout << "\n";
  ldbl z = zmin;
  const ldbl Zeps = dz * 0.01L; // z計算誤差
  for (uint i = 0; i <= RDdef; ++i, z += dz)
  {
    if (fabs(z) < Zeps) z = 0.0L;
    fout << exp(z) << " " << snsur(z - sgn(z) * cf) << "\n";
  }
}

void plotclass::plot_tden(ostream& fout, const option& opt) const
// t分布密度関数出力
{
  constexpr int DFdef = 1;      // 省略時自由度
  constexpr int Ndef = 1;       // 省略時例数
  constexpr ldbl TMINdef = -5.0L; // 省略時t最小値
  constexpr ldbl TMAXdef = 5.0L;  // 省略時t最大値
  constexpr ldbl Mdef = 0.0L;   // 省略時mean
  constexpr ldbl SDdef = 1.0L;  // 省略時SD
  int df = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : DFdef;
  int n = (opt.nipara() > 1 && 0 < opt.ipara(1)) ? opt.ipara(1) : Ndef;
  ldbl mean = (opt.nfpara() > 0) ? (ldbl)(opt.fpara(0)) : Mdef;
  ldbl sd = (opt.nfpara() > 1 && opt.fpara(1) > 0.0) ? (ldbl)(opt.fpara(1)) : SDdef;
  ldbl xmin = TMINdef, xmax = TMAXdef;
  if (opt.rmin() != BLANKf) xmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) xmax = (ldbl)(opt.rmax());
  ldbl dx = (xmax - xmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%t分布密度関数:自由度=" << df << " 例数=" << n << " 平均値=" << mean << " 標準偏差=" << sd << "\n";
  ldbl x = xmin;
  const ldbl Xeps = dx * 0.01L; // x計算誤差
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    if (fabs(x) < Xeps) x = 0.0L;
    ldbl t = (x - mean) / sd;
    fout << x << " " << (ldbl)n * tden(t, df) << "\n";
  }
}

void plotclass::plot_tdis(ostream& fout, const option& opt) const
// t分布関数出力
{
  constexpr int DFdef = 1;      // 省略時自由度
  constexpr int Ndef = 1;       // 省略時例数
  constexpr ldbl TMINdef = -5.0L; // 省略時t最小値
  constexpr ldbl TMAXdef = 5.0L;  // 省略時t最大値
  constexpr ldbl Mdef = 0.0L;     // 省略時mean
  constexpr ldbl SDdef = 1.0L;    // 省略時SD
  int df = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : DFdef;
  int n = (opt.nipara() > 1 && 0 < opt.ipara(1)) ? opt.ipara(1) : Ndef;
  ldbl mean = (opt.nfpara() > 0) ? (ldbl)(opt.fpara(0)) : Mdef;
  ldbl sd = (opt.nfpara() > 1 && opt.fpara(1) > 0.0) ? (ldbl)(opt.fpara(1)) : SDdef;
  ldbl xmin = TMINdef, xmax = TMAXdef;
  if (opt.rmin() != BLANKf) xmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) xmax = (ldbl)(opt.rmax());
  ldbl dx = (xmax - xmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%t分布関数:自由度=" << df << " 例数=" << n << " 平均値=" << mean << " 標準偏差=" << sd << "\n";
  ldbl x = xmin;
  const ldbl Xeps = dx * 0.01L; // x計算誤差
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    if (fabs(x) < Xeps) x = 0.0L;
    ldbl t = (x - mean) / sd;
    fout << x << " " << (ldbl)n * tdis(t, df) << "\n";
  }
}

void plotclass::plot_tpval(ostream& fout, const option& opt) const
// t分布両側確率出力
{
  constexpr int DFdef = 1;      // 省略時自由度
  constexpr ldbl TMAXdef = 5.0L; // 省略時t最大値
  int df = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : DFdef;
  ldbl tmin = 0.0, tmax = TMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) tmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) tmax = (ldbl)(opt.rmax());
  ldbl dt = (tmax - tmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%t分布両側確率:df=" << df << "\n";
  ldbl t = tmin;
  const ldbl Teps = dt * 0.01L; // t計算誤差
  for (uint i = 0; i <= RDdef; ++i, t += dt)
  {
    if (fabs(t) < Teps) t = 0.0L;
    fout << t << " " << tpval(t, df) << "\n";
  }
}

void plotclass::plot_qdis(ostream& fout, const option& opt) const
// Tukeyのq分布関数出力
{
  constexpr int Gdef = 2;       // 省略時群数
  constexpr int DFdef = 0;      // 省略時自由度
  constexpr ldbl QMINdef = -5.0L; // 省略時q最小値
  constexpr ldbl QMAXdef = 5.0L;  // 省略時q最大値
  int g = (opt.nipara() > 0 && 1 < opt.ipara(0)) ? opt.ipara(0) : Gdef;
  int df = (opt.nipara() > 1 && 0 <= opt.ipara(1)) ? opt.ipara(1) : DFdef;
  ldbl qmin = QMINdef, qmax = QMAXdef;
  if (opt.rmin() != BLANKf) qmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf) qmax = (ldbl)(opt.rmax());
  ldbl dq = (qmax - qmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%Tukeyのq分布関数:g=" << g << " df=" << df << "\n";
  ldbl q = qmin;
  const ldbl Qeps = dq * 0.01L; // q計算誤差
  for (uint i = 0; i <= RDdef; ++i, q += dq)
  {
    if (fabs(q) < Qeps) q = 0.0L;
    fout << q << " " << qdis(q, g, df) << "\n";
  }
}

void plotclass::plot_qpval(ostream& fout, const option& opt) const
// Tukeyのq分布両側確率出力
{
  constexpr int Gdef = 2;       // 省略時群数
  constexpr int DFdef = 0;      // 省略時自由度
  constexpr ldbl QMAXdef = 5.0L; // 省略時q最大値
  int g = (opt.nipara() > 0 && 1 < opt.ipara(0)) ? opt.ipara(0) : Gdef;
  int df = (opt.nipara() > 1 && 0 <= opt.ipara(1)) ? opt.ipara(1) : DFdef;
  ldbl qmin = 0.0, qmax = QMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) qmin = (ldbl)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) qmax = (ldbl)(opt.rmax());
  ldbl dq = (qmax - qmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%Tukeyのq分布両側確率:g=" << g << " df=" << df << "\n";
  ldbl q = qmin;
  const ldbl Qeps = q * 0.01L; // q計算誤差
  for (uint i = 0; i <= RDdef; ++i, q += dq)
  {
    if (fabs(q) < Qeps) q = 0.0L;
    fout << q << " " << qpval(q, g, df) << "\n";
  }
}

void plotclass::plot_w1den(ostream& fout, const option& opt) const
// Wilcoxonの符号付順位分布密度関数出力
{
  constexpr int Ndef = 10;      // 省略時例数
  int n = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : Ndef;
  int tmin = 0, tmaxdef = n * (n + 1) / 2;
  int tmax = tmaxdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin() && opt.rmin() < (double)tmaxdef) tmin = (int)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax() && opt.rmax() <= (double)tmaxdef) tmax = (int)(opt.rmax());
  ldbl p0 = 1.0L / pow(2.0L, n);
  fout << "-s1 -drectangle_solid_fill -x1 -y3 -w80\n";
  fout << "%Wilcoxonの符号付順位分布密度関数:n=" << n << "\n";
  for (int t = tmin; t <= tmax; ++t)
  {
    double x = (double)t;
    fout << x - 0.2 << " 0 " << x + 0.2 << " " << p0 * fw1den(n, t) << "\n";
  }
}

void plotclass::plot_w1dis(ostream& fout, const option& opt) const
// Wilcoxonの符号付順位分布関数出力
{
  constexpr int Ndef = 10;      // 省略時例数
  int n = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : Ndef;
  int tmin = 0, tmaxdef = n * (n + 1) / 2;
  int tmax = tmaxdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin() && opt.rmin() < (double)tmaxdef) tmin = (int)(opt.rmin());
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax() && opt.rmax() <= (double)tmaxdef) tmax = (int)(opt.rmax());
  ldbl p0 = 1.0L / pow(2.0L, n), ppre = 0.0L;
  fout << "-s1 -dline -w80\n";
  fout << "%Wilcoxonの符号付順位分布関数:n=" << n << "\n";
  for (int t = tmin; t <= tmax; ++t)
  {
    ldbl p = p0 * fw1dis(n, t); // 分布確率にする
    fout << t << " " << ppre << "\n";
    fout << t << " " << p << "\n";
    ppre = p;
  }
}

void plotclass::plot_w1pval(ostream& fout, const option& opt) const
// Wilcoxonの符号付順位分布両側確率出力
{
  constexpr int Ndef = 10;      // 省略時例数
  int n = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : Ndef;
  int tmindef = n * (n + 1) / 4, tmaxdef = n * (n + 1) / 2;
  int tmin = tmindef, tmax = tmaxdef;
  if (opt.rmin() != BLANKf && (double)tmindef <= opt.rmin() && opt.rmin() < (double)tmaxdef) tmin = (int)(opt.rmin());
  if (opt.rmax() != BLANKf && (double)tmindef < opt.rmax() && opt.rmax() <= (double)tmaxdef) tmax = (int)(opt.rmax());
  ldbl p0 = 1.0L / pow(2.0L, n);
  fout << "-s1 -dline -w80\n";
  fout << "%Wilcoxonの符号付順位分布両側確率:n" << n << "\n";
  for (int t = tmin; t <= tmax; ++t)
  {
    ldbl p = 2.0L * p0 * fw1dis(n, tmaxdef - t); // 両側確率にする
    if (p > 1.0L) p = 1.0L;
    fout << t << " " << p << "\n";
  }
}

void plotclass::plot_w2den(ostream& fout, const option& opt) const
// Wilcoxonの順位和分布密度関数出力
{
  constexpr int Ndef = 5;       // 省略時例数
  int m = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : Ndef;
  int n = (opt.nipara() > 1 && 0 < opt.ipara(1)) ? opt.ipara(1) : Ndef;
  int tmin = m * (m + 1) / 2, umaxdef = m * n, nmax = n, nmin = m;
  if (m > n) { nmax = m;  nmin = n; }
  int umin = 0, umax = umaxdef;
  if (opt.rmin() != BLANKf && 0 <= opt.rmin() && opt.rmin() < (double)umaxdef) umin = (int)(opt.rmin());
  if (opt.rmax() != BLANKf && 0 < opt.rmax() && opt.rmax() <= (double)umaxdef) umax = (int)(opt.rmax());
  ldbl p0 = 1.0L;
  for (int i = 1; i <= nmin; ++i) p0 *= (ldbl)i / (ldbl)(nmax + i);
  fout << "-s1 -drectangle_solid_fill -x1 -y3 -w80\n";
  fout << "%Wilcoxonの順位和分布密度関数:m=" << m << " n=" << n << "\n";
  for (int u = umin; u <= umax; ++u)
  {
    double t = (double)(u + tmin);
    fout << t - 0.2 << " 0 " << t + 0.2 << " " << p0 * fuden(m, n, u) << "\n";
  }
}

void plotclass::plot_w2dis(ostream& fout, const option& opt) const
// Wilcoxonの順位和分布関数出力
{
  constexpr int Ndef = 5;       // 省略時例数
  int m = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : Ndef;
  int n = (opt.nipara() > 1 && 0 < opt.ipara(1)) ? opt.ipara(1) : Ndef;
  int tmin = m * (m + 1) / 2, umaxdef = m * n, nmax = n, nmin = m;
  if (m > n) { nmax = m;  nmin = n; }
  int umin = 0, umax = umaxdef;
  if (opt.rmin() != BLANKf && 0 <= opt.rmin() && opt.rmin() < (double)umaxdef) umin = (int)(opt.rmin());
  if (opt.rmax() != BLANKf && 0 < opt.rmax() && opt.rmax() <= (double)umaxdef) umax = (int)(opt.rmax());
  ldbl p0 = 1.0L, ppre = 0.0L;
  for (int i = 1; i <= nmin; ++i) p0 *= (ldbl)i / (ldbl)(nmax + i);
  fout << "-s1 -dline -w80\n";
  fout << "%Wilcoxonの順位和分布関数:m=" << m << " n=" << n << "\n";
  for (int u = umin; u <= umax; ++u)
  {
    ldbl p = p0 * fudis(m, n, u); // 分布確率にする
    int t = u + tmin;
    fout << t << " " << ppre << "\n";
    fout << t << " " << p << "\n";
    ppre = p;
  }
}

void plotclass::plot_w2pval(ostream& fout, const option& opt) const
// Wilcoxonの順位和分布両側確率出力
{
  constexpr int Ndef = 5;       // 省略時例数
  int m = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : Ndef;
  int n = (opt.nipara() > 1 && 0 < opt.ipara(1)) ? opt.ipara(1) : Ndef;
  int tmin = m * (m + 1) / 2, umaxdef = m * n, nmax = n, nmin = m;
  if (m > n) { nmax = m;  nmin = n; }
  int umin = 0, umax = umaxdef;
  if (opt.rmin() != BLANKf && 0 <= opt.rmin() && opt.rmin() < (double)umaxdef) umin = (int)(opt.rmin()); // 
  if (opt.rmax() != BLANKf && 0 < opt.rmax() && opt.rmax() <= (double)umaxdef) umax = (int)(opt.rmax());
  ldbl p0 = 1.0L;
  for (int i = 1; i <= nmin; ++i) p0 *= (ldbl)i / (ldbl)(nmax + i);
  fout << "-s1 -dline -w80\n";
  fout << "%Wilcoxonの順位和分布両側確率:m=" << m << " n=" << n << "\n";
  for (int u = umax; u >= umin; --u)
  {
    ldbl p = (u < umax) ? 2.0L * p0 * fudis(m, n, u) : 1.0L; // 両側確率にする
    if (p > 1.0) p = 1.0;
    int t = umax - u + tmin;
    fout << t << " " << p << "\n";
  }
}

void plotclass::plot_wbden(ostream& fout, const option& opt) const
// Weibull分布密度関数出力
{
  constexpr int Adef = 1;       // 省略時a
  constexpr ldbl Ldef = 1.0L;   // 省略時λ
  constexpr ldbl TMAXdef = 5.0; // 省略時t最大値
  int a = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : Adef;
  ldbl lamda = (opt.nfpara() > 0 && 0.0 < opt.fpara(0)) ? opt.fpara(0) : Ldef;
  ldbl tmin = 0.0, tmax = TMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) tmin = opt.rmin();
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) tmax = opt.rmax();
  ldbl dt = (tmax - tmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%Weibull分布密度関数:λ=" << lamda << " a=" << a << "\n";
  double t = tmin;
  for (uint i = 0; i <= RDdef; ++i, t += dt)
  { fout << t << " " << wbden(t, lamda, a) << "\n"; }
}

void plotclass::plot_wbdis(ostream& fout, const option& opt) const
// Weibull分布関数出力
{
  constexpr int Adef = 1;         // 省略時a
  constexpr ldbl Ldef = 1.0L;     // 省略時λ
  constexpr ldbl TMAXdef = 5.0L;  // 省略時t最大値
  int a = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : Adef;
  ldbl lamda = (opt.nfpara() > 0 && 0.0 < opt.fpara(0)) ? opt.fpara(0) : Ldef;
  ldbl tmin = 0.0, tmax = TMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) tmin = opt.rmin();
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) tmax = opt.rmax();
  ldbl dt = (tmax - tmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%Weibull分布関数:λ=" << lamda << " a=" << a << "\n";
  ldbl t = tmin;
  for (uint i = 0; i <= RDdef; ++i, t += dt)
  { fout << t << " " << wbdis(t, lamda, a) << "\n"; }
}

void plotclass::plot_wbhaz(ostream& fout, const option& opt) const
// Weibull分布ハザード関数出力
{
  constexpr int Adef = 1;       // 省略時a
  constexpr ldbl Ldef = 1.0L;   // 省略時λ
  constexpr ldbl TMAXdef = 5.0L; // 省略時t最大値
  int a = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : Adef;
  ldbl lamda = (opt.nfpara() > 0 && 0.0 < opt.fpara(0)) ? opt.fpara(0) : Ldef;
  ldbl tmin = 0.0, tmax = TMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) tmin = opt.rmin();
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) tmax = opt.rmax();
  ldbl dt = (tmax - tmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%Weibull分布ハザード関数:λ=" << lamda << " a=" << a << "\n";
  ldbl t = tmin;
  for (uint i = 0; i <= RDdef; ++i, t += dt)
  { fout << t << " " << wbhaz(t, lamda, a) << "\n"; }
}

void plotclass::plot_wbsur(ostream& fout, const option& opt) const
// Weibull分布生存関数出力
{
  constexpr int Adef = 1;        // 省略時a
  constexpr ldbl Ldef = 1.0L;    // 省略時λ
  constexpr ldbl TMAXdef = 5.0L; // 省略時t最大値
  int a = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : Adef;
  ldbl lamda = (opt.nfpara() > 0 && 0.0 < opt.fpara(0)) ? opt.fpara(0) : Ldef;
  ldbl tmin = 0.0, tmax = TMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) tmin = opt.rmin();
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) tmax = opt.rmax();
  ldbl dt = (tmax - tmin) / (ldbl)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%Weibull分布生存関数:λ=" << lamda << " a=" << a << "\n";
  ldbl t = tmin;
  for (uint i = 0; i <= RDdef; ++i, t += dt)
  { fout << t << " " << wbsur(t, lamda, a) << "\n"; }
}

void plotclass::plot_etc(ostream& fout, const option& opt) const
// その他の関数出力
{
  constexpr int Ndef = 11;      // 省略時n
  constexpr double RMAXdef = 10.0; // 省略時r最大値
  int n = (opt.nipara() > 0 && 0 < opt.ipara(0)) ? opt.ipara(0) : Ndef;
  double rmin = 1.0, rmax = RMAXdef;
  if (opt.rmin() != BLANKf && 0.0 <= opt.rmin()) rmin = opt.rmin();
  if (opt.rmax() != BLANKf && 0.0 < opt.rmax()) rmax = opt.rmax();
  double dr = (rmax - rmin) / (double)RDdef;
  fout << "-s1 -dcurve,spline -w80\n";
  fout << "%分散比−自由度関数:φ=(n-1)(1+γ)2乗/(1+γ2乗) n=" << n << "\n";
  double r = rmin;
  for (uint i = 0; i <= RDdef; ++i, r += dr)
  {
    double phi = (double)(n - 1) * (1.0 + r) * (1.0 + r) / (1 + r * r);
    fout << r << " " << phi << "\n";
  }
}

void plotclass::put(option& opt) // 各種分布を出力する
{
  if (opt.check() != 0) return;
  switch (opt.dtype())
  {
  case disas :
    switch (opt.ftype())
    {
    case funden : plot_asden(cout, opt);
      break;
    case fundis : plot_asdis(cout, opt);
      break;
    default : plot_asden(cout, opt);
    }
    break;
  case disb :
    switch (opt.ftype())
    {
    case funden : plot_bden(cout, opt);
      break;
    case fundis : plot_bdis(cout, opt);
      break;
    case funpval : plot_bpval(cout, opt);
      break;
    default : plot_bden(cout, opt);
    }
    break;
  case disbeta :
    switch (opt.ftype())
    {
    case funden : plot_betaden(cout, opt);
      break;
    case fundis : plot_betadis(cout, opt);
      break;
    default : plot_betaden(cout, opt);
    }
    break;
  case disc :
    switch (opt.ftype())
    {
    case funden : plot_cden(cout, opt);
      break;
    case fundis : plot_cdis(cout, opt);
      break;
    case funhaz : plot_chaz(cout, opt);
      break;
    case funpval : plot_cpval(cout, opt);
      break;
    case funsur : plot_csur(cout, opt);
      break;
    default : plot_cden(cout, opt);
    }
    break;
  case disd : plot_delp(cout, opt);
    break;
  case disdun :
    switch (opt.ftype())
    {
    case fundis : plot_ddis(cout, opt);
      break;
    case funpval : plot_dpval(cout, opt);
      break;
    default : plot_dpval(cout, opt);
    }
    break;
  case dise :
    switch (opt.ftype())
    {
    case funden : plot_eden(cout, opt);
      break;
    case fundis : plot_edis(cout, opt);
      break;
    case funhaz : plot_ehaz(cout, opt);
      break;
    case funsur : plot_esur(cout, opt);
      break;
    default : plot_eden(cout, opt);
    }
    break;
  case disexp : plot_exp(cout, opt);
    break;
  case disetc : plot_etc(cout, opt);
    break;
  case disf :
    switch (opt.ftype())
    {
    case funden : plot_fden(cout, opt);
      break;
    case fundis : plot_fdis(cout, opt);
      break;
    case funpval : plot_fpval(cout, opt);
      break;
    default : plot_fden(cout, opt);
    }
    break;
  case disfis :
    switch (opt.ftype())
    {
    case funden : plot_fisden(cout, opt);
      break;
    case fundis : plot_fisdis(cout, opt);
      break;
    case funpval : plot_fispval(cout, opt);
      break;
    default : plot_fisden(cout, opt);
    }
    break;
  case disln :
    switch (opt.ftype())
    {
    case funden : plot_lnden(cout, opt);
      break;
    case fundis : plot_lndis(cout, opt);
      break;
    case funhaz : plot_lnhaz(cout, opt);
      break;
    case funpval : plot_lnpval(cout, opt);
      break;
    case funsur : plot_lnsur(cout, opt);
      break;
    default : plot_lnden(cout, opt);
    }
    break;
  case dislog :
    switch (opt.ftype())
    {
    case funden : plot_logden(cout, opt);
      break;
    case fundis : plot_logdis(cout, opt);
      break;
    default : plot_logden(cout, opt);
    }
    break;
  case disncc :
    switch (opt.ftype())
    {
    case funden : plot_nccden(cout, opt);
      break;
    case fundis : plot_nccdis(cout, opt);
      break;
    default : plot_nccden(cout, opt);
    }
    break;
  case disncf :
    switch (opt.ftype())
    {
    case funden : plot_ncfden(cout, opt);
      break;
    case fundis : plot_ncfdis(cout, opt);
      break;
    default : plot_ncfden(cout, opt);
    }
    break;
  case disnct :
    switch (opt.ftype())
    {
    case funden : plot_nctden(cout, opt);
      break;
    case fundis : plot_nctdis(cout, opt);
      break;
    default : plot_nctden(cout, opt);
    }
    break;
  case disnpv : plot_npv(cout, opt);
    break;
  case disper : plot_perdis(cout, opt);
    break;
  case displog : plot_plogdis(cout, opt);
    break;
  case dispoi :
    switch (opt.ftype())
    {
    case funden : plot_poiden(cout, opt);
      break;
    case fundis : plot_poidis(cout, opt);
      break;
    default : plot_poiden(cout, opt);
    }
    break;
  case dispoly :plot_poly(cout, opt);
    break;
  case disppv : plot_ppv(cout, opt);
    break;
  case dispro :
    switch (opt.ftype())
    {
    case funden : plot_proden(cout, opt);
      break;
    case fundis : plot_prodis(cout, opt);
      break;
    default : plot_proden(cout, opt);
    }
    break;
  case dissn :
    switch (opt.ftype())
    {
    case funden : plot_snden(cout, opt);
      break;
    case fundis : plot_sndis(cout, opt);
      break;
    case funpinv : plot_snpinv(cout, opt);
      break;
    case funpval : plot_snpval(cout, opt);
      break;
    default : plot_snden(cout, opt);
    }
    break;
  case dist :
    switch (opt.ftype())
    {
    case funden : plot_tden(cout, opt);
      break;
    case fundis : plot_tdis(cout, opt);
      break;
    case funpval : plot_tpval(cout, opt);
      break;
    default : plot_tden(cout, opt);
    }
    break;
  case distuk :
    switch (opt.ftype())
    {
    case fundis : plot_qdis(cout, opt);
      break;
    case funpval : plot_qpval(cout, opt);
      break;
    default : plot_qpval(cout, opt);
    }
    break;
  case disw1 :
    switch (opt.ftype())
    {
    case funden : plot_w1den(cout, opt);
      break;
    case fundis : plot_w1dis(cout, opt);
      break;
    case funpval : plot_w1pval(cout, opt);
      break;
    default : plot_w1den(cout, opt);
    }
    break;
  case disw2 :
    switch (opt.ftype())
    {
    case funden : plot_w2den(cout, opt);
      break;
    case fundis : plot_w2dis(cout, opt);
      break;
    case funpval : plot_w2pval(cout, opt);
      break;
    default : plot_w2den(cout, opt);
    }
    break;
  case diswb :
    switch (opt.ftype())
    {
    case funden : plot_wbden(cout, opt);
      break;
    case fundis : plot_wbdis(cout, opt);
      break;
    case funhaz : plot_wbhaz(cout, opt);
      break;
    case funsur : plot_wbsur(cout, opt);
      break;
    default : plot_wbden(cout, opt);
    }
    break;
  default :
    switch (opt.ftype())
    {
    case funden : plot_nden(cout, opt);
      break;
    case fundis : plot_ndis(cout, opt);
      break;
    case funpinv : plot_npinv(cout, opt);
      break;
    case funpval : plot_npval(cout, opt);
      break;
    default : plot_nden(cout, opt);
    }
  }
}

int main(int argc, char* argv[]) // メインルーチン
{
  plotclass plot;               // メインクラス変数
  option opt(argv[0]);          // オプション
  try
  {
    opt.setdef();               // オプションの初期設定
    for (int i = 1; i < argc; i++)
    {
      if (*(argv[i]) == CHRopt) opt.set(argv[i] + 1);
      else throw runtime_error(string("オプション指定文字が違います→") + string(argv[i]));
    }
    plot.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;
}
