// powsim.C  Version 7.43  Copyleft (c) 杉本典夫
// Time-stamp: <2025-11-03 14:40:30 tac>
//
// シミュレーションによって検出力を計算する。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法  :g++ -O -s -W powsim powsim.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 =            // プログラム使用法
{
  "powsim  Version 7.43         Copyleft (Vo\\oV) 2025 杉本典夫",
  "書式:powsim  [<オプション> …]  [<データファイル> …]",
  "機能:シミュレーションによって検出力を計算する。",
  "<オプション>",
  "  -a[=<有意水準>]:有意水準αを指定する。",
  "    (0<有意水準α<1,オプション省略時:-a=0.05)",
  "  -am[={-|+}]:ANOVAが有意の時だけ多重比較実施。",
  "    (-am:-am=+,オプション省略時:-am=-)",
  "  -f=<ファイル名>:指定されたファイルからオプションを読み込む。",
  "  -g[=<群数>]:群の数を指定する。(1<群数,オプション省略時:-g=3)",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示する。(-h:-h=man)",
  "  -m=<平均値1>[,<平均値2>,…]:群の平均値を指定する。",
  "    指定しない群の平均値は0になる。",
  "  -n=<例数1>[,<例数2>,…]:群の例数を指定する。(1<例数)",
  "    指定しない群の例数は100になる。",
  "  -r[=<繰り返し数>]:繰り返し数を指定する。",
  "    (0<繰り返し数,オプション省略時:-r=10000)",
  "  -sd=<標準偏差1>,[<標準偏差2>,…]群の標準偏差を指定する。(0<標準偏差)",
  "    指定しない群の標準偏差は1になる。",
  "  -seed[=<整数>]乱数のタネを指定する。(0<整数,オプション省略時:-seed=1)",
  "  -v[={-|+}]:詳しい情報を出力する。(-v:-v=+,オプション省略時:-v=-)"
};

const string Title = "=== シミュレーションによる検出力の計算 ==="; // 表題
constexpr int Gdef = 3;         // 省略時群数
constexpr double Mdef = 0.0;    // 省略時群の平均値
constexpr int Ndef = 100;       // 省略時群の例数
constexpr int Rdef = 10000;     // 省略時繰り返し数
constexpr double SDdef = 1.0;   // 省略時群の標準偏差
constexpr int SEEDdef = 1;      // 省略時乱数のタネ

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_a;                   // 有意水準
  bool v_am;                    // ANOVA&多重比較フラッグ
  int v_g;                      // 群数
  lstr mlst;                    // 平均値リスト
  vecdbl mvec;                  // 平均値ベクトル
  lstr nlst;                    // 例数リスト
  vecint nvec;                  // 例数ベクトル
  int v_r;                      // 繰り返し数
  lstr sdlst;                   // 標準偏差リスト
  vecdbl sdvec;                 // 標準偏差ベクトル
  int v_seed;                   // 乱数のタネ
  bool v_v;                     // 詳しい情報の出力フラッグ
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  double a() const              // 有意水準を返す
  { return v_a; }
  bool am() const               // ANOVA&多重比較フラッグを返す
  { return v_am; }
  double m(size_t idx) const    // 検出群の平均値を返す
  { return (idx < mvec.size()) ? mvec[idx] : 0.0; }
  int n(size_t idx) const       // 群の例数を返す
  { return (idx < nvec.size()) ? nvec[idx] : 0; }
  int g() const                 // 群数を返す
  { return v_g; }
  int r() const                 // 繰り返し数を返す
  { return v_r; }
  double sd(size_t idx) const   // 群の標準偏差を返す
  { return (idx < sdvec.size()) ? sdvec[idx] : 0.0; }
  int seed() const              // 乱数のタネを返す
  { return v_seed; }
  bool v() const                // 詳しい情報出力フラッグを返す
  { return v_v; }
  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_a(ALdef), v_am(false), v_g(Gdef), mlst(), mvec(), nlst(), nvec(), v_r(Rdef), sdlst(), sdvec(), v_seed(SEEDdef), v_v(false)
