// cltsim.C  Version 7.43  Copyleft (c) 杉本典夫
// Time-stamp: <2025-11-01 11:37:10 tac>
//
// 中心極限定理(CLT)のシミュレーション
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法  :g++ -O -s -W cltsim cltsim.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 =            // プログラム使用法
{
  "cltsim  Version 7.43         Copyleft (Vo\\oV) 2025 杉本典夫",
  "書式:cltsim  [<オプション> …]  [<データファイル> …]",
  "機能:中心極限定理のシミュレーション。",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定する。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -d[={i|ln|nr|r}]:データの種類を指定する。(オプション省略時:-d=nr)",
  "    i:整数  ln:対数正規乱数  nr:正規乱数  r:標準化乱数",
  "  -dnsname[=<出力ファイル名>]:DANSデータを出力するファイル名を指定する。",
  "    オプション省略時:DANSデータ出力せず  -dnsname:-dnsname=cltsim_dns.d",
  "  -f=<ファイル名>:指定されたファイルからオプションを読み込む。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示する。(-h:-h=man)",
  "  -m[={mean|med}]:要約値の種類を指定する。(オプション省略時:-m=mean)",
  "    mean:平均値  med:中央値",
  "  -ngp={g1|g2|mean|se}[,…]:Ngraph用データを指定する。",
  "    g1:歪度  g2:尖度  mean:平均値  se:標準誤差  オプション省略時:全て指定",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データを出力するファイル名を指定する。",
  "    オプション省略時:データ出力せず  -ngpname:-ngpname=cltsim_ngp.dat",
  "  -np[=<母集団の例数>]:母集団の例数を指定する。",
  "    1＜例数  オプション省略時:-np=100000",
  "  -ns[=<標本集団の例数>]:標本集団の例数を指定する。",
  "    0＜例数≦np  オプション省略時:-ns=100",
  "  -r[=<繰り返し数>]:繰り返し数を指定する。",
  "    0＜繰り返し数  オプション省略時:-r=10000",
  "  -res=[-][{def|all|stat}][,…]:解析結果の出力項目を指定する。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    stat:各種統計量"
};

enum dtype                      // データの種類
{
  dnon,                         // 指定なし
  di,                           // 整数
  dln,                          // 対数正規乱数
  dnr,                          // 正規乱数
  dr,                           // 標準化乱数
  detc                          // その他
};

enum mtype                      // 要約値の種類
{
  mnon,                         // 指定なし
  mmean,                        // 平均値
  mmed,                         // 中央値
  metc                          // その他
};

enum ngptype                    // Ngraph用出力データの種類
{
  ngpnon,                       // 出力なし
  ngpdns,                       // データファイル形式
  ngpg1,                        // 歪度分布
  ngpg2,                        // 尖度分布
  ngpmean,                      // 標本要約値の平均値分布
  ngpse,                        // 標準誤差分布
  ngpetc                        // その他
};

