// design.C  Version 7.43  Copyleft (c) 杉本典夫
// Time-stamp: <2025-12-16 08:29:56 tac>
//
// 色々な条件から実験のパラメーターを求める。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法  :g++ -O -s -W design design.C -I/usr/local/include
//               -L/sur/local/lib -ldans -lfun

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

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

const vecstr USAGE =            // プログラム使用法表示用
{
  "design  Version 7.43         Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:design",
  "機能:色々な条件から実験のパラメーターを求めます。",
  "オプション:なし"
};

const string Title = "=== 実験パラメーターの計算 ==="; // 表題
constexpr double Pmin = (double)FDmin; // 最小確率値
constexpr double Pmax = 1.0 - Pmin;    // 最大確率値

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
};

option::option(const char* c_pname) : v_path(), v_pname()
// コンストラクタ
{
  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] == "f")            // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else                          // 間違ったオプション指定
  { err = nonopt(pname(), optstr, USAGE); }
  if (err > 0) throw help_except("ヘルプ表示");
  if (err < 0) throw option_error(string("オプション指定違反→") + optstr);
}

class designclass               // designメインクラス
{
private:
  double get_fdata(const string& message, double low, double up, double def) const;
  // 実数を入力する
  double sqrtn(double x, int n) const; // n次根を返す
  int sel_method(const option& opt, const string& title, int kind, int method) const;
  // 計算メニュー表示と選択
  void n_ttest() const;         // 計量値(t検定)の必要例数計算
  void n_tmcomp() const;        // 計量値(多重比較)の必要例数計算
  void n_ctest(const option& opt) const; // 計数値(χ2乗検定)の必要例数計算
  void n_cmcomp(const option& opt) const; // 計数値(多重比較)の必要例数計算
  void put_tn(int g, double za, double z2b, double sd, double d) const;
  // 計量値(多重比較)の例数出力
  void put_cn(int g, double pi1, double pi2, double za, double z2b) const;
  // 計数値(多重比較)の例数出力
  void n_mcnemar(const option& opt) const; // McNemarの検定の必要例数計算
  void n_cohort() const;
  // 症例対照研究→コホート研究(χ2乗検定)の必要例数計算
  void n_genom(const option& opt) const; // ゲノム疫学(χ2乗検定)の必要例数計算
  void n_w1test() const;   // Wilcoxonの符号付き順位検定の必要例数計算
  void n_utest() const;    // Mann-WhitneyのU検定の必要例数計算
  void n_catest() const;   // Cochran-Armitageの傾向検定の必要例数計算
  void n_cor(const option& opt) const;
  // 相関係数(Pearson、Spearman、四分点相関係数)の必要例数計算
  void n_life(const option& opt) const;
  // 生存時間の比較(単変量生存時間解析)の必要例数計算
  void n_mci() const;  // 計量値の信頼区間の必要例数計算
  void n_rci() const;  // 計数値の信頼区間の必要例数計算
  void n_icc() const;  // 級内相関係数の信頼区間の必要例数計算
  void n_kappa() const;         // 一致係数κの信頼区間の必要例数計算
  void n_hazard() const;        // ハザードλの信頼区間の必要例数計算
  void n_lracl() const;         // 回帰直線の信頼限界の必要例数計算
  void po_ttest() const;        // 計量値(t検定)の検出力計算
  void po_tmcomp() const;       // 計量値(多重比較)の検出力計算
  void po_ctest(const option& opt) const; // 計数値(χ2乗検定)の検出力計算
  void po_cmcomp() const;       // 計数値(多重比較)の検出力計算
  void po_life(const option& opt) const; // 生存時間解析の検出力計算
  void put_tpower(vector<statv<double>>& statvec, double sd, double ta, int df, bool bcput = true) const;
  // 計量値(多重比較)の検出力出力
  void put_cpower(vector<statv<double>>& statvec, double za, bool bcput = true) const;
  // 計数値(多重比較)の検出力出力
  void put_tpowersim(vector<statv<double>>& statvec, double a) const;
  // シミュレーションによる計量値(多重比較)の検出力出力
  void put_cpowersim(vector<statv<double>>& statvec, double a) const;
  // シミュレーションによる計数値(多重比較)の検出力出力
  void diff_ttest() const;      // 計量値(t検定)の検出差計算
  void diff_ctest() const;      // 計数値(χ2乗検定)の検出差計算
  void diff_life() const;       // 生存時間解析の検出差計算
  void absp_ttest() const;      // 計量値の信頼区間の絶対精度計算
  void absp_ctest() const;      // 計数値信頼区間の絶対精度計算
  void absp_lra() const;        // 回帰直線の信頼限界の絶対精度計算
  void aepr_binom() const;      // 有害事象が発現する確率計算
  void aepr_pici() const;       // 母出現率と絶対精度計算
public:
  designclass()                 // コンストラクタ
  {}
  ~designclass()                // デストラクタ
  {}
  int sel(int kind) const;                      // メニュー表示と選択
  void caln(const option& opt, int kind) const; // 実験必要例数の計算
  void calpower(const option& opt, int kind) const; // 検出力の計算
  void caldiff(const option& opt, int kind) const;  // 検出差の計算
  void calabsp(const option& opt, int kind) const;  // 信頼区間の絶対精度の計算
  void calaepr(const option& opt, int kind) const;
  // 有害事象が発現する確率の計算
};

double designclass::get_fdata(const string& message, double low, double up, double def) const
// 実数を入力する
{
  string sline;
  bool ok = false;
  double fdata = 0.0;
  istream* finp = nullptr;
  finp = &cin;
  do
  {
    cout << message << "(" << low << "-" << up << ",[Return]:" << def << "):";
    getsline(sline, *finp);
    strchomp(sline);
    fdata = sline.empty() ? def : strtodf(sline);
    if (fdata < low || up < fdata) beep();
    else ok = true;
  } while (!ok);
  return fdata;
}

double designclass::sqrtn(double x, int n) const // n次根を返す
{
  if (n < 0) throw runtime_error("n次根関数の引数が負だよ!");
  if (n == 0) return 1.0;
  switch (n)
  {
  case 1: return x;
  case 2 : return sqrt(x);
  case 4 : return sqrt(sqrt(x));
  }
  return exp(log(x) / (double)n);
}

int designclass::sel_method(const option& opt, const string& title, int kind, int method) const
// 計算メニュー表示と選択
{
  double fmax = 0.0;
  cout << "--- " << title << " ---\n\n";
  switch (kind)
  {
  case 1 : fmax = 18.0;
    cout << "1:計量値(t検定)\n";
    cout << "2:計量値(多重比較)\n";
    cout << "3:計数値(χ2乗検定)\n";
    cout << "4:計数値(多重比較)\n";
    cout << "5:対応のある2標本計數値(McNemarの検定)\n";
    cout << "6:症例対照研究→コホート研究用(χ2乗検定)\n";
    cout << "7:ゲノム疫学用(χ2乗検定)\n";
    cout << "8:Wilcoxonの符号付き順位検定\n";
    cout << "9:Mann-WhitneyのU検定(勝率=ROC曲線のAUCから)\n";
    cout << "10:Cochran-Armitageの傾向検定\n";
    cout << "11:相関係数(Pearson、Spearman、四分点相関係数)\n";
    cout << "12:生存時間(PHR検定)\n";
    cout << "13:計量値の信頼区間\n";
    cout << "14:計数値の信頼区間\n";
    cout << "15:級内相関係数の信頼区間\n";
    cout << "16:一致係数κの信頼区間\n";
    cout << "17:ハザードλの信頼区間\n";
    cout << "18:回帰直線の信頼限界\n";
    break;
  case 2 : fmax = 5.0;
    cout << "1:計量値(t検定)\n";
    cout << "2:計量値(多重比較)\n";
    cout << "3:計数値(χ2乗検定)\n";
    cout << "4:計数値(多重比較)\n";
    cout << "5:生存時間(PHR検定)\n";
    break;
  case 3 : fmax = 3.0;
    cout << "1:計量値(t検定)\n";
    cout << "2:計数値(χ2乗検定)\n";
    cout << "3:生存時間(PHR検定)\n";
    break;
  case 4 : fmax = 3.0;
    cout << "1:計量値の信頼区間\n";
    cout << "2:計数値の信頼区間\n";
    cout << "3:回帰直線の信頼限界\n";
    break;
  case 5 : fmax = 2.0;
    cout << "1:有害事象が発現する確率\n";
    cout << "2:母出現率と絶対精度\n";
    break;
  default : disperr(opt.pname(), "指定番号不適!");
  }
  cout << "\n0:メインメニューに戻る\n\n";
  return (int)(get_fdata("番号を入力してください", 0.0, fmax, (double)method));
}

void designclass::n_ttest() const // 計量値(t検定)の必要例数計算
{
  cout << "--- 計量値(t検定の場合) ---\n";
  double a = ALdef, po = POdef, d = 1.0, sd = 1.0, r = 1.0;
  while ((a = get_fdata("有意水準(α,0:終了)", 0.0, Pmax, a)) > 0.0)
  {
    po = get_fdata("検出力(1-β)", Pmin, Pmax, po);
    d = get_fdata("検出差(母平均の差δ)", FDmin, FDmax, d);
    sd = get_fdata("標準偏差(σ)", FDmin, FDmax, sd);
    r = get_fdata("2群の場合の例数比(群2/群1)", FDmin, FDmax, r);
    double za = ptoz(a), z2b = ptoz(2.0 * (1.0 - po));
    if (po < 0.5) z2b = -ptoz(2.0 * po);
    double lamda = d / sd;
    double nsqrt = (za + z2b) / lamda;
    double n0 = (1.0 + r) / r * nsqrt * nsqrt;
    cout << "・2群比較(対応無し、群間比較)\n";
    int n1 = (int)(n0 + 0.99);
    int n2 = (int)((double)n1 * r + 0.99);
    double n = n0 + n0 * r;
    double f = (n + 1.0) / (n - 1.0);
    int nt = (int)(n * f + 0.99);
    if (nt < 1) nt = 1;
    double r2 = lamda * lamda / (lamda * lamda + (double)(nt - 2) * (double)nt / (double)n1 / (double)n2);
    double ef = sqrt(r2 / (1.0 - r2));
    cout << "  検出差:δ=" << d << " 標準偏差:σ=" << sd << " 効果量:d=δ/σ=" << lamda << " 効果量f=" << ef << "\n";
    cout << "  相関比:r=" << sqrt(r2) << " 寄与率:r2乗=" << r2 << "\n";
    cout << "  分散既知の場合:群1=" << n1 << " 群2=" << n2 << " (自由度修正済総例数=" << nt << ")\n";
    n1 = (n1 < 2) ? n1 + 1 : (int)(n0 + za * za / 4.0 / (1.0 - 1.0 / (double)n1) + 0.99);
    n2 = (int)(n1 * r + 0.99);
    cout << "  分散未知の場合:群1=" << n1 << " 群2=" << n2 << " (総例数=" << n1 + n2 << ")\n\n";
    cout << "・1群比較(対応有り、群内比較)\n";
    n0 = nsqrt * nsqrt;
    n1 = (int)(n0 + 0.99);
    cout << "  分散既知の場合:例数=" << n1 << "\n";
    n2 = (n1 < 2) ? n1 + 2 : (int)(n0 + za * za / 2.0 / (1.0 - 1.0 / (double)n1) + 0.99);
    cout << "  分散未知の場合:例数=" << n2 << "\n\n";
  }
}