// コンストラクタ
{
  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] == "a")            // 有意水準の指定
  {
    v_a = opt[1].empty() ? ALdef : strtodf(opt[1]);
    if (v_a <= 0 || 1.0 <= v_a) 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] == "am")      // ANOVA&多重比較指定
  {
    if (opt[1].empty() || opt[1] == "+") v_am = true;
    else if (opt[1] == "-") v_am = false;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "g")      // 群数の指定
  {
    char* ptr = nullptr;
    v_g = opt[1].empty() ? Gdef : strtol(opt[1].c_str(), &ptr, 0);
    if (v_g < 2) err = badopt(pname(), optstr, USAGE, 2.0, (double)INT_MAX);
  }
  else if (opt[0] == "m")       // 群の平均値の指定
  {
    if (strsplit(mlst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "n")       // 群の例数の指定
  {
    if (strsplit(nlst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "r")       // 繰り返し数の指定
  {
    char* ptr = nullptr;
    v_r = opt[1].empty() ? Rdef : strtol(opt[1].c_str(), &ptr, 0);
    if (v_r < 1) err = badopt(pname(), optstr, USAGE, 1.0, (double)INT_MAX);
  }
  else if (opt[0] == "sd")      // 群の標準偏差の指定
  {
    if (strsplit(sdlst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "seed")    // 乱数のタネの指定
  {
    char* ptr = nullptr;
    v_seed = opt[1].empty() ? SEEDdef : strtol(opt[1].c_str(), &ptr, 0);
    if (v_seed < 1) err = badopt(pname(), optstr, USAGE, 1.0, (double)RAND_MAX);
  }
  else if (opt[0] == "v")       // 詳しい情報出力指定
  {
    if (opt[1].empty() || opt[1] == "+") v_v = true;
    else if (opt[1] == "-") v_v = 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);
}

int option::check()             // オプションのチェックとセット
{
  mvec.resize(v_g, Mdef);  nvec.resize(v_g, Ndef);  sdvec.resize(v_g, SDdef);
  if (mlst.size() > 0)
  {
    int i = 0;
    lstr_citr pos = mlst.begin();
    while (pos != mlst.end() && i < v_g)
    {
      mvec[i] = strtodf(*pos);
      if (mvec[i] == BLANKf) return disperr(pname(), "平均値指定が不適切！");
      else { ++i;  ++pos; }
    }
  }
  if (nlst.size() > 0)
  {
    int i = 0;
    lstr_citr pos = nlst.begin();
    while (pos != nlst.end() && i < v_g)
    {
      char* ptr = 0;
      nvec[i] = strtol(pos->c_str(), &ptr, 0);
      if (nvec[i] < 1) return disperr(pname(), "例数指定が不適切！");
      else { ++i;  ++pos; }
    }
  }
  if (sdlst.size() > 0)
  {
    int i = 0;
    lstr_citr pos = sdlst.begin();
    while (pos != sdlst.end() && i < v_g)
    {
      sdvec[i] = strtodf(*pos);
      if (sgn(sdvec[i]) <= 0 || sdvec[i] == BLANKf) return disperr(pname(), "標準偏差指定が不適切！");
      else { ++i;  ++pos; }
    }
  }
  return 0;                     // エラーなし
}

class varclass                  // 変数クラス
{
private:
  int v_g;                   // 群数
  int v_nc;                  // 多重比較の比較回数
  int v_nasig;               // 分散分析で有意になった回数
  int v_nfsig;               // Fisher型多重比較で有意になった回数
  int v_nbsig;               // Bonferroni型多重比較で有意になった回数
  int v_nssig;               // Scheffe型多重比較で有意になった回数
  vector<statk<double>> stativec; // 1回ごとの統計量
  vector<statk<double>> stattvec; // 標本平均値の統計量
  mat<double> v_smat;             // t値の相関計算用行列
  int v_nt;                       // t値の例数
  vecdbl tvec;                    // t値ベクトル
  vecdbl v_mevec;                 // t値の平均値ベクトル
  vecdbl v_vrvec;                 // t値の分散ベクトル
  void set_stat(const option&opt);
  // 正規乱数を発生して基礎統計量を集計する
  double cal_anova(double& vr) const; // 分散分析を行う
  bool cal_mcomp(bool& fsig, bool& bsig, bool& ssig, double vr, double alpha);
  // 多重比較を行う
public:
  varclass(const option& opt);  // コンストラクタ
  ~varclass()                   // デストラクタ
  {}
  void set(const option& opt);  // 変数をセットする
  int g() const                 // 群数を返す
  { return v_g; }
  int nc() const                // 多重比較の比較回数を返す
  { return v_nc; }
  int n(size_t idx) const       // 標本平均値の例数を返す
  { return (idx < stattvec.size()) ? stattvec[idx].n() : 0.0; }
  double mean(size_t idx) const // 標本平均値の平均値を返す
  { return (idx < stattvec.size()) ? stattvec[idx].mean() : 0.0; }
  double vr(size_t idx) const   // 標本平均値の不偏分散を返す
  { return (idx < stattvec.size()) ? stattvec[idx].vr() : 0.0; }
  double b1popu(size_t idx) const; // 標本平均値の母歪度を返す
  double b2popu(size_t idx) const; // 標本平均値の母尖度を返す
  int nt() const                // t値の例数を返す
  { return v_nt; }
  const mat<double> smat() const // t値の相関計算用行列を返す
  { return v_smat; }
  const vecdbl& mevec() const   // t値の平均値ベクトルを返す
  { return v_mevec; }
  const vecdbl& vrvec() const   // t値の分散ベクトルを返す
  { return v_vrvec; }
  int nasig() const             // 分散分析で有意になった回数を返す
  { return v_nasig; }
  int nfsig() const             // Fisher型多重比較で有意になった回数を返す
  { return v_nfsig; }
  int nbsig() const             // Bonferroni型多重比較で有意になった回数を返す
  { return v_nbsig; }
  int nssig() const             // Scheffe型多重比較で有意になった回数を返す
  { return v_nssig; }
};

varclass::varclass(const option&opt) : v_g(0), v_nasig(0), v_nfsig(0), v_nbsig(0), v_nssig(0), stativec(), stattvec(), v_smat(), v_nt(0), tvec(), v_mevec(), v_vrvec()
// コンストラクタ
{
  v_g = opt.g();
  v_nc = v_g * (v_g - 1) / 2;
  stativec.resize(v_g + 1);  stattvec.resize(v_g);
  v_smat.resize(v_nc + 1, v_nc + 1);
  tvec.resize(v_nc + 1, 0.0);
  v_mevec.resize(v_nc + 1, 0.0);  v_vrvec.resize(v_nc + 1, 0.0);
}

void varclass::set_stat(const option&opt)
// 正規乱数を発生して基礎統計量を集計する
{
  for (int j = 0; j < g() + 1; ++j) stativec[j].clear();
  for (int j = 0; j < g(); ++j)
  {
    int nj = opt.n(j);
    double mj = opt.m(j), sdj = opt.sd(j);
    for (int k = 0; k < nj; ++k) stativec[j].sigma(randn(mj, sdj));
  }
  for (int j = 0; j < g(); ++j) stativec[g()] += stativec[j];
}

double varclass::cal_anova(double& vr) const // 分散分析を行う
{
  vr = 0.0;
  double sa = 0.0;
  for (int j = 0; j < g(); ++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 = FDeps * fmax(FDmin, st);
  if (sa < Seps || sr < Seps || dfr < 1) return -1.0; // 検定できず
  vr = sr / dfr;
  double f = (sa / dfa) / vr;
  return ftop(f, dfa, dfr);     // 有意確率を返す
}

bool varclass::cal_mcomp(bool& fsig, bool& bsig, bool& ssig, double vr, double alpha)
// 多重比較を行う
{
  int nt = stativec[g()].n();
  int dfa = g() - 1, dfr = nt - g();
  fsig = false;  bsig = false;  ssig = false;
  tvec[0] = 1.0;
  int i = 1;
  for (int j = 0; 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;
      double f = dm * dm / (1.0 / (double)nj + 1.0 / (double)nk) / vr;
      tvec[i] = sgn(dm) * sqrt(f);
      ++i;
      double pr = ftop(f, 1, dfr);
      if (j == 0 && k == 1)     // Fisher型は1度だけ検定
      {
        if (fsig == false && pr < alpha) fsig = true;
      }
      pr = 1.0 - pow(1.0 - pr, nc()); // Bonferroni型
      if (bsig == false && pr < alpha) bsig = true;
      f = dm * dm / dfa / (1.0 / (double)nj + 1.0 / (double)nk) / vr;
      pr = ftop(f, dfa, dfr);   // Scheffe型
      if (ssig == false && pr < alpha) ssig = true;
    }
  }
  return fsig || bsig || ssig;  // どれか1でもtrueならtrue
}

void varclass::set(const option& opt) // 変数をセットする
{
  srand(opt.seed());
  for (int i = 0; i < opt.r(); ++i)
  {
    set_stat(opt);
    for (int j = 0; j < g(); ++j) stattvec[j].sigma(stativec[j].mean());
    double vr = 0.0;
    double pr = cal_anova(vr);
    if (pr < 0.0) continue;     // 検定できず
    if (pr < opt.a()) ++v_nasig;
    else if (opt.am()) continue; // ANOVAが有意の時だけ多重比較実施
    bool fsig = false, bsig = false, ssig = false;
    cal_mcomp(fsig, bsig, ssig, vr, opt.a());
    if (fsig) ++v_nfsig;
    if (bsig) ++v_nbsig;
    if (ssig) ++v_nssig;
    for (int i = 0; i <= nc(); ++i)
    {
      double ti = tvec[i];
      for (int j = 0; j <= nc(); ++j) v_smat.add(i, j, ti * tvec[j]);
    }
  }
  v_nt = (int)(v_smat(0, 0));
  for (int i = 0; i < nc() - 1; ++i)
  {
    for (int j = i + 1; j < nc(); ++j) v_smat.assign(j, i, v_smat(i, j));
  }
  v_smat.sweep(0);
  for (int i = 0; i <= nc(); ++i)
  {
    v_mevec[i] = v_smat(0, i);
    v_vrvec[i] = (v_nt > 1) ? v_smat(i, i) / (double)(v_nt - 1) : 0.0;
  }
}

double varclass::b1popu(size_t idx) const // 標本平均値の母歪度を返す
{
  if (idx >= stattvec.size()) return 0.0;
  if (sgn(stattvec[idx].m2()) <= 0) return 0.0;
  return stattvec[idx].m3() / stattvec[idx].m2() / sqrt(stattvec[idx].m2());
}

double varclass::b2popu(size_t idx) const // 標本平均値の母尖度を返す
{
  if (idx >= stattvec.size()) return 0.0;
  if (sgn(stattvec[idx].m2()) <= 0) return 0.0;
  return stattvec[idx].m4() / stattvec[idx].m2() / stattvec[idx].m2();
}

class powsimclass               // powsimメインクラス
{
private:
  void put_title(ostream& fout, const option& opt); // タイトルを出力する
  void set_sk(string& skew, string& kur, double b1, double b2) const;
  // 母歪度と母尖度を判定して文字列を返す
public:
  powsimclass()                 // コンストラクタ
  {}
  ~powsimclass()                // デストラクタ
  {}
  void put(option& opt);        // シミュレーション結果を出力する
};

void powsimclass::put_title(ostream& fout, const option& opt)
// タイトルを出力する
{
  puttitle(fout, Title);
  fout << "\n";
  fout << "○有意水準=" << opt.a() << " 群数=" << opt.g() << "  繰り返し数=" << opt.r() << "\n";
  if (opt.am()) fout << "  分散分析が有意の時だけ多重比較を実施する\n";
  else fout << "  分散分析と多重比較を独立に実施する\n";
}

void powsimclass::set_sk(string& skew, string& kur, double b1, double b2) const
// 母歪度と母尖度を判定して文字列を返す
{
  skew.erase();  kur.erase();
  if (b1 < -1.0) skew = "右傾斜";
  else if (fabs(b1) <= 1.0) skew = "ほぼ対称";
  else skew = "左傾斜";
  if (b2 < 0.0) kur = "鈍峰";
  else if (0.0 <= b2 && b2 <= 6.0) kur = "ほぼ正規";
  else kur = "鋭峰";
}

void powsimclass::put(option& opt) // シミュレーション結果を出力する
{
  if (opt.check() != 0) throw option_error("実行時オプションエラー！");
  varclass var(opt);            // 変数オブジェクト
  var.set(opt);                 // 変数をセットする
  put_title(cout, opt);
  cout << "\n";
  for (int j = 0; j < opt.g(); ++j)
  {
    cout << "・群" << j << ":" << "例数=" << opt.n(j) << " 平均値=" << opt.m(j) << " 標準偏差=" << opt.sd(j);
    if (opt.v())
    {
      cout << "  ↓標本平均の基礎統計量\n";
      cout << "    ";
      putmsd(cout, var.n(j), var.mean(j), var.vr(j));
      string skew, kur;
      set_sk(skew, kur, var.b1popu(j), var.b2popu(j));
      cout << "    歪度√β1=" << var.b1popu(j) << " " << skew << "  尖度β2=" << var.b2popu(j) << " " << kur;
    }
    cout << "\n";
  }
  if (opt.v())                  // t値の基礎統計量と相関係数出力
  {
    cout << "\n";
    putmsdvec(cout, var.nt(), var.mevec(), var.vrvec());
    putmat(cout, "・t値の相関行列", matr, var.nt(), var.smat());
  }
  double po = (double)var.nasig() / (double)(opt.r());
  cout << "\n・一元配置分散分析の検出力    =" << po << "(" << var.nasig() << "/" << opt.r() <<", β=" << 1.0 - po << ")\n";
  po = (double)var.nfsig() / (double)(opt.r());
  cout << "  Fisher型多重比較の検出力    =" << po << "(" << var.nfsig() << "/" << opt.r() <<", β=" << 1.0 - po << ")\n";
  po = (double)var.nbsig() / (double)(opt.r());
  cout << "  Bonferroni型多重比較の検出力=" << po << "(" << var.nbsig() << "/" << opt.r() <<", β=" << 1.0 - po << ")\n";
  po = (double)var.nssig() / (double)(opt.r());
  cout << "  Scheffe型多重比較の検出力   =" << po << "(" << var.nssig() << "/" << opt.r() <<", β=" << 1.0 - po << ")\n";
  cout << "\n" << endl;
}

int main(int argc, char* argv[]) // メインルーチン
{
  powsimclass powsim;           // メインクラス変数
  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);
    }
    powsim.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;
}