const string Title = "=== 中心極限定理のシミュレーション ==="; // 表題
const string DNSNAMEdef = "cltsim_dns.d"; // 省略時DANSファイル名
const string NGPNAMEdef = "cltsim_ngp.dat"; // 省略時Ngraphファイル名
constexpr dtype Ddef = dnr;     // 省略時データの種類
constexpr mtype Mdef = mmean;   // 省略時要約値の種類
constexpr int NPdef = 100000;   // 省略時母集団の例数
constexpr int NSdef = 100;      // 省略時標本集団の例数
constexpr int Rdef = 10000;     // 省略時繰り返し数

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  dtype v_d;                    // データの種類
  string v_dnsname;             // DANSファイル名
  mtype v_m;                    // 要約値の種類
  int v_np;                     // 母集団の例数
  int v_ns;                     // 標本集団の例数
  int v_r;                      // 繰り返し数
  string v_ngpname;             // Ngraphファイル名
  lstr ngplst;                  // Ngraph出力リスト
  mapsb ngpmap;                 // Ngraph出力の有無マップ
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  double cc() const             // 信頼係数を返す
  { return v_cc; }
  dtype d() const               // データの種類を返す
  { return v_d; }
  mtype m() const               // 要約値の種類を返す
  { return v_m; }
  int np() const                // 母集団の例数を返す
  { return v_np; }
  int ns() const                // 標本集団の例数を返す
  { return v_ns; }
  int r() const                 // 繰り返し数を返す
  { return v_r; }
  const string& dnsname() const // DANSデータファイル名を返す
  { return v_dnsname; }
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  bool ngp(const string& code); // Ngraph出力フラッグを返す
  bool res(const string& code); // 分析結果出力フラッグを返す
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check();                  // オプションのチェックとセット
  int set_ngp();                // -ngpオプションのセット
  int set_res();                // -resオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_cc(CCdef), v_d(Ddef),v_m(Mdef), v_np(NPdef), v_ns(NSdef), v_r(Rdef), v_dnsname(), v_ngpname(), ngplst(), ngpmap(), reslst(), resmap()
// コンストラクタ
{
  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); // プログラム名削除
}

bool option::ngp(const string& code) // Ngraph出力フラッグを返す
{
  mapsb_citr ngppos = ngpmap.find(code);
  if (ngppos == ngpmap.end()) return false; // code無し！
  return ngpmap[code];                      // true/falseを返す
}

bool option::res(const string& code) // 分析結果出力フラッグを返す
{
  mapsb_citr respos = resmap.find(code);
  if (respos == resmap.end()) return false; // code無し！
  return resmap[code];                      // true/falseを返す
}