void designclass::n_tmcomp() const // 計量値(多重比較)の必要例数計算
{
  cout << "--- 計量値(多重比較) ---\n";
  int g = 0, gd = 1;
  double a = ALdef, po = POdef, d = 1.0, sd = 1.0;
  while ((g = (int)(get_fdata("群数(0:終了)", 0.0, 100.0, (double)g))) > 1)
  {
    gd = (int)(get_fdata("平均値が異なる群数(0:不明)", 0.0, (double)(g - 1), (double)gd));
    a = get_fdata("有意水準(α)", Pmin, Pmax, a);
    po = get_fdata("検出力(1-β)", Pmin, Pmax, po);
    d = get_fdata("検出差(母平均の差δ)", FDmin, FDmax, d);
    sd = get_fdata("全群共通標準偏差(σ)", FDmin, FDmax, sd);
    cout << "・Fisher型(比較回数=1 有意になる比較回数=1 α=" << a << " 検出力=" << po << ")\n";
    double za = ptoz(a);
    double z2b = (po > 0.5) ? ptoz(2.0 * (1.0 - po)) : -ptoz(2.0 * po);
    put_tn(g, za, z2b, sd, d);
    int nc = g * (g - 1) / 2;   // 比較回数
    int nd = (gd < 1) ? 1 : gd * (2 * g - gd - 1) / 2; // 有意になる比較回数
    double ac = a / (double)nc; // 多重性調整済αエラー
    double poc = 1.0 - sqrtn(1.0 - po, nd); // 多重性調整済検出力
    double zac = ptoz(ac);
    double z2bc = (poc > 0.5) ? ptoz(2.0 * (1.0 - poc)) : -ptoz(2.0 * poc);
    cout << "・Bonferroni型(比較回数=" << nc << " 有意になる比較回数=" << nd << " αc=" << ac << " 検出力=" << poc << ")\n";
    zac = ptoz(ac);
    put_tn(g, zac, z2bc, sd, d);
    cout << "・Scheffe型(比較回数=" << nc << " 有意になる比較回数=" << nd << " α=" << a << " 検出力=" << poc << ")\n";
    zac = sqrt(ptoc(a, g - 1));
    put_tn(g, zac, z2bc, sd, d);
    ac = sqrtn(a, nd);
    poc = sqrtn(po, nc);
    cout << "・ステップ法(比較回数=" << nc << " 有意になる比較回数=" << nd << " αc=" << ac << " 検出力=" << poc << ")\n";
    zac = ptoz(ac);
    z2bc = (poc > 0.5) ? ptoz(2.0 * (1.0 - poc)) : -ptoz(2.0 * poc);
    put_tn(g, zac, z2bc, sd, d);
  }
}

void designclass::n_ctest(const option& opt) const
// 計数値(χ2乗検定)の必要例数計算
{
  cout << "--- 計数値(χ2乗検定) ---\n";
  double a = ALdef, po = POdef, pi1 = 0.6, pi2 = 0.4, s = 1.0;
  while ((a = get_fdata("有意水準(α,0:終了)", 0.0, Pmax, a)) > 0.0)
  {
    po = get_fdata("検出力(1-β)", Pmin, Pmax, po);
    pi1 = get_fdata("群1の出現率(π1)", 0.0, 1.0, pi1);
    while ((pi2 = get_fdata("群2の出現率または理論出現率(π2)", 0.0, 1.0, pi2)) == pi1) disperr(opt.pname(), "π1=π2では計算できません!");
    s = get_fdata("2群の場合の例数比(群2/群1)", FDmin, FDmax, s);
    double pi = (pi1 + s * pi2) / (1.0 + s), dpi = pi1 - pi2;
    double v = pi * (1.0 - pi);
    double za = ptoz(a), z2b = ptoz(2.0 * (1.0 - po));
    if (po < 0.5) z2b = -ptoz(2.0 * po);
    cout << "・2群の出現率の差の検定(対応無し、群間比較)\n";
    cout << "  群1の出現率:π1=" << pi1 << " 群2の出現率:π2=" << pi2;
    cout << " 出現率の差:(π1-π2)=" << dpi << "\n";
    dpi = fabs(dpi);
    double r2 = dpi / (1.0 + s);
    r2 = s * r2 * r2 / v;
    cout << "  連関係数:V=" << sqrt(r2) << " 寄与率:V2乗=" << r2 << "\n";
    double nsqrt = (za * sqrt((1.0 + s) * v) + z2b * sqrt(s * pi1 * (1.0 - pi1) + pi2 * (1.0 - pi2))) / dpi;
    double n0 = nsqrt * nsqrt / s;
    int n1 = (int)(n0 + 0.99);
    int n2 = (int)(n0 * s + 0.99);
    cout << "  連続修正無し:群1=" << n1 << " 群2=" << n2 << " (総例数=" << n1 + n2 << ")\n";
    nsqrt = 1.0 + sqrt(1.0 + (1.0 + s) / s * 2.0 / n0 / dpi);
    n1 = (int)(n0 / 4.0 * nsqrt * nsqrt + 0.99);
    n2 = (int)((double)n1 * s + 0.99);
    cout << "  連続修正有り:群1=" << n1 << " 群2=" << n2 << " (総例数=" << n1 + n2 << ")\n\n";
    double cs = (s * pi1 * (1.0 - pi1) + pi2 * (1.0 - pi2)) / s;
    nsqrt = (za + z2b) / dpi;
    n0 = cs * nsqrt * nsqrt;
    n1 = (int)(n0 + 0.99);
    n2 = (int)(n0 * s + 0.99);
    cout << "  連続修正無し(簡易式):群1=" << n1 << " 群2=" << n2 << " (総例数=" << n1 + n2 << ")\n";
    nsqrt = (za + z2b) / (dpi - sgn(dpi) * 0.5 / n0 * (1.0 + s) / s);
    n0 = cs * nsqrt * nsqrt;
    n1 = (int)(n0 + 0.99);
    n2 = (int)(n0 * s + 0.99);
    cout << "  連続修正有り(簡易式):群1=" << n1 << " 群2=" << n2 << " (総例数=" << n1 + n2 << ")\n\n";
    cout << "・1群の出現率の検定(対応有り、群内比較)\n";
    cout << "  出現率:π1=" << pi1 << " 理論出現率:π2=" << pi2;
    cout << " 出現率の差:(π1-π2)=" << dpi << "\n";
    double v0 = pi2 * (1.0 - pi2), v1 = pi1 * (1.0 - pi1);
    nsqrt = (za * sqrt(v0) + z2b * sqrt(v1)) / dpi;
    n0 = nsqrt * nsqrt;
    cout << "  連続修正無し:例数=" << (int)(n0 + 0.99) << "\n\n";
    cout << "・リスク比の検定(対数リスク差の検定)\n";
    cout << "  群1の出現率:π1=" << pi1 << " 群2の出現率:π2=" << pi2 << "\n";
    double rr = (pi2 > 0.0) ? pi1 / pi2 : 1.0;
    double lrr = log(rr);
    cout << "  リスク比:RR(π1/π2)=" << rr << " 対数リスク差:ln(RR)=" << lrr << "\n";
    cs = (1.0 - pi1) / pi1 + (1.0 - pi2) / pi2 / s;
    nsqrt = (za + z2b) / lrr;
    n0 = cs * nsqrt * nsqrt;
    n1 = (int)(n0 + 0.99);
    n2 = (int)(n0 * s + 0.99);
    cout << "  群1=" << n1 << " 群2=" << n2 << " (総例数=" << n1 + n2 << ")\n\n";
    cout << "・オッズ比の検定(対数オッズ差の検定)\n";
    cout << "  群1の出現率:π1=" << pi1 << " 群2の出現率:π2=" << pi2 << "\n";
    double oz1 = (pi1 > 0.0) ? pi1 / (1.0 - pi1) : 1.0;
    double oz2 = (pi2 > 0.0) ? pi2 / (1.0 - pi2) : 1.0;
    double lor = log(oz1 / oz2);
    cout << "  オッズ比:OR(群1/群2)=" << oz1 / oz2 << " 対数オッズ差:ln(OR)=" << lor << "\n";
    cs = 1.0 / pi1 / (1.0 - pi1) + 1.0 / pi2 / (1.0 - pi2) / s;
    nsqrt = (za + z2b) / lor;
    n0 = cs * nsqrt * nsqrt;
    n1 = (int)(n0 + 0.99);
    n2 = (int)(n0 * s + 0.99);
    cout << "  群1=" << n1 << " 群2=" << n2 << " (総例数=" << n1 + n2 << ")\n\n";
  }
}

void designclass::n_cmcomp(const option& opt) const
// 計数値(多重比較)の必要例数計算
{
  cout << "--- 計数値(多重比較) ---\n";
  int g = 0, gd = 1;
  double a = ALdef, po = POdef, pi1 = 0.6, pi2 = 0.5;
  while ((g = (int)(get_fdata("群数(0:終了)", 0.0, 100.0, (double)g))) > 1)
  {
    gd = (int)(get_fdata("出現率が異なる群数(0:不明)", 0.0, (double)(g - 1), (double)gd));
    a = get_fdata("有意水準(α)", Pmin, Pmax, a);
    po = get_fdata("検出力(1-β)", Pmin, Pmax, po);
    pi1 = get_fdata("群iの出現率(有効率,πi)", 0.0, 1.0, pi1);
    while ((pi2 = get_fdata("群jの出現率(有効率,πj)", 0.0, 1.0, pi2)) == pi1) disperr(opt.pname(), "πi=πjでは計算できません!");
    cout << "・Fisher型(比較回数=1 有意になる比較回数=1 α=" << a << " 検出力=" << po << ")\n";
    double za = ptoz(a);
    double z2b = (po > 0.5) ? ptoz(2.0 * (1.0 - po)) : -ptoz(2.0 * po);
    put_cn(g, pi1, pi2, za, z2b);
    int nc = g * (g - 1) / 2;   // 比較回数
    int nd = (gd < 1) ? 1 : gd * (2 * g - gd - 1) / 2; // 有意になる比較回数
    double ac = a / (double)nc; // 多重性調整済αエラー
    double poc = 1.0 - sqrtn(1.0 - po, nd); // 多重性調整済検出力
     double zac = ptoz(ac);
    double z2bc = (poc > 0.5) ? ptoz(2.0 * (1.0 - poc)) : -ptoz(2.0 * poc);
    cout << "・Bonferroni型(比較回数=" << nc << " 有意になる比較回数=" << nd << " αc=" << ac << " 検出力=" << poc << ")\n";
    zac = ptoz(ac);
    put_cn(g, pi1, pi2, zac, z2bc);
    cout << "・Scheffe型(比較回数=" << nc << " 有意になる比較回数=" << nd << " α=" << a << " 検出力=" << poc << ")\n";
    zac = sqrt(ptoc(a, g - 1));
    put_cn(g, pi1, pi2, zac, z2bc);
    ac = sqrtn(a, nd);
    poc = sqrtn(po, nc);
    cout << "・ステップ法(比較回数=" << nc << " 有意になる比較回数=" << nd << " αc=" << ac << " 検出力=" << poc << ")\n";
    zac = ptoz(ac);
    z2bc = (poc > 0.5) ? ptoz(2.0 * (1.0 - poc)) : -ptoz(2.0 * poc);
    put_cn(g, pi1, pi2, zac, z2bc);
  }
}

void designclass::put_tn(int g, double za, double z2b, double sd, double d) const
// 計量値(多重比較)の例数出力
{
  double nsqrt = (za + z2b) * sd / d;
  double n0 = 2.0 * nsqrt * nsqrt;
  int n = (int)(n0 + 0.99);
  cout << "  分散既知の場合:1群例数=" << n << " (総例数=" << g * n << ")\n";
  n = (n < 2) ? n + 1 : (int)(n0 + za * za / 4.0 / (1.0 - 1.0 / (double)n) + 0.99);
  cout << "  分散未知の場合:1群例数=" << n << " (総例数=" << g * n << ")\n\n";
}

