// bnmsim.C  Version 7.42  Copyleft (c) 杉本典夫
// Time-stamp: <2025-11-01 11:17:14 tac>
//
// 二項検定のシミュレーション
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法  :g++ -O -s -W bnmsim bnmsim.C -I/usr/local/include
//               -L/usr/local/lib -ldans -lfun

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

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

const vecstr USAGE =            // プログラム使用法
{
  "bnmsim  Version 7.42         Copyleft (Vo\\oV) 2025 杉本典夫",
  "書式:bnmsim  [<オプション> …]  [<データファイル> …]",
  "機能:二項検定のシミュレーション。",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定する。(0≦信頼係数＜1)",
  "  -f=<ファイル名>:指定されたファイルからオプションを読み込む。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示する。(-h:-h=man)",
  "  -n[=<例数>]:例数を指定する。(0＜例数,オプション省略時:-n=10)",
  "  -p[=<母出現率>]:母出現率を指定する。(0＜母出現率＜1,オプション省略時:-p=0.5)",
  "  -two[={-|+}]:両側検定の指定。(-two:two=+)",
  "    -:片側検定  +:両側検定(オプション省略時)"
};

const string Title = "=== 二項検定のシミュレーション ==="; // 表題
constexpr int Ndef = 10;        // 省略時例数
constexpr double Pdef = 0.5;    // 省略時母出現率

double cal_pci(double& pl, double& pu, double p0, int nt, double za, bool cor)
// 正規近似法による出現率の信頼区間とその幅を返す
{
  if (nt < 1 || za <= 0.0) return 0.0;
  double pci = za * sqrt(p0 * (1.0 - p0) / (double)nt);
  if (cor && sgn(pci) != 0) pci += 0.5 / (double)nt;
  pl = p0 - pci;  pu = p0 + pci;
  return pci;                   // 信頼幅を返す
}

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  int v_n;                      // 例数
  double v_p;                   // 母出現率
  bool v_two;                   // 両側検定フラッグ
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  double cc() const             // 信頼係数を返す
  { return v_cc; }
  int n() const                 // 例数を返す
  { return v_n; }
  double p() const              // 母出現率を返す
  { return v_p; }
  bool two() const              // 両側検定フラッグを返す
  { return v_two; }
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
};

option::option(const char* c_pname) : v_path(), v_pname(), v_cc(CCdef), v_n(Ndef), v_p(Pdef), v_two(true)
// コンストラクタ
{
  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] == "cc")           // 信頼係数の指定
  {
    v_cc = opt[1].empty() ? CCdef : strtodf(opt[1]);
    if (v_cc < 0.0 || 1.0 <= v_cc) err = badopt(pname(), optstr, USAGE, 0.0, 1.0);
  }
  else if (opt[0] == "f")       // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else if (opt[0] == "n")       // 例数の指定
  {
    char* ptr = nullptr;
    v_n = opt[1].empty() ? Ndef : strtol(opt[1].c_str(), &ptr, 0);
    if (v_n < 1) err = badopt(pname(), optstr, USAGE, 1.0, (double)INT_MAX);
  }
  else if (opt[0] == "p")       // 母出現率の指定
  {
    v_p = opt[1].empty() ? Pdef : strtodf(opt[1]);
    if (v_p <= 0.0 || 1.0 <= v_p) err = badopt(pname(), optstr, USAGE, 0.0, 1.0);
  }
  else if (opt[0] == "two")     // 両側検定の指定
  {
    if (opt[1].empty() || opt[1] == "+") v_two = true;
    else if (opt[1] == "-") v_two = false;
    else err = badopt(pname(), optstr, USAGE);
  }
  else                          // 間違ったオプション指定
  { err = nonopt(pname(), optstr, USAGE); }
  if (err > 0) throw help_except("ヘルプ表示");
  if (err < 0) throw option_error(string("オプション指定違反→") + optstr);
}

class varclass                  // 変数クラス
{
private:
  vecdbl bnmvec;                // 二項分布ベクトル
  vecdbl cpvec;                 // Clopper-Pearson法のクロス集計表
  vecdbl ndpvec; // 正規近似法(母出現率使用・連続修正有)のクロス集計表
  vecdbl ndcvec; // 正規近似法(連続修正有)のクロス集計表
  vecdbl ndvec;  // 正規近似法(連続修正無)のクロス集計表
  double v_mpl;                 // 出現率下限
  double v_mpu;                 // 出現率上限
  double v_pl;                  // 下側累積確率
  double v_pu;                  // 上側累積確率
public:
  varclass(int n, double p0, double cc0, bool two); // コンストラクタ
  ~varclass()                   // デストラクタ
  {}
  void set(int n, double p0, double cc0, bool two); // 変数をセットする
  const vecdbl& cp() const    // Clopper-Pearson法のクロス集計表を返す
  { return cpvec; }
  const vecdbl& ndp() const
  // 正規近似法(母出現率使用・連続修正有)のクロス集計表を返す
  { return ndpvec; }
  const vecdbl& ndc() const // 正規近似法(連続修正有)のクロス集計表を返す
  { return ndcvec; }
  const vecdbl& nd() const // 正規近似法(連続修正無)のクロス集計表を返す
  { return ndvec; }
  double mpl() const            // 出現率下限を返す
  { return v_mpl; }
  double mpu() const            // 出現率上限を返す
  { return v_mpu; }
  double pl() const             // 下側累積確率を返す
  { return v_pl; }
  double pu() const             // 上側累積確率を返す
  { return v_pu; }
};