void option::set(const string& optstr, char devc)
// オプション指定によってオプションを設定する
{
  string opt[2];                // オプションのopt[0]:項目,opt[1]:内容
  if (getopt(opt, optstr, devc) < 1) return; // オプション指定なし
  int err = 0;
  if (opt[0] == "cc")           // 信頼係数の指定
  {
    v_cc = opt[1].empty() ? CCdef : strtodf(opt[1]);
    if (v_cc < 0.0 || 1.0 <= v_cc) err = badopt(pname(), optstr, USAGE, 0.0, 1.0);
  }
  else if (opt[0] == "d")       // データの種類の指定
  {
    if (opt[1].empty() || opt[1] == "nr") v_d = dnr; // 正規乱数
    else if (opt[1] == "i") v_d = di;   // 整数
    else if (opt[1] == "ln") v_d = dln; // 対数正規乱数
    else if (opt[1] == "r") v_d = dr;   // 標準化乱数
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "dnsname") // DANSファイル名指定
  {
    if (opt[1].empty()) v_dnsname = DNSNAMEdef;
    else v_dnsname = opt[1];
    if (v_dnsname.size() < 3) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "m")       // 要約値の種類の指定
  {
    if (opt[1].empty() || opt[1] == "mean") v_m = mmean;
    else if (opt[1] == "med") v_m = mmed; // 中央値
    else err = badopt(pname(), optstr, USAGE);
  }
  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] == "ngp")     // Ngraph出力項目指定
  {
    if (strsplit(ngplst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "ngpname")     // Ngraphファイル名指定
  {
    if (opt[1].empty()) v_ngpname = NGPNAMEdef;
    else v_ngpname = opt[1];
    if (v_ngpname.size() < 3) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "np")      // 母集団の例数の指定
  {
    char* ptr = nullptr;
    v_np = opt[1].empty() ? NPdef : strtol(opt[1].c_str(), &ptr, 0);
    if (v_np < 1) err = badopt(pname(), optstr, USAGE, 2.0, (double)INT_MAX);
  }
  else if (opt[0] == "ns")      // 標本集団の例数の指定
  {
    char* ptr = nullptr;
    v_ns = opt[1].empty() ? NSdef : strtol(opt[1].c_str(), &ptr, 0);
    if (v_ns < 1) err = badopt(pname(), optstr, USAGE, 1.0, (double)INT_MAX);
  }
  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] == "res")     // 分析結果の出力指定
  {
    if (strsplit(reslst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else                          // 間違ったオプション指定
  { err = nonopt(pname(), optstr, USAGE); }
  if (err > 0) throw help_except("ヘルプ表示");
  if (err < 0) throw option_error(string("オプション指定違反→") + optstr);
}

int option::check()             // オプションのチェックとセット
{
  if (v_np < v_ns) return disperr(pname(), "(標本集団の例数≦母集団の例数)違反！");
  if (set_ngp() < 0) return -2; // -ngpオプションのエラー
  if (set_res() < 0) return -3; // -resオプションのセット
  return 0;                     // エラーなし
}

int option::set_ngp()           // -ngpオプションのセット
{
  ngpmap["g1"] = ngpmap["g2"] = ngpmap["mean"] = ngpmap["se"] = true;
  if (ngplst.size() > 0)
  {
    if (v_ngpname == "") return disperr(pname(), "-ngpnameオプションの指定が必要！");
    ngpmap["g1"] = ngpmap["g2"] = ngpmap["mean"] = ngpmap["se"] = false;
    for (const string& ngpstr : ngplst)
    {
      mapsb_citr ngppos = ngpmap.find(ngpstr);
      if (ngppos == ngpmap.end())
      { return disperr(pname(), "-ngp=" + ngpstr + "オプション指定が不適切！"); }
      else
      { ngpmap[ngpstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["stat"] = true;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["stat"] = false;
    for (const string& resstr : reslst)
    {
      mapsb_citr respos = resmap.find(resstr);
      if (resstr == "-")        // 全て出力せず
      { resmap["res"] = false; }
      else if (resstr[0] == '-' && resstr.size() > 1) // 特定の項目を出力せず
      {
        respos = resmap.find(resstr.substr(1));
        if (respos == resmap.end())
        { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
        else
        { resmap[resstr.substr(1)] = false; }
      }
      else if (resstr == "def" || resstr == "all") // 全て出力
      { resmap["def"] = resmap["stat"] = true; }
      else if (respos == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

class popu                      // 母集団クラス
{
private:
  bool v_smp;                   // データ選択フラッグ
  double v_d;                   // データ
public:
  popu() : v_smp(false), v_d(0.0) // デフォールトコンストラクタ
  {}
  popu(bool smp0, double d0) : v_smp(smp0), v_d(d0) // コンストラクタ
  {}
  popu(const popu& popu0) : v_smp(popu0.v_smp), v_d(popu0.v_d)
                            // コピーコンストラクタ
  {}
  virtual ~popu()               // デストラクタ
  {}
  const popu& operator=(const popu& popu0); // コピー代入演算子
  bool operator<(const popu& popu0) const;  // 比較演算子
  double d() const                          // データを返す
  { return v_d; }
  double d(double d0)           // データを設定する
  { return (v_d = d0); }
  bool smp() const              // データ選択フラッグを返す
  { return v_smp; }
  bool smp(int smp0)            // データ選択フラッグを設定する
  { return (v_smp = smp0); }
};

const popu& popu::operator=(const popu& popu0) // コピー代入演算子
{
  if (&popu0 == this) return *this; // 自分自身への代入を防ぐ
  v_smp = popu0.v_smp;  v_d = popu0.v_d;
  return *this;                 // 自分自身への参照を返す
}

bool popu::operator<(const popu& popu0) const // 比較演算子
{
  return (d() < popu0.d());
}

using vecpopu = vector<popu>;   // 母集団データベクトルの短名定義

double cal_mv(vecpopu& tmpvec, statk<double>& statstmp, mtype mttmp)
// 要約値を返す
{
  if (mttmp == mmean) return statstmp.mean(); // 平均値を返す
  int n = statstmp.n();
  vecpopu tvec(tmpvec);         // 一時的なデータベクトル
  sort(tvec.begin(), tvec.end());
  int i = n / 2;                // 中央値
  double mv = ((n % 2) == 0) ? mv = (tvec[i - 1].d() + tvec[i].d()) * 0.5 : tvec[i].d();
  return mv;                    // 中央値を返す
}

class varclass                  // 変数クラス
{
private:
  vecpopu popuvec;              // 母集団のデータベクトル
  statk<double> statspopu;      // 母集団の基礎統計量
  statk<double> statssmp;       // 標本要約値の基礎統計量
  vecdbl msmpvec;              // 標本要約値ベクトル
  vecdbl zsmpvec;               // 標準化標本要約値ベクトル
  vecdbl zlsmpvec;              // 非心標準化要約値ベクトル
  ftable v_ftabps;              // 母集団と標本集団のクロス集計表
  double v_mpopu;               // 母数
  double v_mpl;                 // 母集団の信頼限界下限
  double v_mpu;                 // 母集団の信頼限界上限
  double v_semp;                // 標本要約値の理論標準誤差
  double cal_msmp(double& vm, int ns, mtype mty, int seed);
  // 標本集団の要約値とその分散を求めて返す
public:
  varclass(dtype dy, int np, int r, mtype mty); // コンストラクタ
  ~varclass()                   // デストラクタ
  {}
  void set(int ns, mtype mty, double cc0); // 変数をセットする
  double dpopu(size_t idx) const // 母集団のデータを返す
  { return (idx < popuvec.size()) ? popuvec[idx].d() : BLANKf; }
  double msmp(size_t idx) const // 標本要約値を返す
  { return (idx < msmpvec.size()) ? msmpvec[idx] : BLANKf; }
  double zsmp(size_t idx) const // 標準化標本要約値を返す
  { return (idx < zsmpvec.size()) ? zsmpvec[idx] : BLANKf; }
  double zlsmp(size_t idx) const // 非心標準化標本要約値を返す
  { return (idx < zlsmpvec.size()) ? zlsmpvec[idx] : BLANKf; }
  const ftable& ftabps() const // 母集団と標本集団のクロス集計表を返す
  { return v_ftabps; }
  int npopu() const             // 母集団の例数を返す
  { return statspopu.n(); }
  double meanpopu() const       // 母平均値を返す
  { return statspopu.mean(); }
  double vrpopu() const         // 母分散を返す
  { return statspopu.vrs(); }
  double b1popu() const         // 母歪度を返す
  { return (sgn(statspopu.m2()) > 0) ? statspopu.m3() / statspopu.m2() / sqrt(statspopu.m2()) : 0.0; }
  double b2popu() const         // 母尖度を返す
  { return (sgn(statspopu.m2()) > 0) ? statspopu.m4() / statspopu.m2() / statspopu.m2() : 0.0; }
  int nsmp() const              // 標本集団の例数を返す
  { return statssmp.n(); }
  double meansmp() const        // 標本平均値を返す
  { return statssmp.mean(); }
  double vrsmp() const          // 不偏分散を返す
  { return statssmp.vr(); }
  double b1smp() const          // 標本歪度を返す
  { return (sgn(statssmp.k2()) > 0) ? statssmp.k3() / statssmp.k2() / sqrt(statssmp.k2()) : 0.0; }
  double b2smp() const          // 標本尖度を返す
  { return (sgn(statssmp.k2()) > 0) ? statssmp.k4() / statssmp.k2() / statssmp.k2() : 0.0; }
  double mpopu() const          // 母数を返す
  { return v_mpopu; }
  double mpl() const            // 出現率下限を返す
  { return v_mpl; }
  double mpu() const            // 出現率上限を返す
  { return v_mpu; }
  double semp() const           // 標本要約値の理論標準誤差を返す
  { return v_semp; }
};

varclass::varclass(dtype dty, int np, int r, mtype mty) : popuvec(np), statspopu(), statssmp(), msmpvec(r, 0.0), zsmpvec(r, 0.0), zlsmpvec(r, 0.0), v_ftabps(2, 2), v_mpopu(0.0), v_mpl(0.0), v_mpu(0.0), v_semp(0.0)
// コンストラクタ
{
  double d = 0.0, m = 0.0;
  statk<double> statstmp;       // 標準化用一時的基礎統計量
  switch (dty)                  // データの種類
  {
  case di :                     // 整数
    for (int i = 0; i < np; ++i)
    {
      d = (double)(i + 1);
      statspopu.sigma(d);  popuvec[i].d(d);
    }
    break;
  case dln : srand(np);         // 対数正規乱数
    for (int i = 0; i < np; ++i)
    {
      d = randn(0.0, 1.0);      // N(0,1)の正規乱数
      statstmp.sigma(d);  popuvec[i].d(d);
    }
    for (int i = 0; i < np; ++i) // 標準化と指数変換
    {
      d = exp((popuvec[i].d() - statstmp.mean()) / statstmp.sds());
      statspopu.sigma(d);  popuvec[i].d(d);
    }
    break;
  case dr : srand(np);          // 乱数
    for (int i = 0; i < np; ++i)
    {
      d = (double)(rand()) * 0.00000001;
      statstmp.sigma(d);  popuvec[i].d(d);
    }
    for (int i = 0; i < np; ++i) // 標準化
    {
      d = (popuvec[i].d() - statstmp.mean()) / statstmp.sds();
      statspopu.sigma(d);  popuvec[i].d(d);
    }
    break;
  default : srand(np);          // 正規乱数
    for (int i = 0; i < np; ++i)
    {
      d = randn(0.0, 1.0);      // N(0,1)の正規乱数
      statstmp.sigma(d);  popuvec[i].d(d);
    }
    for (int i = 0; i < np; ++i) // 標準化
    {
      d = (popuvec[i].d() - statstmp.mean()) / statstmp.sds();
      statspopu.sigma(d);  popuvec[i].d(d);
    }
  }
  v_mpopu = cal_mv(popuvec, statspopu, mty); // 母数
}

double varclass::cal_msmp(double& vm, int ns, mtype mty, int seed)
// 標本集団の要約値とその分散を求めて返す
{
  size_t np = popuvec.size();
  for (size_t i = 0; i < np; ++i) popuvec[i].smp(false);
  srand(seed);                  // 乱数の種をセットする
  vecpopu smpvec(ns);           // 標本集団のデータベクトル
  statk<double> statstmp;       // 標本集団の基礎統計量
  int i = 0;
  while (i < ns)
  {
    int idx = rand() % np;           // 乱数の下位np桁
    if (popuvec[idx].smp() == false) // 未選択データ
    {
      statstmp.sigma(popuvec[idx].d());
      smpvec[i].d(popuvec[idx].d());
      popuvec[idx].smp(true);  ++i;
    }
  }
  double msmp = cal_mv(smpvec, statstmp, mty);
  vm = statstmp.vr() / (double)ns; // 要約値の分散
  return msmp;                    // 要約値を返す
}

void varclass::set(int ns, mtype mty, double cc0) // 変数をセットする
{
  v_ftabps.hed("母集団＼標本集団");
  v_ftabps.ccode(1, "  信頼区間");  v_ftabps.ccode(2, "  棄却域");
  v_ftabps.rcode(1, "信頼区間");  v_ftabps.rcode(2, "棄却域");
  double cf = (double)(popuvec.size() - ns) / (double)(popuvec.size() - 1);
  if (mty == mmed) cf *= M_PI * 0.5; // 中央値の分散修正項
  v_semp = sqrt(statspopu.vrs() / (double)ns * cf); // 理論標準誤差
  double alpha = 1.0 - cc0;
  double ta = ptot(alpha, ns - 1), za = ptoz(alpha);
  double ci = za * v_semp;
  v_mpl = v_mpopu - ci;  v_mpu = v_mpopu + ci; // 母集団から求めた信頼区間
  double lamda = statspopu.sd() / sqrt((double)ns); // 非心度
  size_t nocip = 0, nocis = 0;
  for (size_t i = 0; i < msmpvec.size(); ++i)
  {
    double vm = 0.0;
    double msmp = cal_msmp(vm, ns, mty, i + 1);
    double sem = sqrt(vm);
    msmpvec[i] = msmp;
    double z = (sem > 0.0) ? (msmp - v_mpopu) / sem : 0.0;
    zsmpvec[i] = z;
    zlsmpvec[i] = (sem > 0.0) ? (msmp - v_mpopu + lamda) / sem : 0.0;
    statssmp.sigma(msmp);
    nocip = (v_mpl <= msmp && msmp <= v_mpu) ? 1 : 2; // 母集団の信頼区間チェック
    double ci = ta * sem;
    double ml = msmp - ci, mu = msmp + ci; // 標本集団から求めた信頼区間
    nocis = (ml <= v_mpopu && v_mpopu <= mu) ? 1 : 2; // 標本集団の信頼区間チェック
    v_ftabps.f(nocip, nocis, v_ftabps.f(nocip, nocis) + 1);
  }
  v_ftabps.set();               // 度数分布表の整理
}

class cltsimclass               // cltsimメインクラス
{
private:
  void put_stat(ostream& fout, option& opt, const varclass& var) const;
  // 母集団と標本要約値の計算結果を出力する
  void set_sk(string& skew, string& kur, double b1, double b2, bool popu) const;
  // 歪度と尖度を判定して文字列を返す
  void put_ftab(ostream& fout, option& opt, double mpopu, double mpl, double mpu, const ftable& ftabtmp) const;
  // 分割表と診断率を出力する
  void put_dns(option& opt, const varclass& var) const;
  // 母集団と標本要約値をデータファイル形式で出力する
  void plot_ngp(option& opt, const varclass& var) const;
  // Ngraph用各種データを出力する
  void plot_ngpid(ofstream& fout, option& opt) const;
  // Ngraph用データのタイトルを出力する
  void plot_ngpdata(ofstream& fout, option& opt, const varclass& var, ngptype ngp) const;
  // Ngraph用プロットデータを出力する
public:
  cltsimclass()                 // コンストラクタ
  {}
  ~cltsimclass()                // デストラクタ
  {}
  void put(option& opt);        // シミュレーション結果を出力する
};

void cltsimclass::put_stat(ostream& fout, option& opt, const varclass& var) const
// 母集団と標本要約値の計算結果を出力する
{
  puttitle(fout, Title);
  bool putres = false;
  if (opt.res("stat"))          // 各種統計量を出力する
  {
    fout << "\n・母集団のデータ:";
    switch (opt.d())
    {
    case di : fout << "整数";
      break;
    case dln : fout << "対数正規乱数";
      break;
    case dr : fout << "一様乱数";
      break;
    default : fout << "正規乱数";
    }
    fout << " 母数=" << var.mpopu() << "\n";
    putmsd(fout, var.npopu(), var.meanpopu(), var.vrpopu());
    double b1 = var.b1popu(), b2 = var.b2popu(); // 母尖度と母歪度
    string skew, kur;
    set_sk(skew, kur, b1, b2, true);
    fout << "歪度:√β1=" << b1 << " " << skew << "  尖度:β2=" << b2 << " " << kur << "\n";
    fout << "\n・標本要約値(m):";
    if (opt.m() == mmed) fout << "中央値";
    else fout << "平均値";
    fout << " 例数=" << opt.ns() << " 標準誤差理論値=" << var.semp() << " 繰り返し数=" << opt.r() << "\n";
    putmsd(fout, var.nsmp(), var.meansmp(), var.vrsmp());
    b1 = var.b1smp();  b2 = var.b2smp(); // 標本尖度と標本歪度
    set_sk(skew, kur, b1, b2, false);
    fout << "歪度:g1=" << b1 << " " << skew << "  尖度:g2=" << b2 << " " << kur << "\n";
    putres = true;
  }
  if (putres) fout << "\n";
  put_ftab(fout, opt, var.mpopu(), var.mpl(), var.mpu(), var.ftabps());
}

void cltsimclass::put_ftab(ostream& fout, option& opt, double mpopu, double mpl, double mpu, const ftable& ftabtmp) const
// 分割表と診断率を出力する
{
  int df = opt.ns() - 1;
  double pro = opt.cc() * 100.0;
  fout << "・母集団の" << pro << "%信頼区間:正規近似 ml=" << mpl << "〜mu=" << mpu << "\n";
  fout << "  標本集団の" << pro << "%信頼区間:t分布(自由度=" << df << ") 母数=" << mpopu << "が入った回数\n";
//  ftabtmp.set();
  putftable(fout, ftabtmp, tabfprt);
  int ncip = ftabtmp.frtotal(1), ncis = ftabtmp.fctotal(1);
  int nrejp = opt.r() - ncip, nrejs = opt.r() - ncis;
  pro = (ftabtmp.frtotal(2) > 0) ? (double)(ftabtmp.f(2, 2)) / (double)(ftabtmp.frtotal(2)) : 0.0;
  fout << "検定:感度=" << ftabtmp.f(2, 2) << "/" << ftabtmp.frtotal(2) << "=" << pro;
  pro = (ftabtmp.frtotal(1) > 0) ? (double)(ftabtmp.f(1, 1)) / (double)(ftabtmp.frtotal(1)) : 0.0;
  fout << " 特異度=" << ftabtmp.f(1, 1) << "/" << ftabtmp.frtotal(1) << "=" << pro;
  int nac = ftabtmp.f(1, 1) + ftabtmp.f(2, 2);
  pro = (double)nac / (double)(ftabtmp.ftotal());
  fout << " 正診率=" << nac << "/" << ftabtmp.ftotal() << "=" << pro << "\n";
}

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

void cltsimclass::put_dns(option& opt, const varclass& var) const
// 母集団と標本要約値をデータファイル形式で出力する
{
  ofstream fout;
  if (ofopen(fout, (opt.dnsname()).c_str()) == false) return; // ファイルオープン失敗
  constexpr int Nitem = 4;
  idd id("中心極限定理のシミュレーション", Nitem);
  ostringstream ostr;
  switch (opt.d())
  {
  case di : ostr << "整数";
    break;
  case dln : ostr << "対数正規乱数";
    break;
  case dr : ostr << "一様乱数";
    break;
  default : ostr << "正規乱数";
  }
  id.assign_item(0, "母集団のデータ:" + ostr.str());
  ostr.str("");
  ostr << opt.ns();
  if (opt.m() == mmed) id.assign_item(1, "中央値:例数=" + ostr.str());
  else id.assign_item(1, "平均値:例数=" + ostr.str());
  id.assign_item(2, "標準化データ:例数=" + ostr.str());
  id.assign_item(3, "非心標準化データ:例数=" + ostr.str());
  id.put(fout, formtag);
  datac cdata(Nitem);
  for (int i = 0; i < opt.np() || i < opt.r(); ++i)
  {
    cdata.assign(0, var.dpopu(i));  cdata.assign(1, var.msmp(i));
    cdata.assign(2, var.zsmp(i));  cdata.assign(3, var.zlsmp(i));
    cdata.put(fout);
  }
}

void cltsimclass::plot_ngp(option& opt, const varclass& var) const
// Ngraph用各種データを出力する
{
  ofstream fout;
  if (ofopen(fout, (opt.ngpname()).c_str()) == false) return; // ファイルオープン失敗
  bool delml = false;
  if (opt.ngp("mean"))          // 平均値
  {
    fout << "-s1 -dline -w80\n";
    plot_ngpid(fout, opt);
    plot_ngpdata(fout, opt, var, ngpmean);
    delml = true;
  }
  if (opt.ngp("se"))            // 標準誤差
  {
    if (delml)                  // 区切り行を出力する
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline -w80\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_ngpdata(fout, opt, var, ngpse);
  }
  if (opt.ngp("g1"))            // 歪度
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline -w80\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_ngpdata(fout, opt, var, ngpg1);
    delml = true;
  }
  if (opt.ngp("g2"))            // 尖度
  {
    if (delml) fout << NGPdelm << "\n";
    else fout << "-s1 -dline -w80\n";
    plot_ngpid(fout, opt);
    plot_ngpdata(fout, opt, var, ngpg2);
  }
}

void cltsimclass::plot_ngpid(ofstream& fout, option& opt) const
// Ngraph用データのタイトルを出力する
{
  fout << "%中心極限定理のシミュレーション\n";
  fout << "%母集団のデータ:";
  switch (opt.d())
  {
  case di : fout << "整数";
    break;
  case dln : fout << "対数正規乱数";
    break;
  case dr : fout << "一様乱数";
    break;
  default : fout << "正規乱数";
  }
  fout << "  例数=" << opt.np() << "\n";
  fout << "%標本集団の要約値:";
  string tstr;
  if (opt.m() == mmed) tstr = "中央値";
  else tstr = "平均値";
  fout << tstr << "  例数=" << opt.ns() << "  繰り返し数=" << opt.r() << "\n";
}

void cltsimclass::plot_ngpdata(ofstream& fout, option& opt, const varclass& var, ngptype ngp) const
// Ngraph用プロットデータを出力する
{
  constexpr int RDdef = 50;     // グラフの描画範囲分割数省略時値
  switch(ngp)
  {
  case ngpg1 : fout << "%要約値の歪度g1";
    break;
  case ngpg2 : fout << "%要約値の尖度g2";
    break;
  case ngpse : fout << "%要約値の標準誤差";
    break;
  default : fout << "%要約値の平均値";
  }
  fout << "\n";
  int rd = opt.r() / RDdef;     // rdデータごとに出力する
  statk<double> statstmp;
  for (int i = 0; i < opt.r(); ++i)
  {
    int x = i + 1;
    statstmp.sigma(var.msmp(i));
    if (rd < 2 || (x % rd) == 0)
    {
      double b1 = (sgn(statstmp.k2()) > 0) ? statstmp.k3() / statstmp.k2() / sqrt(statstmp.k2()) : 0.0;
      double b2 = (sgn(statstmp.k2()) > 0) ? statstmp.k4() / statstmp.k2() / statstmp.k2() : 0.0;
      fout << x << " ";
      switch(ngp)
      {
      case ngpg1 : fout << b1;
        break;
      case ngpg2 : fout << b2;
        break;
      case ngpse : fout << statstmp.sd();
        break;
      default : fout << statstmp.mean();
      }
      fout << "\n";
    }
  }
}

void cltsimclass::put(option& opt) // シミュレーション結果を出力する
{
  if (opt.check() != 0) throw option_error("実行時オプションエラー！");
  varclass var(opt.d(), opt.np(), opt.r(), opt.m()); // 変数オブジェクト
  var.set(opt.ns(), opt.m(), opt.cc());              // 変数をセットする
  if (opt.res("res"))           // 計算結果出力
  {
    put_stat(cout, opt, var);
    cout << "\n" << endl;
  }
  if (opt.dnsname() != "") put_dns(opt, var); // DNAデータ出力
  if (opt.ngpname() != "") plot_ngp(opt, var); // Ngraph用データ出力
}

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