void designclass::put_cn(int g, double pi1, double pi2, double za, double z2b) const
// 計数値(多重比較)の例数出力
{
  double pi = (pi1 + pi2) * 0.5, dpi = fabs(pi1 - pi2);
  double v = pi * (1.0 - pi);
  double nsqrt = (za * sqrt(2.0 * v) + z2b * sqrt(pi1 * (1.0 - pi1) + pi2 * (1.0 - pi2))) / dpi;
  double n0 = nsqrt * nsqrt;
  int n = (int)(n0 + 0.99);
  cout << "  連続修正無し:1群例数=" << n << " (総例数=" << g * n << ")\n";
  nsqrt = 1.0 + sqrt(1.0 + 4.0 / n0 / dpi);
  n = int(n0 / 4.0 * nsqrt * nsqrt + 0.99);
  cout << "  連続修正有り:1群例数=" << n << " (総例数=" << g * n << ")\n\n";
}

void designclass::n_mcnemar(const option& opt) const
// McNemarの検定の必要例数計算
{
  cout << "--- 対応のある2標本計數値(McNemarの検定) ---\n";
  double a = ALdef, po = POdef, pi1 = 0.6, pi2 = 0.5;
  while ((a = get_fdata("有意水準(α,0:終了)", 0.0, Pmax, a)) > 0.0)
  {
    po = get_fdata("検出力(1-β)", Pmin, Pmax, po);
    pi1 = get_fdata("項目1の分類1出現率(π1)", 0.0, 1.0, pi1);
    while ((pi2 = get_fdata("項目2の分類1出現率(π2)", 0.0, 1.0, pi2)) == pi1) disperr(opt.pname(), "π1=π2では計算できません!");
    double pix = pi1 + pi2 - 2.0 * pi1 * pi2;
    double dpi = fabs(pi1 - pi2) / pix;
    double pip = 0.5 + dpi * 0.5, pim = 0.5 - dpi * 0.5;
    double za = ptoz(a), z2b = ptoz(2.0 * (1.0 - po));
    if (po < 0.5) z2b = -ptoz(2.0 * po);
    double n0 = (za + z2b * 2.0 * sqrt(pip * pim)) / (pip - pim);
    n0 = n0 * n0;
    double n = (int)(n0 / pix + 0.99);
    cout << "・対応のある2標本計數値(McNemarの検定)\n";
    cout << "  項目1の分類1出現率:π1=" << pi1 << " 項目2の分類1出現率:π2=" << pi2 << "\n";
    cout << "  π+=" << pip << " π-=" << pim << " 理論出現率=0.5\n";
    cout << "  McNemarの検定の対象例数(連続修正無し)=" << int(n0 + 0.99) << "\n";
    cout << "  全体の例数=" << n << "\n\n";
  }
}

void designclass::n_cohort() const
// 症例対照研究→コホート研究(χ2乗検定)の必要例数計算
{
  cout << "--- 症例対照研究→コホート研究(χ2乗検定) ---\n";
  double a = ALdef, po = POdef, pid = 0.1, sn = 0.7, sp = 0.8, s = 1.0;
  while ((a = get_fdata("有意水準(α,0:終了)", 0.0, Pmax, a)) > 0.0)
  {
    po = get_fdata("検出力(1-β)", Pmin, Pmax, po);
    pid = get_fdata("疾患の発症率(πd)", 0.0, 1.0, pid);
    sn = get_fdata("検査項目の感度(sn)", 0.0, 1.0, sn);
    sp = get_fdata("検査項目の特異度(sp)", 0.0, 1.0, sp);
    s = get_fdata("2群の場合の例数比(群2/群1)", FDmin, FDmax, s);
    double pi1 = pid * sn / (pid * sn + (1.0 - pid) * (1.0 - sp));
    double pi2 = 1.0 - (1.0 - pid) * sp / ((1.0 - pid) * sp + pid * (1.0 - sn));
    double rr = pi1 / pi2;
    double pi = (pi1 + s * pi2) / (1.0 + s), dpi = fabs(pi1 - pi2);
    double v = pi * (1.0 - pi);
    double za = ptoz(a), z2b = ptoz(2.0 * (1.0 - po));
    if (po < 0.5) z2b = -ptoz(2.0 * po);
    double r = sqrt(s) / (1.0 + s) * dpi / sqrt(v);
    cout << "・群1の出現率:π1=" << pi1 << " 群2の出現率:π2=" << pi2 << "\n";
    cout << "  相対危険度=" << rr << " 連関係数:V=" << r << " 寄与率:V2乗=" << r * r << "\n";
    cout << "・2群比較(対応無し、群間比較)\n";
    double nsqrt = (za * sqrt((1.0 + s) * v) + z2b * sqrt(s * pi1 * (1.0 - pi1) + pi2 * (1.0 - pi2))) / dpi;
    double n0 = nsqrt * nsqrt / s;
    int n1 = (int)(n0 + 0.99);
    int n2 = (int)((double)n1 * s + 0.99);
    cout << "  連続修正無し:群1=" << n1 << " 群2=" << n2 << " (総例数=" << n1 + n2 << ")\n";
    nsqrt = 1.0 + sqrt(1.0 + (1.0 + s) / s * 2.0 / n0 / dpi);
    n1 = (int)(n0 / 4.0 * nsqrt * nsqrt + 0.99);
    n2 = (int)((double)n1 * s + 0.99);
    cout << "  連続修正有り:群1=" << n1 << " 群2=" << n2 << " (総例数=" << n1 + n2 << ")\n\n";
    cout << "・1群比較(対応有り、群内比較)\n";
    pi = (pi1 + pi2) / 2.0;
    nsqrt = (za * sqrt(2.0 * pi) + z2b * sqrt(2.0 * pi1 * pi2 / pi)) / dpi;
    n0 = nsqrt * nsqrt;
    cout << "  連続修正有り:例数=" << (int)(n0 + 0.99) << "\n\n";
  }
}

void designclass::n_genom(const option& opt) const
// ゲノム疫学(χ2乗検定)の必要例数計算
{
  cout << "--- ゲノム疫学(χ2乗検定) ---\n";
  double a = ALdef, po = POdef, pid = 0.1, pig = 0.5, rr = 2.0, s = 1.0;
  while ((a = get_fdata("有意水準(α,0:終了)", 0.0, Pmax, a)) > 0.0)
  {
    po = get_fdata("検出力(1-β)", Pmin, Pmax, po);
    pid = get_fdata("疾患の発症率(πd)", 0.0, 1.0, pid);
    pig = get_fdata("遺伝子の頻度(πg)", 0.0, 1.0, pig);
    while ((rr = get_fdata("相対リスク(π1/π2)", FDmin, FDmax, rr)) == 1.0) disperr(opt.pname(), "相対リスク=1では計算できません!");
    s = get_fdata("2群の場合の例数比(群2/群1)", FDmin, FDmax, s);
    double pi2 = pid / (pig * (rr - 1.0) + 1.0);
    double pi1 = rr * pi2;
    double odr = pi1 * (1.0 - pi2) / pi2 / (1.0 - pi1);
    double pi = (pi1 + s * pi2) / (1.0 + s), dpi = fabs(pi1 - pi2);
    double v = pi * (1.0 - pi);
    double za = ptoz(a), z2b = ptoz(2.0 * (1.0 - po));
    if (po < 0.5) z2b = -ptoz(2.0 * po);
    double r = sqrt(s) / (1.0 + s) * dpi / sqrt(v);
    cout << "・群1の出現率:π1=" << pi1 << " 群2の出現率:π2=" << pi2 << "\n";
    cout << "  オッズ比:OR=" << odr << " 連関係数:V=" << r << " 寄与率:V2乗=" << r * r << "\n";
    cout << "・2群比較(対応無し、群間比較)\n";
    double nsqrt = (za * sqrt((1.0 + s) * v) + z2b * sqrt(s * pi1 * (1.0 - pi1) + pi2 * (1.0 - pi2))) / dpi;
    double n0 = nsqrt * nsqrt / s;
    int n1 = (int)(n0 + 0.99);
    int n2 = (int)((double)n1 * s + 0.99);
    cout << "  連続修正無し:群1=" << n1 << " 群2=" << n2 << " (総例数=" << n1 + n2 << ")\n";
    nsqrt = 1.0 + sqrt(1.0 + (1.0 + s) / s * 2.0 / n0 / dpi);
    n1 = (int)(n0 / 4.0 * nsqrt * nsqrt + 0.99);
    n2 = (int)((double)n1 * s + 0.99);
    cout << "  連続修正有り:群1=" << n1 << " 群2=" << n2 << " (総例数=" << n1 + n2 << ")\n\n";
    cout << "・1群比較(対応有り、群内比較)\n";
    pi = (pi1 + pi2) / 2.0;
    nsqrt = (za * sqrt(2.0 * pi) + z2b * sqrt(2.0 * pi1 * pi2 / pi)) / dpi;
    n0 = nsqrt * nsqrt;
    cout << "  連続修正有り:例数=" << (int)(n0 + 0.99) << "\n\n";
  }
}

void designclass::n_w1test() const // Wilcoxnの符号付き順位検定の必要例数計算
{
  cout << "--- Wilcoxnの符号付き順位検定 ---\n";
  double a = ALdef, po = POdef, pi = 0.2, kr = 1.0;
  while ((a = get_fdata("有意水準(α,0:終了)", 0.0, Pmax, a)) > 0.0)
  {
    po = get_fdata("検出力(1-β)", Pmin, Pmax, po);
    pi = get_fdata("符号付き順位平均の割合", FDmin, 1.0, pi);
    kr = get_fdata("同位の補正K", FDmin, FDmax, kr);
    double za = ptoz(a), z2b = ptoz(2.0 * (1.0 - po));
    if (po < 0.5) z2b = -ptoz(2.0 * po);
    double z = za + z2b;
    int n = (int)(kr * z * z / 3.0 / pi / pi + 0.99);
    double ts = pi * (double)n * (double)n;
    cout << "・Wilcoxnの符号付き順位検定\n";
    cout << "  符号付き順位平均の割合=" << pi << "\n";
    cout << "  符号付き順位平均=" << pi * (double)n << "  符号付き順位和=" << int(ts + 0.5) << "\n";
    cout << "  例数(正群の例数＋負群の例数)=" << n << "\n\n";
  }
}

void designclass::n_utest() const // Mann-WhitneyのU検定の必要例数計算(勝率から)
{
  cout << "--- Mann-WhitneyのU検定 ---\n";
  double a = ALdef, po = POdef, pi1 = 0.7, s = 1.0, kr = 1.0;
  while ((a = get_fdata("有意水準(α,0:終了)", 0.0, Pmax, a)) > 0.0)
  {
    po = get_fdata("検出力(1-β)", Pmin, Pmax, po);
    pi1 = get_fdata("群1の勝率(U/(2n))", 0.5 + DBL_EPSILON, 1.0, pi1);
    kr = get_fdata("同位の補正K", FDmin, FDmax, kr);
    s = get_fdata("2群の例数比(群2/群1)", FDmin, FDmax, s);
    double dp = pi1 - 0.5;
    double za = ptoz(a), z2b = ptoz(2.0 * (1.0 - po));
    if (po < 0.5) z2b = -ptoz(2.0 * po);
    double z = za + z2b;
    int n1 = (int)((s + 1.0) * kr * z * z / s / 12.0 / dp / dp + 0.99);
    int n2 = (int)((double)n1 * s + 0.99);
    double u1 = pi1 * (double)n1 * (double)n2;
    double r2 = 3.0 * dp * dp / kr;
    cout << "・Mann-WhitneyのU検定(ROC曲線のAUC)\n";
    cout << "  群1の勝率:U/(2n)=" << pi1 << "  U=" << (int)(u1 + 0.5) << "\n";
    cout << "  順位相関比:r=" << sqrt(r2) << "  寄与率:r2乗=" << r2 << "\n";
    cout << "  群1=" << n1 << " 群2=" << n2 << " (総例数=" << n1 + n2 << ")\n\n";
  }
}