varclass::varclass(int n, double p0, double cc0, bool two) : bnmvec(n + 1, 0.0), cpvec(4, 0.0), ndpvec(4, 0.0), ndcvec(4, 0.0), ndvec(4, 0.0), v_mpl(0.0), v_mpu(0.0), v_pl(0.0), v_pu(0.0)
// コンストラクタ
{
  double alow = two ? (1.0 - cc0) * 0.5 : 0.0;
  double aup = cc0 + alow;
  double cipr = 0.0;
  int flow = 0, fup = n;
  bool lowsch = two ? true : false, upsch = true;
  v_pl = 0.0;
  for (int i = 0; i <= n; ++i)
  {
    double fpr = (double)(bnmden(i, n, (ldbl)p0));
    bnmvec[i] = fpr;
    if (lowsch)                 // 信頼限界下限探索
    {
      if (v_pl + fpr <= alow) v_pl += fpr;
      else { flow = i;  lowsch = false; }
    }
    if (upsch)                  // 信頼限界上限探索
    {
      cipr += fpr;
      if (cipr >= aup) { fup = i;  upsch = false; }
    }
  }
  v_mpl = (double)flow / (double)n; // 出現率下限
  v_mpu = (double)fup / (double)n;  // 出現率上限
  v_pu = (fup < n) ? 1.0 - cipr : 0.0; // 下側累積確率
}

void varclass::set(int n, double p0, double cc0, bool two) // 変数をセットする
{
  size_t nocip = 0, nocis = 0;
  double alpha = 1.0 - cc0;
  if (two == false) alpha *= 2.0;
  double cc = 1.0 - alpha, za = ptoz(alpha);
  for (int i = 0; i <= n; ++i)
  {
    double m = (double)i / (double)n;
    if (two) nocip = (mpl() <= m && m <= mpu()) ? 0 : 1;
    else nocip = (m <= mpu()) ? 0 : 1; // 母集団の信頼区間チェック
    double ml = 0.0, mu = 1.0;
    double pci = calbexci(ml, mu, i, n, cc); // Clopper-Pearson法
    if (two) nocis = (ml <= p0 && p0 <= mu) ? 0 : 1;
    else nocis = (ml <= p0) ? 0 : 1;
    size_t idx = nocip * 2 + nocis;
    cpvec[idx] += bnmvec[i];
    pci = cal_pci(ml, mu, p0, n, za, true); // 母出現率を使用した正規近似法
    ml = m - pci;  mu = m + pci;
    if (two) nocis = (ml <= p0 && p0 <= mu) ? 0 : 1;
    else nocis = (ml <= p0) ? 0 : 1;
    idx = nocip * 2 + nocis;
    ndpvec[idx] += bnmvec[i];
    pci = cal_pci(ml, mu, m, n, za, true); // 正規近似法(連続修正有)
    if (two) nocis = (ml <= p0 && p0 <= mu) ? 0 : 1;
    else nocis = (ml <= p0) ? 0 : 1;
    idx = nocip * 2 + nocis;
    ndcvec[idx] += bnmvec[i];
    pci = cal_pci(ml, mu, m, n, za, false); // 正規近似法(連続修正無)
    if (two) nocis = (ml <= p0 && p0 <= mu) ? 0 : 1;
    else nocis = (ml <= p0) ? 0 : 1;
    idx = nocip * 2 + nocis;
    ndvec[idx] += bnmvec[i];
  }
}

class bnmsimclass               // bnmsimメインクラス
{
private:
  void put_bnm(ostream& fout, const option& opt, const varclass& var) const;
  // 検定結果と推定結果のクロス集計表を出力する
  void put_table(ostream& fout, const option& opt, const vecdbl& tmpvec) const;
  // 分割表と診断率を出力する
public:
  bnmsimclass()                 // コンストラクタ
  {}
  ~bnmsimclass()                // デストラクタ
  {}
  void put(const option& opt);  // シミュレーション結果を出力する
};