void designclass::n_catest() const // Cochran-Armitageの傾向検定の必要例数計算
{
  cout << "--- Cochran-Armitageの傾向検定 ---\n";
  int g = 0;
  double a = ALdef, po = POdef;
  while ((g = (int)(get_fdata("群数(0:終了)", 0.0, 100.0, (double)g))) > 1.0)
  {
    vecdbl xvec(g, 0.0), pvec(g, 0.0);
    a = get_fdata("有意水準(α)", Pmin, Pmax, a);
    po = get_fdata("検出力(1-β)", Pmin, Pmax, po);
    double p = 0.0, sx = 0.0, sx2 = 0.0;
    for (int i = 0; i < g; i++)
    {
      if (xvec[i] <= 0.0) xvec[i] = i;
      if (pvec[i] <= 0.0) pvec[i] = (double)i * 0.1;
      ostringstream ostr;
      ostr << "群" << i + 1 << "の説明変数x" << i + 1;
      xvec[i] = get_fdata(ostr.str(), -FDmax, FDmax, xvec[i]);
      sx += xvec[i];  sx2 += xvec[i] * xvec[i];
      ostr.str("");
      ostr << "群" << i + 1 << "の反応率p" << i + 1;
      pvec[i] = get_fdata(ostr.str(), 0.0, 1.0, pvec[i]);
      p += pvec[i];
    }
    p /= (double)g;
    double u = 0.0, v = (double)g * sx2 - sx * sx;
    cout << "・Cochran-Armitageの傾向検定\n";
    for (int i = 0; i < g; i++)
    {
      cout << "  群" << i + 1 << "のx" << i + 1 << "=" << xvec[i];
      cout << " p" << i + 1 << "=" << pvec[i] << "\n";
      u += xvec[i] * (pvec[i] - p);
    }
    if (p <= 0.0 || 1.0 <= p || sgn(u) == 0 || sgn(v) == 0)
    { cout << "  総例数計算できず！\n\n"; }
    else
    {
      double r2 = u * u / p / (1.0 - p) / v;
      double za = ptoz(a), z2b = ptoz(2.0 * (1.0 - po));
      if (po < 0.5) z2b = -ptoz(2.0 * po);
      double z = za + z2b;
      int n = (int)(z * z / r2) + 1;
      int ni = (int)((double)n / (double)g + 0.99);
      cout << "  相関比:r=" << sqrt(r2) << "  寄与率:r2乗=" << r2 << "\n";
      cout << "  総例数=" << n << "  1群の例数=" << ni << "\n\n";
    }
  }
}

void designclass::n_cor(const option& opt) const
// 相関係数(Pearson、Spearman、四分点相関係数)の必要例数計算
{
  cout << "--- 相関係数(Pearson、Spearman、四分点相関係数) ---\n";
  double a = ALdef, po = POdef, r = 0.5, s = 1.0;
  while ((a = get_fdata("有意水準(α,0:終了)", 0.0, Pmax, a)) > 0.0)
  {
    po = get_fdata("検出力(1-β)", Pmin, Pmax, po);
    while ((r = get_fdata("相関係数(ρ)", 0.0, 1.0, r)) == 0.0 || r == 1.0) disperr(opt.pname(), "ρ=0,ρ=1では計算できません!");
    s = get_fdata("2群の場合の例数比(群2/群1)", FDmin, FDmax, s);
    cout << "・相関係数:ρ=" << r << " 寄与率:ρ2乗=" << r * r << "\n";
    double pi2 = (2.0 - r * (1.0 + s)) / 4.0 / sqrt(s);
    double pi1 = 1.0 - pi2;
    cout << "  群1の出現率目安:π1=" << pi1 << " 群2の出現率目安:π2=" << pi2 << "\n\n";
    double z = log((1.0 + r) / (1.0 - r));
    double za = ptoz(a), z2b = ptoz(2.0 * (1.0 - po));
    if (po < 0.5) z2b = -ptoz(2.0 * po);
    double z0 = za + z2b;
    int n = (int)(4.0 * z0 / z * z0 / z + 0.99) + 3;
    cout << "・Pearsonの相関係数の時:総例数=" << n << "\n";
    int n1 = (int)((double)n / (1.0 + s) + 0.99);
    int n2 = (int)((double)n1 * s + 0.99);
    cout << "  相関比の時:群1=" << n1 << " 群2=" << n2 << " (総例数=" << n1 + n2 << ")\n";
    n = (int)(z0 / r * z0 / r + 0.99) + 1;
    cout << "  Spearmanの順位相関係数の時:総例数=" << n << "\n";
    n1 = (int)((double)n / (1.0 + s) + 0.99);
    n2 = (int)((double)n1 * s + 0.99);
    cout << "  順位相関比の時:群1=" << n1 << " 群2=" << n2 << " (総例数=" << n1 + n2 << ")\n";
    cout << "  四分点相関係数または連関係数の時:群1=" << n1 << " 群2=" << n2 << " (総例数=" << n1 + n2 << ")\n\n";
  }
}

void designclass::n_life(const option& opt) const
// 生存時間の比較(単変量生存時間解析)の必要例数計算
{
  cout << "--- 生存時間の比較(単変量生存時間解析) ---\n";
  double a = ALdef, po = POdef, lamda1 = 0.2, lamda2 = 0.1, tmax = 0.0, s = 1.0;
  while ((a = get_fdata("有意水準(α,0:終了)", 0.0, Pmax, a)) > 0.0)
  {
    po = get_fdata("検出力(1-β)", Pmin, Pmax, po);
    lamda1 = get_fdata("群1のハザード(λ1)", Pmin, 1.0, lamda1);
    while ((lamda2 = get_fdata("理論ハザードまたは群2のハザード(λ2)", Pmin, 1.0, lamda2)) == lamda1) disperr(opt.pname(), "λ1=λ2では計算できません！");
    tmax = get_fdata("観測期間(0:無限大)", 0.0, FDmax, tmax);
    s = get_fdata("2群の例数比(群1/群2)", FDmin, FDmax, s);
    double hd = lamda1 - lamda2, hr = lamda1 / lamda2;
    double b1 = log(hr);
    double za = ptoz(a), z2b = ptoned(po);
    cout << "・1群のハザード比の検定\n";
    cout << "  群1のハザード:λ1=" << lamda1 << " 理論ハザード:λ2=" << lamda2 << "\n";
    cout << "  ハザード比:HR(λ1/λ2)=" << hr << " 対数ハザード差:ln(HR)=" << b1 << " 観測期間:t=";
    if (tmax > 0.0) cout << tmax << "\n";
    else cout << "∞\n";
    double zt = (za  + z2b) / b1;
    double zt2 = zt * zt;
    double exp1 = (tmax > 0.0) ? 1.0 - exp(-lamda1 * tmax) : 1.0;
    int nd = (int)(zt2 + 0.99), n1 = (int)(zt2 / exp1 + 0.99);
    cout << "  群1=" << n1 << "  (総死亡例数=" << nd << ")\n\n";
    cout << "・パラメトリック・ハザード比検定(群間比較)\n";
    cout << "  群1のハザード:λ1=" << lamda1 << " 群2のハザード:λ2=" << lamda2 << "\n";
    cout << "  ハザード比:HR(λ1/λ2)=" << hr << " 対数ハザード差:ln(HR)=" << b1 << " 観測期間:t=";
    if (tmax > 0.0) cout << tmax << "\n";
    else cout << "∞\n";
    double exp2 = (tmax > 0.0) ? 1.0 - exp(-lamda2 * tmax) : 1.0;
    double n0 = (s * exp1 + exp2) / s / exp1 / exp2 * zt2;
    n1 = (int)(n0 * s + 0.99);
    int n2 = (int)(n0 + 0.99);
    nd = (tmax > 0.0) ? (int)(s * n0 * exp1 + n0 * exp2 + 0.99) : n1 + n2;
    cout << "  群1=" << n1 << " 群2=" << n2 << " (総例数=" << n1 + n2 << "  総死亡例数=" << nd << ")\n\n";
    cout << "・パラメトリック・ハザード差検定(群間比較)\n";
    cout << "  群1のハザード:λ1=" << lamda1 << " 群2のハザード:λ2=" << lamda2 << "\n";
    cout << "  ハザード差:HD(λ1-λ2)=" << hd << " 観測期間:t=";
    if (tmax > 0.0) cout << tmax << "\n";
    else cout << "∞\n";
    double lamdam = (s * lamda1 + lamda2) / (1.0 + s);
    double w1 = lamda1 * lamda1, w2 = lamda2 * lamda2, wm = lamdam * lamdam;
    if (tmax > 0.0)
    {
      w1 /= exp1;  w2 /= exp2;
      wm /= (1.0 - exp(-lamdam * tmax));
    }
    zt = (za * sqrt((1.0 + s) * wm) + z2b * sqrt(s * w1 + w2)) / hd;
    n0 = zt * zt / s;
    n1 = (int)(n0 * s + 0.99);
    n2 = (int)(n0 + 0.99);
    cout << "  群1=" << n1 << " 群2=" << n2 << " (総例数=" << n1 + n2 << ")\n\n";
  }
}

void designclass::n_mci() const // 計量値の信頼区間の必要例数計算
{
  cout << "--- 計量値の信頼区間 ---\n";
  double cc = CCdef, e = 1.0, sd = 1.0, po = 0.0, d = 0.0;
  int nt = 0;
  while ((cc = get_fdata("信頼係数(1-α,0:終了)", 0.0, Pmax, cc)) > 0.0)
  {
    sd = get_fdata("標準偏差(σ)", FDmin, FDmax, sd);
    e = get_fdata("絶対精度(信頼区間幅の半分)", FDmin, FDmax, e);
    po = get_fdata("検出力(1-β,0:検出力考慮せず)", 0.0, Pmax, po);
    if (po > 0.0) d = get_fdata("偏り(0:偏りなし)", 0.0, e - DBL_EPSILON, d);
    nt = get_fdata("母集団の例数(0:無限大)", 0.0, (double)INT_MAX, nt);
    double a = 1.0 - cc;
    double za = ptoz(a), de = e - d;
    double n0 = za * za * sd * sd / de / de;
    if (nt > 0) n0 = (double)nt / ((double)(nt - 1) / n0 + 1.0); // 母集団が無限大ではない時
    else n0 = (int)(n0 + 1.9);
    int n = (int)n0, df = (int)n0 - 1;
    double po0 = 0.0, t0 = ptot(a, df), se = sd / sqrt(n0);
    if (po > 0.0)
    {
      double b = (d == 0.0) ? 1.0 - po : 1.0 - po * de / e;
      double zab = za + ptoz(b);
      n0 = zab * zab * sd * sd / de / de;
      if (nt > 0) n0 = (double)nt / ((double)(nt - 1) / n0 + 1.0);
      else n0 = (int)(n0 + 0.9);
      n = (int)n0;
      df = n - 1;  t0 = ptot(a, df);  se = sd / sqrt(n0);
      double tl = (- e - d) / se + t0, tu = (e - d) / se - t0;
      po0 = ncttolp(tu, df) - ncttolp(tl, df);
      while (po0 < po)
      {
        n++;
        df = n - 1;
        t0 = ptot(a, df);  se = sd / sqrt((double)n);
        tl = (-e - d) / se + t0;  tu = (e - d) / se - t0;
        po0 = ncttolp(tu, df) - ncttolp(tl, df);
      }
    }
    cout << "・信頼係数=" << cc << " 標準偏差=" << sd << " 絶対精度=" << e;
    if (nt > 0) cout << " 母集団の例数=" << nt;
    if (po > 0.0)
    {
      cout << " 偏り=" << d << " 検出力=" << po;
      cout << "\n  例数=" << n << " 実際の検出力=" << po0;
    }
    else
    { cout << "\n  例数=" << n; }
    cout << " 実際の絶対精度=" << t0 * se << "\n\n";
  }
}

void designclass::n_rci() const // 計数値の信頼区間の必要例数計算
{
  cout << "--- 計数値の信頼区間 ---\n";
  double cc = CCdef, p0 = 0.5, pci = 0.25;
  int nt = 0;
  while ((cc = get_fdata("信頼係数(1-α,0:終了)", 0.0, Pmax, cc)) > 0.0)
  {
    p0 = get_fdata("母出現率(π)", 0.0, 1.0, p0);
    pci = get_fdata("絶対精度(信頼区間幅の半分)", DBL_EPSILON, 1.0, pci);
    nt = get_fdata("母集団の例数(0:無限大)", 0.0, (double)INT_MAX, nt);
    double za = ptoz(1.0 - cc);
    double n0 = za * za * p0 * (1.0 - p0) / pci / pci;
    cout << "・信頼係数=" << cc << " 母出現率=" << p0 << " 絶対精度=" << pci;
    if (nt > 0)                 // 母集団が無限大ではない時
    {
      n0 = (double)nt / ((double)(nt - 1) / n0 + 1.0);
      cout << " 母集団の例数=" << nt;
    }
    int nt = (int)(n0 + 0.99);
    int n = (int)((n0 + 0.99) * p0 + 0.5);
    double plow = 0.0, pup = 1.0;
    calbexci(plow, pup, n, nt, cc);
    cout << "\n  例数=" << nt << " 正確な信頼区間:下限=" << plow << " 上限=" << pup << "\n\n";
  }
}

void designclass::n_icc() const // 級内相関係数の信頼区間の必要例数計算
{
  cout << "--- 級内相関係数の信頼区間 ---\n";
  constexpr int Nmax = 10000;
  double cc = CCdef, icc = 0.95, iccl = 0.9, b = 2.0;
  while ((cc = get_fdata("信頼係数(1-α,0:終了)", 0.0, Pmax, cc)) > 0.0)
  {
    b = get_fdata("評価者数(b)", 2.0, (double)INT_MAX, b);
    icc = get_fdata("級内相関係数ICC(3,1)", 0.0 + DBL_EPSILON, 1.0 - DBL_EPSILON, icc);
    iccl = get_fdata("級内相関係数下限ICC(3,1)L", 0.0 + DBL_EPSILON, icc - DBL_EPSILON, iccl);
    double ah = (1.0 - cc) * 0.5;
    double fo = (1.0 + (b - 1.0) * icc) / (1.0 - icc);
    double fl = (1.0 + (b - 1.0) * iccl) / (1.0 - iccl);
    double fol = fo / fl;
    int n = 2, b1 = (int)b - 1;
    int df1 = n - 1, df2 = (n - 1) * ((int)b - 1);
    double f = ptof(ah, df1, df1 * b1);
    while (f > fol && n < Nmax)
    {
      df1 = (++n) - 1;
      f = ptof(ah, df1, df1 * b1);
    }
    fl = fo / f;
    iccl = (fl - 1.0) / (fl + b - 1.0);
    double fu = fo / ptof(1.0 - ah, df1, df1 * b1);
    double iccu = (fu - 1.0) / (fu + b - 1.0);
    cout << "・信頼係数=" << cc << " 評価者数=" << b << " ICC(3,1)=" << icc << "\n";
    cout << "  ICC(3,1)下限=" << iccl << " ICC(3,1)上限=" << iccu << " 例数=" << n << "\n\n";
  }
}

void designclass::n_lracl() const // 回帰直線の信頼限界の必要例数計算
{
  cout << "--- 回帰直線の信頼限界 ---\n";
  double cc = CCdef, e = 1.0, sd = 2.0, r2 = 0.5;
  while ((cc = get_fdata("信頼係数(1-α,0:終了)", 0.0, Pmax, cc)) > 0.0)
  {
    sd = get_fdata("yの標準偏差(SD)", FDmin, FDmax, sd);
    r2 = get_fdata("回帰直線の寄与率(ρ2乗)", 0.0, 1.0, r2);
    e = get_fdata("y推測値の絶対精度(信頼限界幅の半分)", FDmin, FDmax, e);
    double a = 1.0 - cc;
    double za = ptoz(a);
    double vr = sd * sd * (1.0 - r2);
    double n0 = za * za * vr / e / e;
    n0 = (int)(n0 + 2.9);
    int n = (int)n0, df = (int)n0 - 2;
    double t0 = ptot(a, df), se = sqrt(vr / n0);
    cout << "・信頼係数=" << cc << "  yの標準偏差=" << sd << " 残差分散=" << vr << " 絶対精度=" << e;
    cout << "\n  例数:xが平均値の時=" << n << " xが平均値±SDの時=" << 2.0 * n;
    cout << "  実際の絶対精度=" << t0 * se << "\n\n";
  }
}

void designclass::n_kappa() const // 一致係数κの信頼区間の必要例数計算
{
  cout << "--- 一致係数κの信頼区間 ---\n";
  constexpr int Cmax = 1000;
  double cc = CCdef, kappa = 0.8, e = 0.1, c = 2.0;
  while ((cc = get_fdata("信頼係数(1-α,0:終了)", 0.0, Pmax, cc)) > 0.0)
  {
    c = get_fdata("分類数(c)", 2.0, (double)Cmax, c);
    kappa = get_fdata("一致係数κ", 0.0 + DBL_EPSILON, 1.0 - DBL_EPSILON, kappa);
    e = get_fdata("絶対精度(信頼区間幅の半分)", FDmin, 1.0, e);
    double a = 1.0 - cc;
    double za = ptoz(a);
    double pc = 1.0 / c;
    double p0 = kappa * (1.0 - pc) + pc, pc1 = 1.0 - pc;
    double n0 = za * za * p0 * (1.0 - p0) / e / e / pc1 / pc1;
    int n = (int)(n0 + 0.99);
    cout << "・信頼係数=" << cc << " 一致係数:κ=" << kappa << " 絶対精度=" << e;
    cout << " 例数=" << n << "\n\n";
  }
}

void designclass::n_hazard() const // ハザードλの信頼区間の必要例数計算
{
  cout << "--- ハザードλの信頼区間 ---\n";
  double cc = CCdef, lam = 0.1, e = 0.05, t = 0.0;
  while ((cc = get_fdata("信頼係数(1-α,0:終了)", 0.0, Pmax, cc)) > 0.0)
  {
    lam = get_fdata("母ハザード(λ)", 0.0, FDmax, lam);
    e = get_fdata("絶対精度(信頼区間幅の半分)", DBL_EPSILON, lam, e);
    t = get_fdata("観測期間(0:∞)", 0.0, FDmax, t);
    double a = 1.0 - cc;
    double za = ptoz(a);
    double ex = (t > 0.0) ? 1.0 - exp(-lam * t) : 1.0;
    double n0 = za * za * lam * lam / e / e / ex;
    cout << "・信頼係数=" << cc << " 母ハザード:λ=" << lam << " 絶対精度=" << e;
    cout << " 観測期間:t=";
    if (t > 0.0) cout << t;
    else cout << "∞";
    cout << "\n  平均生存時間:1/λ=" << 1.0 / lam << " 下限=" << 1.0 / (lam + e) << " 上限=" << 1.0 / (lam - e);
    int n = (int)(n0 + 0.99);
    cout << " 例数=" << n << "\n\n";
  }
}

void designclass::po_ttest() const // 計量値(t検定)の検出力計算
{
  cout << "--- 計量値(t検定の場合) ---\n";
  double a = ALdef, n1 = 2.0, m1 = 1.0, sd1 = 1.0, n2 = 2.0, m2 = 0.0, sd2 = 1.0;
  while ((a = get_fdata("有意水準(α,0:終了)", 0.0, Pmax, a)) > 0.0)
  {
    n1 = get_fdata("群1の例数", 2.0, (double)INT_MAX, n1);
    m1 = get_fdata("群1の平均値", -FDmax, FDmax, m1);
    sd1 = get_fdata("群1の標準偏差", FDmin, FDmax, sd1);
    n2 = get_fdata("群2の例数(0:1群比較,群内比較)", 0.0, (double)INT_MAX, n2);
    int df = (int)n1 - 1;
    double dm = fabs(m1), sd = sd1, ta = 0.0, t2b = 0.0;
    double d = dm / sd;         // 効果量
    double lamda = d * sqrt(n1); // 非心度
    int nt = (int)(n1 + n2);
    if (n2 >= 1.0)              // 2群比較
    {
      m2 = get_fdata("群2の平均値", -FDmax, FDmax, m2);
      sd2 = get_fdata("群2の標準偏差", FDmin, FDmax, sd2);
      cout << "・2群比較(対応無し、群間比較)\n";
      cout << "  群1:例数=" << n1 << " 標準偏差=" << sd1 << "\n";
      cout << "  群2:例数=" << n2 << " 標準偏差=" << sd2 << "\n";
      df = nt - 2;
      dm = fabs(m1 - m2);
      sd = sqrt((sd1 * sd1 * (n1 - 1.0) + sd2 * sd2 * (n2 - 1.0)) / (double)df);
      d = dm / sd;
      lamda = d / sqrt(1.0 / n1 + 1.0 / n2);
      ta = ptot(a, df);
      t2b = lamda - ta;
      double r2 = d * d / (d * d + (double)(nt - 2) * (double)nt / n1 / n2);
      cout << "  相関比:r=" << sqrt(r2) << " 寄与率:r2乗=" << r2 << "\n";
    }
    else
    {
      cout << "・1群比較(対応有り、群内比較)\n";
      cout << "  例数=" << n1 << " 標準偏差=" << sd1 << "\n";
      ta = ptot(a, df);
      t2b = lamda - ta;
    }
    double b = (t2b > 0.0) ? ttop(t2b, df) * 0.5 : 1.0 - ttop(-t2b, df) * 0.5;
    double b0 = ncttolp(ta, df, lamda);
    cout << "  検出差:d=" << dm << " 共通標準偏差:SD=" << sd << " 効果量:d/SD=" << d << " 非心度:λ=" << lamda << "\n";
    cout << "  検出力(正規近似)=" << 1.0 - b << "(β=" << b << ")\n";
    cout << "  検出力(非心t分布)=" << 1.0 - b0 << "(β=" << b0 << ")\n\n";
  }
}