void bnmsimclass::put_bnm(ostream& fout, const option& opt, const varclass& var) const
// 検定結果と推定結果のクロス集計表を出力する
{
  puttitle(fout, Title);
  fout << "\n○検定:二項検定  母出現率=" << opt.p() << "  例数=" << opt.n();
  if (opt.two()) fout << "  両側検定\n";
  else fout << "  片側検定\n";
  double pro = opt.cc() * 100.0;
  fout << "  " << pro << "%信頼区間:出現率=" << var.mpl() << "〜" << var.mpu() << "\n";
  fout << "  " << 100.0 - pro << "%棄却域:下側累積確率=" << var.pl() << "  上側累積確率=" << var.pu() << "\n";
  double alpha = 1.0 - opt.cc();
  if (opt.two() == false) alpha *= 2.0;
  double cc = 1.0 - alpha, za = ptoz(alpha);
  double ml = 0.0, mu = 1.0;
  int fp = (int)((opt.n()) * opt.p() + 0.111);
  double pci = calbexci(ml, mu, fp, opt.n(), cc);
  if (opt.two() == false) mu = 1.0;
  fout << "\n・推定:Clopper-Pearson法\n";
  fout << "  出現率=" << opt.p() << "の時の" << pro << "%信頼区間:出現率=" << ml << "〜" << mu << "\n";
  put_table(fout, opt, var.cp());
  pci = cal_pci(ml, mu, opt.p(), opt.n(), za, true);
  if (opt.two() == false) mu = 1.0;
  fout << "\n・推定:正規近似法(母出現率使用)\n";
  fout << "  出現率=" << opt.p() << "の時の" << pro << "%信頼区間:出現率=" << ml << "〜" << mu << "\n";
  put_table(fout, opt, var.ndp());
  pci = cal_pci(ml, mu, opt.p(), opt.n(), za, true);
  if (opt.two() == false) mu = 1.0;
  fout << "\n・推定:正規近似法(連続修正有)\n";
  fout << "  出現率=" << opt.p() << "の時の" << pro << "%信頼区間:出現率=" << ml << "〜" << mu << "\n";
  put_table(fout, opt, var.ndc());
  pci = cal_pci(ml, mu, opt.p(), opt.n(), za, false);
  if (opt.two() == false) mu = 1.0;
  fout << "\n・推定:正規近似法(連続修正無)\n";
  fout << "  出現率=" << opt.p() << "の時の" << pro << "%信頼区間:出現率=" << ml << "〜" << mu << "\n";
  put_table(fout, opt, var.nd());
}

void bnmsimclass::put_table(ostream& fout, const option& opt, const vecdbl& tmpvec) const
// 分割表と診断率を出力する
{
  double r1total = tmpvec[0] + tmpvec[1], r2total = tmpvec[2] + tmpvec[3];
  double c1total = tmpvec[0] + tmpvec[2], c2total = tmpvec[1] + tmpvec[3];
  double pro = opt.cc() * 100.0;
  fout << "検定結果＼推定結果    信頼区間      棄却域        合計\n";
  int wlen = 54, w = 12;
  putline(fout, wlen);
  fout << "信頼区間          " << setw(w) << tmpvec[0] << setw(w) << tmpvec[1] << setw(w) << r1total << "\n";
  fout << "棄却域            " << setw(w) << tmpvec[2] << setw(w) << tmpvec[3] << setw(w) << r2total << "\n";
  putline(fout, wlen);
  fout << "合計              " << setw(w) << c1total << setw(w) << c2total << setw(w) << 1 << "\n";
  pro = (r1total > 0.0) ? tmpvec[0] / r1total : 0.0;
  fout << "感度=" << tmpvec[0] << "/" << r1total << "=" << pro;
  pro = (r2total > 0.0) ? tmpvec[3] / r2total : 0.0;
  fout << " 特異度=" << tmpvec[3] << "/" << r2total << "=" << pro << "\n";
  pro = tmpvec[0] + tmpvec[3];
  fout << "正診率=" << tmpvec[0] << "+" << tmpvec[3] << "=" << pro << "\n";
}

void bnmsimclass::put(const option& opt) // シミュレーション結果を出力する
{
  varclass var(opt.n(), opt.p(), opt.cc(), opt.two()); // 変数オブジェクト
  var.set(opt.n(), opt.p(), opt.cc(), opt.two()); // 変数をセットする
  put_bnm(cout, opt, var);
  cout << "\n" << endl;
}

int main(int argc, char* argv[]) // メインルーチン
{
  bnmsimclass bnmsim;           // メインクラス変数
  option opt(argv[0]);          // オプション
  try                           // 例外監視
  {
    opt.setdef();               // オプションの初期設定
    bool fok = true;            // ファイル処理モードフラッグ
    for (int i = 1; i < argc; ++i)
    {
      if (*(argv[i]) == CHRopt) opt.set(argv[i] + 1);
    }
    bnmsim.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;
}