void designclass::po_tmcomp() const // 計量値(多重比較)の検出力計算
{
  cout << "--- 計量値(多重比較) ---\n";
  int g = 0, gold = 0;
  double a = ALdef;
  statv<double> dummys;
  dummys.set(10, 0.0, 1.0);
  vector<statv<double>> statvec;
  while ((g = (int)(get_fdata("群数(0:終了)", 0.0, 100.0, (double)g))) > 1)
  {
    a = get_fdata("有意水準(α)", Pmin, Pmax, a);
    if (statvec.size() > 0 && gold != g) statvec.clear();
    if (statvec.empty()) statvec.resize(g, dummys);
    int nt = 0;
    double sst = 0.0;
    for (int i = 0; i < g; ++i)
    {
      ostringstream ostr;
      ostr << "群" << i + 1 << "の例数";
      int n = (int)(get_fdata(ostr.str(), 2.0, (double)INT_MAX, (double)(statvec[i].n())));
      ostr.str("");
      ostr << "群" << i + 1 << "の平均値";
      double mean = get_fdata(ostr.str(), -FDmax, FDmax, statvec[i].mean());
      ostr.str("");
      ostr << "群" << i + 1 << "の標準偏差";
      double sd = get_fdata(ostr.str(), FDmin, FDmax, statvec[i].sd());
      statvec[i].set(n, mean, sd);
      nt += n;
      sst += statvec[i].ss();
    }
    int df = nt - g;
    sst = sqrt(sst / (double)df);
    cout << "・Fisher型(α=" << a << ")\n";
    put_tpower(statvec, sst, ptot(a, df), df, false);
    int nc = g * (g - 1) / 2;
    double ac = a / (double)nc; // 多重性調整済αエラー
    double tac = ptot(ac, df);
    cout << "・Bonferroni型(αc=" << ac << ")\n";
    put_tpower(statvec, sst, tac, df);
    cout << "・Scheffe型(α=" << a << ")\n";
    int df1 = g - 1;
    tac = sqrt(ptof(a, df1, df) * df1);
    put_tpower(statvec, sst, tac, df);
    put_tpowersim(statvec, a);
    gold = g;
  }
}

void designclass::po_ctest(const option& opt) const
// 計数値(χ2乗検定)の検出力計算
{
  cout << "--- 計数値(χ2乗検定) ---\n";
  double n1 = 2.0, n2 = 2.0, a = ALdef, pi1 = 0.7, pi2 = 0.5;
  while ((a = get_fdata("有意水準(α,0:終了)", 0.0, Pmax, a)) > 0.0)
  {
    n1 = get_fdata("群1の例数", 2.0, (double)INT_MAX, n1);
    pi1 = get_fdata("群1の出現率(有効率,π1)", 0.0, 1.0, pi1);
    n2 = get_fdata("群2の例数(0:1群比較,群内比較)", 0.0, (double)INT_MAX, n2);
    if (n2 >= 1.0)              // 2群比較
    {
      while ((pi2 = get_fdata("群2の出現率(有効率,π2)", 0.0, 1.0, pi2)) == pi1) disperr(opt.pname(), "群1の出現率=群2の出現率では計算できません!");
      double dpi = fabs(pi1 - pi2);
      double s = n2 / n1;
      double pi = (pi1 + s * pi2) / (s + 1.0);
      double v = pi * (1.0 - pi);
      double za = ptoz(a);
      double tmp1 = za * sqrt((s + 1.0) * v);
      double tmp2 = sqrt(s * pi1 * (1.0 - pi1) + pi2 * (1.0 - pi2));
      double z2b = (sqrt(s * dpi * dpi * n1) - tmp1) / tmp2;
      double po = nedp(z2b);
      double r = sqrt(s) / (1.0 + s) * dpi / sqrt(v);
      cout << "・2群比較(対応無し、群間比較)\n";
      cout << "  群1:例数=" << n1 << " 出現率=" << pi1 << "\n";
      cout << "  群2:例数=" << n2 << " 出現率=" << pi2 << "\n";
      cout << "  連関係数:V=" << r << " 寄与率:V2乗=" << r * r << "\n";
      cout << "  連続修正無し:検出力=" << po << "(β=" << 1.0 - po << ")\n";
      double n1t = n1 * 0.9;
      double tmp3 = 1.0 + sqrt(1.0 + 2.0 * (s + 1.0) / s / n1t / dpi);
      n1t = 4.0 * n1 / tmp3 / tmp3;
      z2b = (sqrt(s * dpi * dpi * n1t) - tmp1) / tmp2;
      po = nedp(z2b);
      cout << "  連続修正有り:検出力=" << po << "(β=" << 1.0 - po << ")\n\n";
    }
    else                        // 1群比較
    {
      while ((pi2 = get_fdata("理論出現率(有効率,π0)", 0.0, 1.0, pi2)) == pi1) disperr(opt.pname(), "群1の出現率=理論出現率では計算できません!");
      double dpi = fabs(pi1 - pi2);
      double za = ptoz(a);
      double z2b = (dpi * sqrt(n1) - za * sqrt(pi2 * (1.0 - pi2))) / sqrt(pi1 * (1.0 - pi1));
      double po = nedp(z2b);
      cout << "・1群比較(群内比較)\n";
      cout << "  例数=" << n1 << " 出現率=" << pi1 << " 理論出現率=" << pi2 << "\n";
      cout << "  検出力=" << po << "(β=" << 1.0 - po << ")\n\n";
    }
  }
}

void designclass::po_cmcomp() const // 計数値(多重比較)の検出力計算
{
  cout << "--- 計数値(多重比較) ---\n";
  int g = 0, gold = 0;
  double a = ALdef;
  statv<double> dummys(10, 5.0, 0.0);
  vector<statv<double>> statvec;
  while ((g = (int)(get_fdata("群数(0:終了)", 0.0, 100.0, (double)g))) > 1)
  {
    a = get_fdata("有意水準(α)", Pmin, Pmax, a);
    if (statvec.size() > 0 && gold != g) statvec.clear();
    if (statvec.empty()) statvec.resize(g, dummys);
    int nt = 0;
    for (int i = 0; i < g; ++i)
    {
      ostringstream ostr;
      ostr << "群" << i + 1 << "の例数";
      int n = (int)(get_fdata(ostr.str(), 2.0, (double)INT_MAX, (double)(statvec[i].n())));
      ostr.str("");
      ostr << "群" << i + 1 << "の出現率(π" << i + 1 << ")";
      double mean = get_fdata(ostr.str(), 0.0, 1.0, statvec[i].mean());
      statvec[i].set(n, mean, 0.0);
      nt += n;
    }
    cout << "・Fisher型(α=" << a << ")\n";
    put_cpower(statvec, ptoz(a), false);
    int nc = g * (g - 1) / 2;
    double ac = a / (double)nc; // 多重性調整済αエラー
    double zac = ptoz(ac);
    cout << "・Bonferroni型(αc=" << ac << ")\n";
    put_cpower(statvec, zac);
    cout << "・Scheffe型(α=" << a << ")\n";
    zac = sqrt(ptoc(a, g - 1));
    put_cpower(statvec, zac);
    put_cpowersim(statvec, a);
    gold = g;
  }
}

void designclass::po_life(const option& opt) const // 生存時間解析の検出力計算
{
  cout << "--- 生存時間解析(PHR検定) ---\n";
  double n1 = 2.0, n2 = 2.0, nd1 = 1.0, nd2 = 1.0;
  double a = ALdef, lamda1 = 0.1, lamda2 = 0.2;
  while ((a = get_fdata("有意水準(α,0:終了)", 0.0, Pmax, a)) > 0.0)
  {
    n1 = get_fdata("群1の例数", 2.0, (double)INT_MAX, n1);
    nd1 = get_fdata("群1の死亡例数", 1.0, (double)INT_MAX, nd1);
    lamda1 = get_fdata("群1のハザード(λ1)", FDmin, 1.0, lamda1);
    n2 = get_fdata("群2の例数", 2.0, (double)INT_MAX, n2);
    nd2 = get_fdata("群2の死亡例数", 1.0, (double)INT_MAX, nd2);
    while ((lamda2 = get_fdata("群2のハザード(λ2)", Pmin, 1.0, lamda2)) == lamda1) disperr(opt.pname(), "λ1=λ2では計算できません！");
    {
      double hr = lamda2 / lamda1, za = ptoz(a);
      double b1 = log(hr), seb1 = sqrt(1.0 / nd1 + 1.0 / nd2);
      double z2b = fabs(b1) / seb1 - za;
      double po = nedp(z2b);
      cout << "・PHR検定(対応無し、2群間比較)\n";
      cout << "  群1:例数=" << n1 << " 死亡例数=" << nd1 << " ハザード:λ1=" << lamda1 << "\n";
      cout << "  群2:例数=" << n2 << " 死亡例数=" << nd2 << " ハザード:λ2=" << lamda2 << "\n";
      cout << "  ハザード比:HR(λ2/λ1)=" << hr << " 対数ハザード差:ln(HR)=" << b1 << "\n";
      cout << "  検出力=" << po << "(β=" << 1.0 - po << ")\n\n";
    }
  }
}

void designclass::put_tpower(vector<statv<double>>& statvec, double sd, double ta, int df, bool bcput) const
// 計量値(多重比較)の検出力出力
{
  int g = statvec.size();
  double bc = 1.0, bc0 = 1.0;
  for (int i = 0; i < g - 1; ++i)
  {
    for (int j = i + 1; j < g; ++j)
    {
      double d = fabs(statvec[i].mean() - statvec[j].mean());
      double lamda = d / sd / sqrt(1.0 / (double)(statvec[i].n()) + 1.0 / (double)(statvec[j].n()));
      double t2b = lamda - ta;
      double b = (t2b > 0.0) ? ttop(t2b, df) * 0.5 : 1.0 - ttop(-t2b, df) * 0.5;
      double b0 = ncttolp(ta, df, lamda);
      cout << "  群" << i + 1 << "-群" << j + 1 << ":" << "検出差=" << d << " 共通標準偏差=" << sd << " 非心度=" << lamda << "\n";
      cout << "  検出力(正規近似)=" << 1.0 - b << "(β=" << b << ")\n";
      cout << "  検出力(非心t分布)=" << 1.0 - b0 << "(β=" << b0 << ")\n";
      bc *= b;  bc0 *= b0;
    }
  }
  if (bcput && g > 2)
  {
    cout << "  全体の検出力(正規近似)=" << 1.0 - bc << "(β=" << bc << ")\n";
    cout << "  全体の検出力(非心t分布)=" << 1.0 - bc0 << "(β=" << bc0 << ")\n";
  }
  cout << "\n";
}

void designclass::put_cpower(vector<statv<double>>& statvec, double za, bool bcput) const
// 計数値(多重比較)の検出力出力
{
  int g = statvec.size();
  double bc = 1.0, bcc = 1.0;
  for (int i = 0; i < g - 1; ++i)
  {
    for (int j = i + 1; j < g; ++j)
    {
      double pi1 = statvec[i].mean(), pi2 = statvec[j].mean();
      if (sgndiff(pi1, pi2) == 0)
      {
        cout << "  群" << i + 1 << "-群" << j + 1 << ":2群の出現率が等しいので検出力を計算できない！\n";
      }
      else
      {
        double dpi = fabs(pi1 - pi2);
        double n1 = statvec[i].n();
        double r = (double)(statvec[j].n()) / n1;
        double pi = (pi1 + r * pi2) / (r + 1.0);
        double v = pi * (1.0 - pi);
        double tmp1 = za * sqrt((r + 1.0) * v);
        double tmp2 = sqrt(r * pi1 * (1.0 - pi1) + pi2 * (1.0 - pi2));
        double z2b = (sqrt(r * dpi * dpi * n1) - tmp1) / tmp2;
        double po = nedp(z2b);
        double b = 1.0 - po;
        cout << "  群" << i + 1 << "-群" << j + 1 << ":検出力=" << po << "(β=" << b << ")";
        bc *= b;
        double n1t = n1 * 0.9;
        double tmp3 = 1.0 + sqrt(1.0 + 2.0 * (r + 1.0) / r / n1t / dpi);
        n1t = 4.0 * n1 / tmp3 / tmp3;
        z2b = (sqrt(r * dpi * dpi * n1t) - tmp1) / tmp2;
        po = nedp(z2b);
        b = 1.0 - po;
        cout << " 修正有検出力=" << po << "(β=" << b << ")\n";
        bcc *= b;
      }
    }
  }
  if (bcput && g > 2)
  {
    cout << "  全体の検出力=" << 1.0 - bc << "(β=" << bc << ") 修正有検出力=" << 1.0 - bcc << "(β=" << bcc << ")\n";
  }
  cout << endl;
}

void designclass::put_tpowersim(vector<statv<double>>& statvec, double a) const
// シミュレーションによる計量値(多重比較)の検出力出力
{
  int g = statvec.size();
  vector<statv<double>> stativec(g + 1); // 1回1回の基礎統計量
  int nasig = 0, nfsig = 0, nbsig = 0, nssig = 0; // 有意になった回数
  int nc = g * (g - 1) / 2;     // 多重比較の比較回数
  srand(time((long *)0));       // 乱数のタネをまく
  constexpr int Rdef = 10000;   // 繰り返し回数
  for (int i = 0; i < Rdef; ++i)
  {
    for (int j = 0; j < g + 1; ++j) stativec[j].clear();
    for (int j = 0; j < g; ++j)
    {
      int nj = statvec[j].n();
      double mj = statvec[j].mean(), sdj = statvec[j].sd();
      for (int k = 0; k < nj; ++k) stativec[j].sigma(randn(mj, sdj));
    }
    double sa = 0.0;
    for (int j = 0; j < g; ++j)
    {
      stativec[g] += stativec[j];
      int nj = stativec[j].n();
      double mj = stativec[j].mean();
      sa += nj * mj * mj;
    }
    double st = stativec[g].ss();
    int nt = stativec[g].n();
    double mt = stativec[g].mean();
    sa -= nt * mt * mt;
    double sr = st - sa;
    int dft = nt - 1, dfa = g - 1, dfr = nt - g;
    const double Seps = DBL_EPSILON * fmax(FDmin, st);
    if (sa < Seps || sr < Seps || dfr < 1) continue; // 検定できず
    double vr = sr / dfr;
    double f = (sa / dfa) / vr;
    double pr = ftop(f, dfa, dfr);
    if (pr < a) ++nasig;
    bool fsig = false, bsig = false, ssig = false;
    for (int j = 0; (bsig == false || ssig == false) && j < g - 1; ++j)
    {
      int nj = stativec[j].n();
      double mj = stativec[j].mean();
      for (int k = j + 1; (bsig == false || ssig == false) && k < g; ++k)
      {
        int nk = stativec[k].n();
        double mk = stativec[k].mean();
        double dm = mj - mk;
        f = dm * dm / (1.0 / (double)nj + 1.0 / (double)nk) / vr;
        pr = ftop(f, 1, dfr);
        if (j == 0 && k == 1)   // Fisher型は1度だけ検定
        {
          if (pr < a) fsig = true;
        }
        pr *= nc;
        if (pr < a) bsig = true;
        f = dm * dm / dfa / (1.0 / (double)nj + 1.0 / (double)nk) / vr;
        pr = ftop(f, dfa, dfr);
        if (pr < a) ssig = true;
      }
    }
    if (fsig) ++nfsig;
    if (bsig) ++nbsig;
    if (ssig) ++nssig;
  }
  double po = (double)nasig / (double)Rdef;
  cout << "・シミュレーションによる検出力(繰り返し回数:" << Rdef << ")\n";
  cout << "  一元配置分散分析の検出力 =" << po << "(β=" << 1.0 - po << ")\n";
  po = (double)nfsig / (double)Rdef;
  cout << "  Fisher型の検出力(群1-群2)=" << po << "(β=" << 1.0 - po << ")\n";
  po = (double)nbsig / (double)Rdef;
  cout << "  Bonferroni型の全体検出力 =" << po << "(β=" << 1.0 - po << ")\n";
  po = (double)nssig / (double)Rdef;
  cout << "  Scheffe型の全体検出力    =" << po << "(β=" << 1.0 - po << ")\n\n";
}

void designclass::put_cpowersim(vector<statv<double>>& statvec, double a) const
// シミュレーションによる計数値(多重比較)の検出力出力
{
  int g = statvec.size();
  ftable ftabi(g, 2);           // 1回1回の度数分布表
  int nasig = 0, nfsig = 0, nbsig = 0, nssig = 0; // 有意になった回数
  int nc = g * (g - 1) / 2;     // 多重比較の比較回数
  srand(time((long *)0));       // 乱数のタネをまく
  constexpr int Rdef = 10000;   // 繰り返し回数
  for (int i = 0; i < Rdef; ++i)
  {
    ftabi.fclear();
    for (int j = 0; j < g; ++j)
    {
      int nj = statvec[j].n();
      double pij = statvec[j].mean();
      double mj = (double)nj * pij, sdj = sqrt((double)nj * pij * (1.0 - pij));
      int np = (int)(randn(mj, sdj)); // 正規乱数を利用して度数を設定する
      ftabi.assign(j + 1, 1, np);  ftabi.assign(j + 1, 2, nj - np);
    }
    ftabi.set();
    double csqr = 0.0;
    int df = 0;
    ctest(csqr, df, ftabi, false);
    if (sgn(csqr) <= 0) continue; // 検定できず
    double pr = ctop(csqr, df);
    if (pr < a) ++nasig;
    int nt = ftabi.ftotal(), dff = ftabi.nc() - 1;
    bool fsig = false, bsig = false, ssig = false;
    for (int j = 1; (bsig == false || ssig == false) && j < g; ++j)
    {
      double xj = (double)(ftabi.frtotal(j));
      for (int k = j + 1; (bsig == false || ssig == false) && k <= g; ++k)
      {
        double xk = (double)(ftabi.frtotal(k));
        csqr = 0.0;
        for (int l = 1; l <= (int)(ftabi.nc()); ++l)
        {
          double pd = (double)(ftabi.f(j, l)) / xj - (double)(ftabi.f(k, l)) / xk;
          csqr += pd * pd / (double)(ftabi.fctotal(l));
        }
        csqr *= nt * xj * xk / (xj + xk);
        pr = ctop(csqr, dff);
        if (j == 1 && k == 2)   // Fisher型は1度だけ検定
        {
          if (pr < a) fsig = true;
        }
        pr *= nc;
        if (pr < a) bsig = true;
        pr = ctop(csqr, df);
        if (pr < a) ssig = true;
      }
    }
    if (fsig) ++nfsig;
    if (bsig) ++nbsig;
    if (ssig) ++nssig;
  }
  double po = (double)nasig / (double)Rdef;
  cout << "・シミュレーションによる検出力(繰り返し回数:" << Rdef << ")\n";
  cout << "  χ2乗検定(" << g << "×2)の検出力   =" << po << "(β=" << 1.0 - po << ")\n";
  po = (double)nfsig / (double)Rdef;
  cout << "  Fisher型の検出力(群1:群2)=" << po << "(β=" << 1.0 - po << ")\n";
  po = (double)nbsig / (double)Rdef;
  cout << "  Bonferroni型の全体検出力 =" << po << "(β=" << 1.0 - po << ")\n";
  po = (double)nssig / (double)Rdef;
  cout << "  Scheffe型の全体検出力    =" << po << "(β=" << 1.0 - po << ")\n\n";
}

void designclass::diff_ttest() const // 計量値(t検定)の検出差計算
{
  cout << "--- 計量値(t検定の場合) ---\n";
  double a = ALdef, po = POdef, n1 = 2.0, sd1 = 1.0, n2 = 2.0, sd2 = 1.0;
  while ((a = get_fdata("有意水準(α,0:終了)", 0.0, Pmax, a)) > 0.0)
  {
    po = get_fdata("検出力(1-β)", Pmin, Pmax, po);
    n1 = get_fdata("群1の例数", 2.0, (double)INT_MAX, n1);
    sd1 = get_fdata("群1の標準偏差", FDmin, FDmax, sd1);
    int df = (int)n1 - 1;
    double sd = sd1, ta = 0.0, t2b = 0.0, lamda = 0.0;
    n2 = get_fdata("群2の例数(0:1群比較,群内比較)", 0.0, (double)INT_MAX, n2);
    int nt = (int)(n1 + n2);
    if (n2 >= 1.0)              // 2群比較
    {
      sd2 = get_fdata("群2の標準偏差", FDmin, FDmax, sd2);
      cout << "・2群比較(対応無し、群間比較)\n";
      cout << "  群1:例数=" << n1 << " 標準偏差=" << sd1 << "\n";
      cout << "  群2:例数=" << n2 << " 標準偏差=" << sd2 << "\n";
      df = nt - 2;
      sd = sqrt((sd1 * sd1 * (n1 - 1.0) + sd2 * sd2 * (n2 - 1.0)) / (double)df);
      double s = n2 / n1;
      ta = ptot(a, df);  t2b = ptot(2.0 * (1.0 - po), df);
      double nc = ta * ta / 4.0 / (1.0 - 1.0 / (double)nt);
      if (n1 <= nc) nc = 0.0;
      lamda = (ta + t2b) * sqrt((1.0 + s) / s / (n1 - nc));
      double r2 = lamda * lamda / (lamda * lamda + (double)(nt - 2) * (double)nt / n1 / n2);
      cout << "  相関比:r=" << sqrt(r2) << " 寄与率:r2乗=" << r2 << "\n";
    }
    else
    {
      cout << "・1群比較(対応有り、群内比較)\n";
      cout << "  例数=" << n1 << " 標準偏差=" << sd1 << "\n";
      ta = ptot(a, df), t2b = ptot(2.0 * (1.0 - po), df);
      double nc = (n1 > 1.0) ? ta * ta / 2.0 / (1.0 - 1.0 / n1) : 0.0;
      if (n1 <= nc) nc = 0.0;
      lamda = (ta + t2b) / sqrt(n1 - nc);
    }
    double d = lamda * sd;
    cout << "  検出差:d=" << d << " 共通標準偏差:SD=" << sd << " 効果量:d/SD=" << lamda << "\n\n";
  }
}

void designclass::diff_ctest() const // 計数値(χ2乗検定)の検出差計算
{
  cout << "--- 計数値(χ2乗検定) ---\n";
  double n1 = 2.0, n2 = 2.0, a = ALdef, po = POdef, pi = 0.5;
  while ((a = get_fdata("有意水準(α,0:終了)", 0.0, Pmax, a)) > 0.0)
  {
    po = get_fdata("検出力(1-β)", Pmin, Pmax, po);
    n1 = get_fdata("群1の例数", 2.0, (double)INT_MAX, n1);
    n2 = get_fdata("群2の例数(0:1群比較,群内比較)", 0.0, (double)INT_MAX, n2);
    if (n2 >= 1.0)              // 2群比較
    {
      pi = get_fdata("全体の出現率(有効率,π0)", 0.0, 1.0, pi);
      double s = n2 / n1;
      double v = pi * (1.0 - pi);
      double za = ptoz(a), z2b = ptoz(2.0 * (1.0 - po));
      double dpi = (za + z2b) * sqrt((1.0 + s) / s * v / n1);
      double pi1 = pi - dpi * s / (1.0 + s);
      double pi2 = pi1 + dpi;
      double r = sqrt(s) / (1.0 + s) * dpi / sqrt(v);
      cout << "・2群比較(対応無し、群間比較)\n";
      cout << "  連続修正無し:検出差=" << dpi << "\n";
      cout << "  連関係数:V=" << r << " 寄与率:V2乗=" << r * r << "\n";
      cout << "  全体の出現率=" << pi << "\n";
      cout << "  群1:例数=" << n1 << " 出現率=" << pi1 << "\n";
      cout << "  群2:例数=" << n2 << " 出現率=" << pi2 << "\n";
      cout << "  群2の出現率/群1の出現率=" << pi2 / pi1 << "\n\n";
    }
    else                        // 1群比較
    {
      double za = ptoz(a), z2b = ptoz(2.0 * (1.0 - po));
      double dpi = (za + z2b) / sqrt(n1);
      double pi1 = 0.5 + dpi / 2.0, pi2 = 0.5 - dpi / 2.0;
      cout << "・1群比較(群内比較)\n";
      cout << "  連続修正無し:検出差=" << dpi << "\n";
      cout << "  有の理論出現率=0.5\n";
      cout << "  例数=" << n1 << " 有の出現率=" << pi1 << " 無の出現率=" << pi2 << "\n";
      cout << "  有の出現率/無の出現率=" << pi1 / pi2 << "\n\n";
    }
  }
}

void designclass::diff_life() const // 生存時間解析の検出差計算
{
  cout << "--- 生存時間解析(PHR検定) ---\n";
  double nd1 = 1.0, nd2 = 1.0, a = ALdef, po = POdef;
  while ((a = get_fdata("有意水準(α,0:終了)", 0.0, Pmax, a)) > 0.0)
  {
    po = get_fdata("検出力(1-β)", Pmin, Pmax, po);
    nd1 = get_fdata("群1の死亡例数", 1.0, (double)INT_MAX, nd1);
    nd2 = get_fdata("群2の死亡例数", 1.0, (double)INT_MAX, nd2);
    double beta = 1.0 - po;
    double za = ptoz(a), z2b = ptoz(beta * 2.0);
    double seb1 = sqrt(1.0 / nd1 + 1.0 / nd2);
    double b1 = (za + z2b) * seb1;
    cout << "・PHR検定(対応無し、2群間比較)\n";
    cout << "  群1:死亡例数=" << nd1 << "  群2:死亡例数=" << nd2 << "\n";
    cout << "  検出ハザード比:HR(λ2/λ1)=" << exp(b1) << "  対数ハザード差:ln(HR)=" << b1 << "\n\n";
  }
}

void designclass::absp_ttest() const // 計量値の信頼区間の絶対精度計算
{
  cout << "--- 計量値の信頼区間 ---\n";
  double cc = CCdef, sd = 1.0;
  int n = 20;
  while ((cc = get_fdata("信頼係数(1-α,0:終了)", 0.0, Pmax, cc)) > 0.0)
  {
    sd = get_fdata("標準偏差(σ)", FDmin, FDmax, sd);
    n = get_fdata("例数", 2.0, (double)INT_MAX, (double)n);
    double al = 1.0 - cc;
    int df = n - 1;
    double t0 = ptot(al, df);
    double e = t0 * sd / sqrt((double)n);
    cout << "・信頼係数=" << cc << " 標準偏差=" << sd << " 例数=" << n;
    cout << " 絶対精度=" << e << "\n\n";
  }
}

void designclass::absp_ctest() const // 計数値の信頼区間の絶対精度計算
{
  cout << "--- 計数値の信頼区間 ---\n";
  double cc = CCdef, p0 = 0.5;
  int nt = 20;
  while ((cc = get_fdata("信頼係数(1-α,0:終了)", 0.0, Pmax, cc)) > 0.0)
  {
    p0 = get_fdata("母出現率(π)", DBL_EPSILON, 1.0, p0);
    nt = get_fdata("例数", 2.0, (double)INT_MAX, (double)nt);
    double plow = 0.0, pup = 1.0;
    double pci = calbnmci(plow, pup, p0, nt, cc);
    cout << "・信頼係数=" << cc << " 母出現率=" << p0 << " 例数=" << nt << " 絶対精度=" << pci << "\n";
    int n = (int)((double)nt * p0 + 0.5);
    pci = calbexci(plow, pup, n, nt, cc);
    cout << "  正確な信頼区間:下限=" << plow << " 上限=" << pup << " 信頼区間幅=" << pci << "\n\n";
  }
}

void designclass::absp_lra() const // 回帰直線の信頼限界の絶対精度計算
{
  cout << "--- 回帰直線の信頼限界 ---\n";
  double cc = CCdef, vr = 1.0;
  int n = 20;
  while ((cc = get_fdata("信頼係数(1-α,0:終了)", 0.0, Pmax, cc)) > 0.0)
  {
    vr = get_fdata("残差分散(VR)", FDmin, FDmax, vr);
    n = get_fdata("例数", 2.0, (double)INT_MAX, (double)n);
    double al = 1.0 - cc;
    int df = n - 2;
    double t0 = ptot(al, df);
    double e = t0 * sqrt(vr / (double)n);
    cout << "・信頼係数=" << cc << " 残差分散=" << vr << " 例数=" << n;
    cout << "\n  絶対精度:xが平均値の時=" << e << "  xが平均値±SDxの時≒" << e * sqrt(2.0) << "\n\n";
  }
}

void designclass::aepr_binom() const // 有害事象が発現する確率計算
{
  cout << "--- 有害事象が発現する確率 ---\n";
  double p0 = 0.1;
  int nt = 20;
  while ((p0 = get_fdata("母出現率(π,0:終了)", 0.0, 1.0, p0)) > 0.0)
  {
    nt = get_fdata("例数", 2.0, (double)INT_MAX, (double)nt);
    double pr = binom(0, nt, p0);
    cout << "・母出現率=" << p0 << " 例数=" << nt << "\n";
    cout << "  有害事象0例の時の片側有意確率=" << pr << "\n";
    cout << "  有害事象が1例以上発現する確率=" << 1.0 - pr << "\n\n";
  }
}

void designclass::aepr_pici() const // 母出現率と絶対精度計算
{
  cout << "--- 母出現率と絶対精度 ---\n";
  double cc = CCdef, rci = 1.0;
  int nt = 20;
  while ((cc = get_fdata("信頼係数(1-α,0:終了)", 0.0, Pmax, cc)) > 0.0)
  {
    rci = get_fdata("絶対精度／母出現率比", FDmin, FDmax, rci);
    nt = get_fdata("例数", 2.0, (double)INT_MAX, (double)nt);
    double za = ptoz(1.0 - cc);
    double p0 = za * za / (rci * rci * (double)nt + za * za);
    if (p0 >= 1.0) p0 = 1.0 - DBL_EPSILON;
    double plow = 0.0, pup = 1.0;
    double pci = calbnmci(plow, pup, p0, nt, cc, false); // 連続補正なし
    cout << "・信頼係数=" << cc << " 絶対精度／母出現率比=" << rci << " 例数=" << nt << "\n";
    cout << "  母出現率=" << p0 << " 絶対精度=" << pci << " 信頼区間:下限=" << plow << " 上限=" << pup << "\n";
    int n = (int)((double)nt * p0 + 0.5);
    pci = calbexci(plow, pup, n, nt, cc);
    cout << "  正確な信頼区間:下限=" << plow << " 上限=" << pup << " 信頼区間幅=" << pci << "\n\n";
  }
}

int designclass::sel(int kind) const // メニュー表示と選択
{
  puttitle(cout, Title);
  cout << "\n";
  cout << "1:必要例数の計算\n";
  cout << "2:検出力の計算(βエラーの計算)\n";
  cout << "3:検出差の計算\n";
  cout << "4:信頼区間の絶対精度の計算\n";
  cout << "5:有害事象の確率計算\n\n";
  cout << "9:マニュアル表示\n";
  cout << "0:コマンド終了\n\n";
  return (int)(get_fdata("番号を入力してください", 0.0, 9.0, (double)kind));
}

void designclass::caln(const option& opt, int kind) const // 実験必要例数の計算
{
  int method = 0;
  while ((method = sel_method(opt, "実験必要例数の計算", kind, method)) != 0)
  {
    switch (method)
    {
      case 1 : n_ttest();
        break;
      case 2 : n_tmcomp();
        break;
      case 3 : n_ctest(opt);
        break;
      case 4 : n_cmcomp(opt);
        break;
      case 5 : n_mcnemar(opt);
        break;
      case 6 : n_cohort();
        break;
      case 7 : n_genom(opt);
        break;
      case 8 : n_w1test();
        break;
      case 9 : n_utest();
        break;
      case 10 : n_catest();
        break;
      case 11 : n_cor(opt);
        break;
      case 12 : n_life(opt);
        break;
      case 13 : n_mci();
        break;
      case 14 : n_rci();
        break;
      case 15 : n_icc();
        break;
      case 16 : n_kappa();
        break;
      case 17 : n_hazard();
        break;
      case 18 : n_lracl();
        break;
      default : disperr(opt.pname(), "指定番号不適！");
    }
  }
}

void designclass::calpower(const option& opt, int kind) const // 検出力の計算
{
  int method = 0;
  while ((method = sel_method(opt, "検出力の計算", kind, method)) != 0)
  {
    switch (method)
    {
      case 1 : po_ttest();
        break;
      case 2 : po_tmcomp();
        break;
      case 3 : po_ctest(opt);
        break;
      case 4 : po_cmcomp();
        break;
      case 5 : po_life(opt);
        break;
      default : disperr(opt.pname(), "指定番号不適！");
    }
  }
}

void designclass::caldiff(const option& opt, int kind) const // 検出差の計算
{
  int method = 0;
  while ((method = sel_method(opt, "検出差の計算", kind, method)) != 0)
  {
    switch (method)
    {
      case 1 : diff_ttest();
        break;
      case 2 : diff_ctest();
        break;
      case 3 : diff_life();
        break;
      default : disperr(opt.pname(), "指定番号不適！");
    }
  }
}

void designclass::calabsp(const option& opt, int kind) const // 信頼区間の絶対精度の計算
{
  int method = 0;
  while ((method = sel_method(opt, "信頼区間の絶対精度の計算", kind, method)) != 0)
  {
    switch (method)
    {
      case 1 : absp_ttest();
        break;
      case 2 : absp_ctest();
        break;
      case 3 : absp_lra();
        break;
      default : disperr(opt.pname(), "指定番号不適！");
    }
  }
}

void designclass::calaepr(const option& opt, int kind) const // 有害事象の確率計算
{
  int method = 0;
  while ((method = sel_method(opt, "有害事象の確率計算", kind, method)) != 0)
  {
    switch (method)
    {
      case 1 : aepr_binom();
        break;
      case 2 : aepr_pici();
        break;
      default : disperr(opt.pname(), "指定番号不適！");
    }
  }
}

int main(int argc, char* argv[]) // メインルーチン
{
  designclass design;           // メインクラス変数
  option opt(argv[0]);          // オプション
  try
  {
    opt.setdef();               // オプションの初期設定
    int kind = 0;
    while ((kind = design.sel(kind)) != 0)
    {
      switch (kind)
      {
      case 1 : design.caln(opt, kind);
        break;
      case 2 : design.calpower(opt, kind);
        break;
      case 3 : design.caldiff(opt, kind);
        break;
      case 4 : design.calabsp(opt, kind);
        break;
      case 5 : design.calaepr(opt, kind);
        break;
      case 9 : help(opt.pname(), USAGE);
        break;
      default : disperr(opt.pname(), "指定番号不適！");
      }
    }
    disperr(opt.pname(), "終了！");
  }
  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;
}
