// statist.C Version 7.45  Copyleft (c) 杉本典夫
// Time-stamp: <2026-01-03 16:37:56 tac>
//
// 実数をデータを集計し、基礎統計量を出力する。複数項目を一度に集計可能。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -w statist statist.C
//               -I/usr/local/include -L/usr/local/lib -ldans -lfun

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

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

const vecstr USAGE =            // プログラム使用法表示用
{
  "statist  Version 7.45        Copyleft(Vo\\oV) 2026 杉本典夫",
  "書式:statist  [<オプション> …]  [<データファイル> …]",
  "機能:各項目の基礎統計量を計算するのだ。",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定するのだ。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -ch=<階級下限>,<級間>:度数分布データの階級下限値と級間を指定する。",
  "  -dnsname[=<出力ファイル名>]:DANSデータを出力するファイル名を指定する。",
  "    オプション省略時:DANSデータ出力せず  -dnsname:-dnsname=statist_dns.d",
  "    出力項目:項目名と基礎統計量",
  "  -f=<ファイル名>:指定されたファイルからオプションを読み込む。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示する。(-h:-h=man)",
  "  -ngp={box|cd|cdp|d|ncd|nd|nrcd|nrd|rcd|rcdp|rd|vnrd|vrd}[,…]",
  "      :Ngraph用データの指定。",
  "    box:箱ヒゲプロット  cd:累積度数分布  cdp:累積度数分布プロット  d:度数分布",
  "    ncd:理論累積度数分布  nd:理論度数分布  nrcd:理論相対累積度数分布",
  "    nrd:理論相対度数分布  rcd:相対累積度数分布  rcdp:相対累積度数分布プロット",
  "    rd:相対度数分布  vnrd:仮想的理論相対度数分布  vrd:仮想的相対度数分布",
  "    オプション省略時:全て指定",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データを出力するファイル名を指定する。",
  "    オプション省略時:Ngraph用データ出力せず  -ngpname:-ngpname=statist_ngp.dat",
  "  -pt=<％tile>[,…]:パーセンタイルを出力する。(0＜％tile≦100)",
  "  -res=[-][{def|all|norm|stat}][,…]:解析結果の出力項目を指定する。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    norm:正規性の指標  stat:基礎統計量",
  "  -sep[={-|+}]:項目の集計法を指定する。(-sep:-sep=+)"
};

const string Title = "=== 基礎統計量 ==="; // 表題
constexpr uint RDdef = 50;      // グラフの描画範囲分割数省略時値
const string DNSNAMEdef = "statist_dns.d"; // 省略時DANSファイル名
const string NGPNAMEdef = "statist_ngp.dat"; // 省略時Ngraphファイル名
const double Ddef = 0.1;        // 仮想的相対度数分布の微分係数

using listdbl = list<double>;   // 倍精度実数リストの別名定義

enum ngptype                    // Ngraph用出力データの種類
{
  ngpnon,                       // 出力なし
  ngpbox,                       // 箱ヒゲプロット
  ngpcd,                        // 累積度数分布
  ngpcdp,                       // 累積度数分布のプロット
  ngpd,                         // 度数分布
  ngpdns,                       // データファイル形式
  ngpncd,                       // 理論累積度数分布
  ngpnd,                        // 理論度数分布
  ngpnrcd,                      // 理論相対累積度数分布
  ngpnrd,                       // 理論相対度数分布
  ngprcd,                       // 相対累積度数分布
  ngprcdp,                      // 相対累積度数分布のプロット
  ngprd,                        // 相対度数分布
  ngpvnrd,                      // 仮想的理論相対度数分布
  ngpvrd,                       // 仮想的相対度数分布
  ngpetc                        // その他
};

class option                    // オプション集合体
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  double v_cmin;                // 度数分布の階級下限
  double v_h;                   // 度数分布の級間
  string v_dnsname;             // DANSファイル名
  string v_ngpname;             // Ngraphファイル名
  lstr ngplst;                  // Ngraph出力リスト
  mapsb ngpmap;                 // Ngraph出力の有無マップ
  lstr ptlst;                   // パーセンタイルリスト
  vecdbl ptvec;                 // パーセンタイルベクトル
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
  bool v_sep;                   // 項目ごと集計フラッグ
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; }
  double cmin() const           // 度数分布の階級下限を返す
  { return v_cmin; }
  double h() const              // 度数分布の級間を返す
  { return v_h; }
  const string& dnsname() const // DANSデータファイル名を返す
  { return v_dnsname; }
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  bool ngp(const string& code); // Ngraph出力フラッグを返す
  size_t ptsize() const         // パーセンタイル数を返す
  { return ptlst.size(); }
  double pt(size_t idx) const   // パーセンタイルを返す
  { return (idx < ptvec.size()) ? ptvec[idx] : BLANKf; }
  bool res(const string& code); // 分析結果出力フラッグを返す
  bool sep() const              // 項目ごと集計フラッグを返す
  { return v_sep; }
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check();                  // オプションのセットとチェック
  int set_pt();                 // -ptオプションのセット
  int set_ngp();                // -ngpオプションのセット
  int set_res();                // -resオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_cc(CCdef), v_cmin(BLANKf), v_h(BLANKf), v_dnsname(), v_ngpname(), ngplst(), ngpmap(), ptlst(), ptvec(), reslst(), resmap(), v_sep(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); // プログラム名削除
}

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] == "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] == "f")       // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "ch")      // 階級下限と級間の指定
  {
    datac cdata(2);
    cdata.getdata(opt[1]);
    v_cmin = cdata.fdata(0);  v_h = cdata.fdata(1);
    if (v_cmin == BLANKf || v_h == BLANKf) err = badopt(pname(), optstr, USAGE);
  }
  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] == "pt")      // パーセンタイルの指定
  {
    if (strsplit(ptlst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "res")     // 分析結果の出力指定
  {
    if (strsplit(reslst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "sep")     // 項目ごとの集計指定
  {
    if (opt[1].empty() || opt[1] == "+") v_sep = true;
    else if (opt[1] == "-") v_sep = 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()             // オプションのセットとチェック
{
  if (set_pt() < 0) return -1;  // -ptオプションのセット
  if (set_ngp() < 0) return -2; // -ngpオプションのセット
  if (set_res() < 0) return -3; // -resオプションのセット
  return 0;                     // エラーなし
}

int option::set_pt()            // -ptオプションのセット
{
  if (ptlst.empty()) return 0;  // -ptオプション非指定時
  ptvec.resize(ptlst.size());
  size_t i = 0;
  for (auto ptstr : ptlst)
  {
    double pti = strtodf(ptstr);
    if (pti <= 0.0 || 100.0 < pti) return disperr(pname(), "-pt=" + ptstr + "←パーセンタイルが不適切！");
    if (i > 0 && pti <= ptvec[i - 1]) return disperr(pname(), "-pt=" + ptstr + "パーセンタイルが前の値以下！");
    ptvec[i] = pti;
    ++i;
  }
  return 0;                     // エラーなし
}

int option::set_ngp()           // -ngpオプションのセット
{
  ngpmap["box"] = ngpmap["cd"] = ngpmap["cdp"] = ngpmap["d"] = ngpmap["ncd"] = ngpmap["nd"] = ngpmap["nrcd"] = ngpmap["nrd"] = ngpmap["rcd"] = ngpmap["rcdp"] = ngpmap["rd"] = ngpmap["vrd"] = ngpmap["vnrd"] = true;
  if (ngplst.size() > 0)
  {
    if (v_ngpname == "") return disperr(pname(), "-ngpnameオプションの指定が必要！");
    ngpmap["box"] = ngpmap["cd"] = ngpmap["cdp"] = ngpmap["d"] = ngpmap["ncd"] = ngpmap["nd"] = ngpmap["nrcd"] = ngpmap["nrd"] = ngpmap["rcd"] = ngpmap["rcdp"] = ngpmap["rd"] = ngpmap["vrd"] = ngpmap["vnrd"] = 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["norm"] = resmap["stat"] = true;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["norm"] = 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["norm"] = resmap["stat"] = true; }
      else if (respos == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

class stats                     // 基礎統計量+k統計量+四分位値クラス
{
private:
  statk<double> v_fdstatk;      // 基礎統計量+k統計量クラス
  double v_min;                 // 最小値
  double v_q1;                  // 第1四分位点(25パーセンタイル)
  double v_med;                 // 第1四分位点(中央値)
  double v_q3;                  // 第3四分位点(75パーセンタイル)
  double v_max;                 // 最大値
  listdbl v_fdlst;              // 実数データリスト
public:
  stats() : v_fdstatk(), v_min(BLANKf), v_q1(BLANKf), v_med(BLANKf), v_q3(BLANKf), v_max(BLANKf), v_fdlst()
  // デフォールトコンストラクタ
  {}
  stats(int n0, double sig0, double sig20, double sig30, double sig40, double min0, double q10, double med0, double q30, double max0, listdbl& fdlst0) : v_fdstatk(n0, sig0, sig20, sig30, sig40), v_min(min0), v_q1(q10), v_med(med0), v_q3(q30), v_max(max0), v_fdlst(fdlst0)
  // コンストラクタ
  {}
  stats(const stats& stats0) : v_fdstatk(stats0.v_fdstatk), v_min(stats0.v_min), v_q1(stats0.v_q1), v_med(stats0.v_med), v_q3(stats0.v_q3), v_max(stats0.v_max), v_fdlst(stats0.v_fdlst)
  // コピーコンストラクタ
  {}
  virtual ~stats()              // デストラクタ
  {}
  const stats& operator=(const stats& stats0); // コピー代入演算子
  void clear();                                // データをクリアする
  int sigma(double x);                         // データの累積集計
  int sort();                                  // データをソートする
  const statk<double>& fdstatk() const // 基礎統計量+k統計量クラスを返す
  { return v_fdstatk; }
  const listdbl& fdlst() const  // 実数データリストを返す
  { return v_fdlst; }
  int n() const                 // 例数を返す
  { return v_fdstatk.n(); }
  double mean() const           // 平均値を返す
  { return v_fdstatk.mean(); }
  double vr() const             // 不偏分散を返す
  { return v_fdstatk.vr(); }
  double sd() const             // 不偏標準偏差を返す
  { return v_fdstatk.sd(); }
  double se() const             // 標準誤差を返す
  { return v_fdstatk.se(); }
  double k2() const             // k2統計量を返す
  { return v_fdstatk.k2(); }
  double k3() const             // k3統計量を返す
  { return v_fdstatk.k3(); }
  double k4() const             // k4統計量を返す
  { return v_fdstatk.k4(); }
  double min() const            // 最小値を返す
  { return v_min; }
  double q1() const             // 第1四分位値を返す
  { return v_q1; }
  double med() const            // 中央値を返す
  { return v_med; }
  double q3() const             // 第3四分位値を返す
  { return v_q3; }
  double max() const            // 最大値を返す
  { return v_max; }
  double ptv(double pt);        // パーセンタイル値を返す
};

const stats& stats::operator=(const stats& stats0) // コピー代入演算子
{
  if (&stats0 == this) return *this; // 自分自身への代入を防ぐ
  v_fdstatk = stats0.v_fdstatk;
  v_min = stats0.v_min;  v_q1 = stats0.v_q1;  v_med = stats0.v_med;
  v_q3 = stats0.v_q3;  v_max = stats0.v_max;  v_fdlst = stats0.v_fdlst;
  return *this;
}

void stats::clear()             // データをクリアする
{
  v_fdstatk.clear();
  v_min = v_q1 = v_med = v_q3 = v_max = BLANKf;  v_fdlst.clear();
}

int stats::sigma(double x)      // データの累積集計
{
  if (x == BLANKf) return 0;    // 累積せず
  v_fdstatk.sigma(x);  v_fdlst.push_back(x);
  return 1;                     // 累積
}

int stats::sort()               // データをソートする
{
  if (v_fdlst.empty()) return 0; // データ無し
  if (v_min != BLANKf) return 1; // ソート済み
  v_fdlst.sort();
  listdbl::const_iterator fditr = v_fdlst.begin();
  v_min = v_q1 = v_med = *fditr;
  fditr = v_fdlst.end();
  --fditr;
  v_max = v_q3 = *fditr;
  if (v_fdlst.size() < 2) return 1;
  double fidx = (double)(v_fdlst.size() - 1) * 0.25;
  uint idx = (uint)fidx;
  fditr = v_fdlst.begin();
  for (uint i = 0; i < idx; ++i) ++fditr;
  double xpre = *fditr;
  ++fditr;
  v_q1 = xpre + (*fditr - xpre) * (fidx - (double)idx);
  fidx = (double)(v_fdlst.size() - 1) * 0.5;
  idx = (uint)(fidx);
  fditr = v_fdlst.begin();
  for (uint i = 0; i < idx; ++i) ++fditr;
  xpre = *fditr;
  ++fditr;
  v_med = xpre + (*fditr - xpre) * (fidx - (double)idx);
  fidx = (double)(v_fdlst.size() - 1) * 0.75;
  idx = (uint)fidx;
  fditr = v_fdlst.begin();
  for (uint i = 0; i < idx; ++i) ++fditr;
  xpre = *fditr;
  ++fditr;
  v_q3 = xpre + (*fditr - xpre) * (fidx - (double)idx);
  return 1;
}

double stats::ptv(double pt)    // パーセンタイルを返す
{
  if (v_fdlst.empty()) return BLANKf; // データ無し
  if (v_min == BLANKf) this->sort(); // ソートする
  listdbl::const_iterator fditr = v_fdlst.begin();
  if (v_fdlst.size() < 2) return *fditr; // データが1つだけの時
  if (pt >= 100.0) return v_max;
  double fidx = (double)(v_fdlst.size() - 1) * pt * 0.01;
  uint idx = (uint)fidx;
  fditr = v_fdlst.begin();
  for (uint i = 0; i < idx; ++i) ++fditr;
  double xpre = *fditr;
  ++fditr;
  return xpre + (*fditr - xpre) * (fidx - (double)idx);
}

class statistclass              // statistメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout, option& opt);
  // 分析結果のタイトルとデータ識別情報を出力する
  void put_stat(ostream& fout, option& opt, stats& statsj);
  // 基礎統計量を出力する
  void put_mci(ostream& fout, option& opt, const string& itemj, stats& statsj) const;
  // 各種の信頼区間を出力する
  void put_norm(ostream& fout, option& opt, const string& itemj, stats& statsj);
  // 正規性の指標を出力する
  void cal_sk(double& g1, double& v1, string& skew, double& g2, double& v2, string& kur, stats& statsj) const;
  // 歪度と尖度を計算する
  void put_dns(ofstream& fout, option& opt, const string& itemj, stats& statsj);
  // 計算結果をデータファイル形式で出力する
  void put_dnsid(ofstream& fout);
  // DANSデータのタイトルとデータ識別情報を出力する
  void plot_ngp(ofstream& fout, option& opt, const string& itemj, stats& statsj, bool delml);
  // Ngraph用データを出力する
  void plot_ngpid(ofstream& fout, option& opt);
  // Ngraph用データのタイトルとデータ識別情報を出力する
  void plot_box(ostream& fout, option& opt, const string& itemj, stats& statsj) const;
  // Ngraph用箱ヒゲプロットデータを出力する
  void plot_drd(ostream& fout, option& opt, ngptype ngp, const string& itemj, stats& statsj) const;
  // Ngraph用度数分布／相対度数分布データを出力する
  void plot_cdrcd(ostream& fout, option& opt, ngptype ngp, const string& itemj, stats& statsj) const;
  // Ngraph用累積度数分布＆プロット／相対累積度数分布＆プロットデータを出力する
  void plot_ndrd(ostream& fout, option& opt, ngptype ngp, const string& itemj, stats& statsj) const;
  // Ngraph用理論度数分布／理論相対度数分布データを出力する
  void plot_ncdrcd(ostream& fout, option& opt, ngptype ngp, const string& itemj, stats& statsj) const;
  // Ngraph用理論累積度数分布／理論相対累積度数分布データを出力する
  void plot_vrd(ostream& fout, option& opt, ngptype ngp, const string& itemj, stats& statsj) const;
  // Ngraph用仮想的相対度数分布データを出力する
  void plot_vnrd(ostream& fout, option& opt, ngptype ngp, const string& itemj, stats& statsj) const;
  // Ngraph用仮想的理論相対度数分布データを出力する
public:
  statistclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~statistclass()               // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込む
  void put(option& opt);                     // 基礎統計量を出力する
};

void statistclass::put_resid(ostream& fout, option& opt)
// 分析結果のタイトルとデータ識別情報を出力する
{
  puttitle(fout, Title);
  fout << "\n";
  id.put_title(fout);
  id.put_stra(fout);
  if (opt.sep() == false) id.put_item(fout);
}

void statistclass::put_stat(ostream& fout, option& opt, stats& statsj)
// 基礎統計量を出力する
{
  constexpr uint Npt = 4;       // 1行あたりのパーセンタイルの個数
  int n = statsj.n();
  double m = statsj.mean();
  double cv = (sgn(m) == 0) ? 0.0 : statsj.sd() / m * 100.0;
  putmsd(fout, n, m, statsj.vr());
  fout << "最小値=" << statsj.min() << "  最大値=" << statsj.max();
  fout << "  範囲=" << statsj.max() - statsj.min() << "  変動係数=" << cv << "%\n";
  fout << "  四分位数:25%値=" << statsj.q1() << "  中央値=" << statsj.med() << "  75%値=" << statsj.q3();
  if (opt.ptsize() > 0)       // パーセンタイル指定時
  {
    fout << "\n  %tile:";
    for (size_t i = 0; i < opt.ptsize(); ++i)
    {
      if (i > 0)
      {
        if ((i % Npt) == 0) fout << "\n      ";
        fout << "  ";
      }
      fout << opt.pt(i) << "%値=" << statsj.ptv(opt.pt(i));
    }
  }
  fout << endl;
}

void statistclass::put_mci(ostream& fout, option& opt, const string& itemj, stats& statsj) const
// 各種の信頼区間を出力する
{
  int n = statsj.n();
  double m = statsj.mean();
  double t = (double)(ptot((ldbl)(1.0 - opt.cc()), n - 1));
  double ci = t * statsj.se();
  fout << "  平均値の" << opt.cc() * 100.0 << "%信頼区間:" << m - ci << "〜" << m + ci;
  fout << " (" <<  m << "±" << ci << ")\n";
  ci = t * sqrt((1.0 + 1.0 / (double)n) * statsj.vr());
  fout << "  データの" << opt.cc() * 100.0 << "%予測区間:" << m - ci << "〜" << m + ci;
  fout << " (" <<  m << "±" << ci << ")\n";
  t = (double)(ptoz((ldbl(1.0 - opt.cc()))));
  uint r = (uint)((double)(n + 1) * 0.5 - t * sqrt((double)n) * 0.5);
  uint i = 1;
  listdbl::const_iterator fditr = (statsj.fdlst()).begin();
  while (i < r && fditr != (statsj.fdlst()).end()) { ++i;  ++fditr; }
  double medlow = (fditr != (statsj.fdlst()).end()) ? *fditr : statsj.min();
  while (i < n - r + 1 && fditr != (statsj.fdlst()).end()) { ++i;  ++fditr; }
  double medup = (fditr != (statsj.fdlst()).end()) ? *fditr : statsj.max();
  fout << "  中央値の" << opt.cc() * 100.0 << "%信頼区間:" << medlow << "〜" << medup;
  fout << " (" << r << "番目〜" << n - r + 1 << "番目のデータ)\n";
  double ah = (1.0 - opt.cc()) * 50.0;
  medlow = statsj.ptv(ah);
  medup = statsj.ptv(100.0 - ah);
  fout << "  データの" << ah << "%点〜" << 100.0 - ah << "%点:" << medlow << "〜" << medup << endl;
}

void statistclass::put_norm(ostream& fout, option& opt, const string& itemj, stats& statsj)
// 正規性の指標を出力する
{
  double g1 = 0.0, v1 = 0.0, g2 = 0.0, v2 = 0.0;
  string skew, kur;
  cal_sk(g1, v1, skew, g2, v2, kur, statsj);
  double z = (v1 > 0.0) ? g1 / sqrt(v1) : 0.0;
  double zc = 0.0, ci = 0.0;
  ldbl pr = ztop(z);
  fout << "歪度(skewness):g1=" << left << setw(10) << g1;
  fout << "正規検定:z=" << setw(10) << z << "p=" << setw(12) << pr;
  fout << setw(4) << pind(pr) << right << skew << "\n";
  z = (v2 > 0.0) ? g2 / sqrt(v2) : 0.0;
  pr = ztop((ldbl)z);
  fout << "尖度(kurtosis):g2=" << left << setw(10) << g2;
  fout << "正規検定:z=" << setw(10) << z << "p=" << setw(12) << pr;
  fout << setw(4) << pind(pr) << right << kur << "\n";
  double jb = (double)(statsj.n()) * (g1 * g1 + g2 * g2 / 4.0) / 6.0;
  pr = ctop(jb, 2);
  fout << "Jarque-Bera検定:JB=χ2乗=" <<  jb << "  自由度=2  p=" << pr << pind(pr) << "\n";
  zc = (double)(ptoz((ldbl)(1.0 - opt.cc())));
  ci = zc * sqrt(v1);
  fout << "  g1の" << opt.cc() * 100.0 << "%信頼区間:" << g1 - ci << "〜" << g1 + ci;
  fout << " (" <<  g1 << "±" << ci << ")\n";
  ci = zc * sqrt(v2);
  fout << "  g2の" << opt.cc() * 100.0 << "%信頼区間:" << g2 - ci << "〜" << g2 + ci;
  fout << " (" <<  g2 << "±" << ci << ")\n";
}

void statistclass::cal_sk(double& g1, double& v1, string& skew, double& g2, double& v2, string& kur, stats& statsj) const
// 歪度と尖度を計算する
{
  g1 = v1 = g2 = v2 = 0.0;  skew.erase();  kur.erase();
  int n = statsj.n();
  if (n < 4) return;
  double k2 = statsj.k2(), k3 = statsj.k3(), k4 = statsj.k4();
  if (k2 > 0.0) { g1 = k3 / k2 / sqrt(k2);  g2 = k4 / k2 / k2; }
  v1 = 6.0 * (double)n * (double)(n - 1) / (double)(n - 2) / (double)(n + 1) / (double)(n + 3);
  v2 = 24.0 * (double)n * (double)(n - 1) * (double)(n - 1) / (double)(n - 3) / (double)(n - 2) / (double)(n + 3) / (double)(n + 5);
  if (g1 < -1.0) skew = "右傾斜";
  else if (fabs(g1) <= 1.0) skew = "ほぼ対称";
  else skew = "左傾斜";
  if (g2 < -3.0) kur = "鈍峰";
  else if (fabs(g2) <= 3.0) kur = "ほぼ正規";
  else kur = "鋭峰";
}

void statistclass::put_dns(ofstream& fout, option& opt, const string& itemj, stats& statsj)
// 計算結果をデータファイル形式で出力する
{
  datac cdata(18);
  if (opt.sep() || id.ni() == 1) cdata.assign(0, itemj);
  else cdata.assign(0, "全項目");
  cdata.assign(1, (double)(statsj.n()));
  cdata.assign(2, statsj.mean());
  cdata.assign(3, statsj.sd());
  cdata.assign(4, statsj.se());
  cdata.assign(5, statsj.min());
  cdata.assign(6, statsj.q1());
  cdata.assign(7, statsj.med());
  cdata.assign(8, statsj.q3());
  cdata.assign(9, statsj.max());
  cdata.assign(10, statsj.max() - statsj.min());
  double cv = (statsj.mean() == 0.0) ? 0.0 : statsj.sd() / statsj.mean() * 100.0;
  cdata.assign(11, cv);
  double g1 = 0.0, v1 = 0.0, g2 = 0.0, v2 = 0.0;
  string skew, kur;
  cal_sk(g1, v1, skew, g2, v2, kur, statsj);
  cdata.assign(12, g1);
  double z = (v1 > 0.0) ? g1 / sqrt(v1) : 0.0L;
  cdata.assign(13, z);
  cdata.assign(14, (double)(ztop(z)));
  cdata.assign(15, g2);
  z = (v2 > 0.0) ? g2 / sqrt(v2) : 0.0;
  cdata.assign(16, z);
  cdata.assign(17, (double)(ztop(z)));
  cdata.put(fout);
}

void statistclass::put_dnsid(ofstream& fout)
// DANSデータのタイトルとデータ識別情報を出力する
{
  idd tid(id.title(), 18);
  id.put_title(fout, formtag);
  id.put_stra(fout, formtag);
  tid.assign_item(0, "項目名");
  tid.assign_item(1, "例数");
  tid.assign_item(2, "平均値");
  tid.assign_item(3, "標準偏差");
  tid.assign_item(4, "標準誤差");
  tid.assign_item(5, "最小値");
  tid.assign_item(6, "25%値");
  tid.assign_item(7, "中央値");
  tid.assign_item(8, "75%値");
  tid.assign_item(9, "最大値");
  tid.assign_item(10, "範囲");
  tid.assign_item(11, "変動係数(%)");
  tid.assign_item(12, "歪度");
  tid.assign_item(13, "歪度のz値");
  tid.assign_item(14, "歪度のp値");
  tid.assign_item(15, "尖度");
  tid.assign_item(16, "尖度のz値");
  tid.assign_item(17, "尖度のp値");
  tid.put_item(fout, formtag);
}

void statistclass::plot_ngp(ofstream& fout, option& opt, const string& itemj, stats& statsj, bool delml)
// Ngraph用データを出力する
{
  if (opt.ngp("box"))           // 箱ひげプロット
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -drectangle -x1 -y3 -w80\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_box(fout, opt, itemj, statsj);
  }
  if (opt.ngp("cd"))            // 累積度数分布
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline -w80\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_cdrcd(fout, opt, ngpcd, itemj, statsj);
  }
  if (opt.ngp("cdp"))           // 累積度数分布プロット
  {
    if (delml)
    { fout << NGPdelm << "\n";}
    else
    {
      fout << "-s1 -dmark,00\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_cdrcd(fout, opt, ngpcdp, itemj, statsj);
  }
  if (opt.ngp("d"))             // 度数分布
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -drectangle -x1 -y3 -w80\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_drd(fout, opt, ngpd, itemj, statsj);
  }
  if (opt.ngp("ncd"))           // 理論累積度数分布
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline -w80\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_ncdrcd(fout, opt, ngpncd, itemj, statsj);
  }
  if (opt.ngp("nd"))            // 理論度数分布
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline -w80\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_ndrd(fout, opt, ngpnd, itemj, statsj);
  }
  if (opt.ngp("nrcd"))          // 理論相対累積度数分布
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline -w80 -zy0,1,0.2\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_ncdrcd(fout, opt, ngpnrcd, itemj, statsj);
  }
  if (opt.ngp("nrd"))           // 理論相対度数分布
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline -w80\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_ndrd(fout, opt, ngpnrd, itemj, statsj);
  }
  if (opt.ngp("rcd"))           // 相対累積度数分布
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline -w80 -zy0,1,0.2\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_cdrcd(fout, opt, ngprcd, itemj, statsj);
  }
  if (opt.ngp("rcdp"))          // 相対累積度数分布プロット
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dmark,00\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_cdrcd(fout, opt, ngprcdp, itemj, statsj);
  }
  if (opt.ngp("rd"))            // 相対度数分布
  {
    if (delml) fout << NGPdelm << "\n";
    else
    {
      fout << "-s1 -drectangle -x1 -y3 -w80\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_drd(fout, opt, ngprd, itemj, statsj);
  }
  if (opt.ngp("vrd"))           // 仮想的相対度数分布
  {
    if (delml) fout << NGPdelm << "\n";
    else
    {
      fout << "-s1 -dmark,00\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_vrd(fout, opt, ngprd, itemj, statsj);
  }
  if (opt.ngp("vnrd"))          // 仮想的理論相対度数分布
  {
    if (delml) fout << NGPdelm << "\n";
    else
    {
      fout << "-s1 -dline -w80 -zy0,1,0.2\n";
      delml = true;
    }
    plot_ngpid(fout, opt);
    plot_vnrd(fout, opt, ngprd, itemj, statsj);
  }
}

void statistclass::plot_ngpid(ofstream& fout, option& opt)
// Ngraph用データのタイトルとデータ識別情報を出力する
{
  fout << "%";  puttitle(fout, Title);
  fout << "%\n";
  id.put_title(fout, formngp);
  id.put_stra(fout, formngp);
  if (opt.sep() == false) id.put_item(fout, formngp);
}

void statistclass::plot_box(ostream& fout, option& opt, const string& itemj, stats& statsj) const
// Ngraph用箱ヒゲプロットデータを出力する
{
  if (opt.sep()) fout << "%集計項目:" << itemj << "\n";
  fout << "%箱ヒゲプロット:箱\(-s1 zx0.5,1.5,0.5 -drectangle -x1 -y3 -w80)\n";
  double x = 1.0;
  fout << x - 0.1 << " " << statsj.q1() << " " << x + 0.1 << " " << statsj.med() << "\n";
  fout << x - 0.1 << " " << statsj.med() << " " << x + 0.1 << " " << statsj.q3() << "\n";
  double hw = (statsj.q3() - statsj.q1()) * 1.5;
  double hlow = statsj.q1() - hw, hup = statsj.q3() + hw;
  double low = statsj.min(), up = statsj.max();
  bool shclow = true, shcup = true;
  listdbl::const_iterator fditr = (statsj.fdlst()).begin();
  while ((shclow || shcup) && fditr != (statsj.fdlst()).end())
  {
    if (shclow && hlow <= *(fditr)) { shclow = false;  low = *(fditr); }
    if (shcup && hup < *(fditr)) shcup = false;
    else up = *(fditr);
    ++fditr;
  }
  fout << NGPdelm2 << "\n";
  fout << "%箱ヒゲプロット:ヒゲ(-zx0.5,1.5,0.5 -derrorbar_y -x1 -y2 -w40 -o500)\n";
  fout << x << " " << statsj.q1() << " " << low - statsj.q1() << " 0\n";
  fout << x << " " << statsj.q3() << " 0 " << up - statsj.q3() << "\n";
  fout << NGPdelm2 << "\n";
  fout << "%箱ヒゲプロット:平均値(-s1 -zx0.5,1.5,0.5 -dmark,72 -w80 -o400)\n";
  fout << x << " " << statsj.mean() << "\n";
  fout << NGPdelm2 << "\n";
  fout << "%箱ヒゲプロット:外れ値(-s1 -zx0.5,1.5,0.5 -dmark,2 -w80 -o400)\n";
  fditr = (statsj.fdlst()).begin();
  while (fditr != (statsj.fdlst()).end())
  {
    if (*(fditr) < low) fout << x << " " << *(fditr) << endl;
    else if (up < *(fditr)) fout << x << " " << *(fditr) << endl;
    ++fditr;
  }
}

void statistclass::plot_drd(ostream& fout, option& opt, ngptype ngp, const string& itemj, stats& statsj) const
// Ngraph用度数分布／相対度数分布データを出力する
{
  if (opt.sep())
  {
    if (ngp == ngpd) fout << "%度数分布";
    else fout << "%相対度数分布";
    fout << " 集計項目:" << itemj << "\n";
  }
  double cmin = opt.cmin();
  if (cmin == BLANKf) cmin = statsj.min();
  double h = opt.h();
  if (h == BLANKf) h = statsj.sd() / 3.0; // Charlierの級間基準
  double xlow = cmin;
  int i = 0;
  listdbl::const_iterator fditr = (statsj.fdlst()).begin();
  while (fditr != (statsj.fdlst()).end())
  {
    double xup = xlow + h;
    int ilow = i;
    while (fditr != (statsj.fdlst()).end() && *fditr < xup)
    { ++i;  ++fditr; }
    int f = i - ilow;
    if (f > 0)
    {
      double y = (ngp == ngpd) ? (double)f : (double)f / (double)(statsj.n());
      fout << xlow << " 0 " << xup << " " << y << "\n";
    }
    xlow = xup;
  }
}

void statistclass::plot_cdrcd(ostream& fout, option& opt, ngptype ngp, const string& itemj, stats& statsj) const
// Ngraph用累積度数分布＆プロット／相対累積度数分布＆プロットデータを出力する
{
  if (opt.sep())
  {
    if (ngp == ngpcd) fout << "%累積度数分布";
    else if (ngp == ngpcdp) fout << "%累積度数分布のプロット";
    else if (ngp == ngprcd) fout << "%相対累積度数分布";
    else fout << "%相対累積度数分布のプロット";
    fout << " 集計項目:" << itemj << "\n";
  }
  double range = statsj.max() - statsj.min();
  double xmin = statsj.min() - range * 0.2, xmax = statsj.max() + range * 0.2;
  int n = statsj.n(), ni = 0;
  int di = (n > (int)RDdef) ? n / (int)RDdef : 0;
  double ypre = 0.0;
  if (ngp == ngpncd || ngp == ngpnrcd) fout << xmin << " 0\n";
  for (double fdata : statsj.fdlst())
  {
    double x = fdata;
    ++ni;
    double y = (double)ni;
    if (di == 0 || ni % di == 0) // プロット出力
    {
      if (ngp == ngpcd || ngp == ngprcd) fout << x << " " << ypre << "\n";
      if (ngp == ngprcd || ngp == ngprcdp) y /= (double)n;
      fout << x << " " << y << "\n";
      ypre = y;
    }
  }
  if (ngp == ngpncd || ngp == ngpnrcd) fout << xmax << " " << ypre << endl;
}

void statistclass::plot_ndrd(ostream& fout, option& opt, ngptype ngp, const string& itemj, stats& statsj) const
// Ngraph用理論度数分布／理論相対度数分布データを出力する
{
  if (opt.sep())
  {
    if (ngp == ngpnd) fout << "%理論度数分布";
    else fout << "%理論相対度数分布";
    fout << " 集計項目:" << itemj << "\n";
  }
  double sd = statsj.sd();
  double h = opt.h();
  if (h == BLANKf) h = sd / 3.0; // Charlierの級間基準
  double range = sd * 8.0;
  double xmin = statsj.mean() - sd * 4.0, xmax = statsj.mean() + sd * 4.0;
  double dx = (xmax - xmin) / (double)RDdef, x = xmin;
  fout << xmin << " 0\n";
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    ldbl pr = sndis((ldbl)((x + h * 0.5 - statsj.mean()) / sd)) - sndis((ldbl)((x - h * 0.5 - statsj.mean()) / sd));
    ldbl y = (ngp == ngpnd) ? pr * (ldbl)(statsj.n()) : pr;
    fout << x << " " << y << "\n";
  }
}

void statistclass::plot_ncdrcd(ostream& fout, option& opt, ngptype ngp, const string& itemj, stats& statsj) const
// Ngraph用理論累積度数分布／理論相対累積度数分布データを出力する
{
  if (opt.sep())
  {
    if (ngp == ngpncd) fout << "%理論累積度数分布";
    else fout << "%理論相対累積度数分布";
    fout << " 集計項目:" << itemj << "\n";
  }
  double mean = statsj.mean(), sd = statsj.sd(), ndbl = (double)statsj.n();
  double range = sd * 8.0;
  double xmin = mean - sd * 4.0, xmax = mean + sd * 4.0;
  double dx = range / (double)RDdef, x = xmin;
  ldbl y = 0.0L;
  fout << xmin << " 0\n";
  for (uint i = 0; i < RDdef; ++i, x += dx)
  {
    ldbl pr = sndis((x - mean) / sd);
    y = (ngp == ngpncd) ? pr * ndbl : pr;
    fout << x << " " << y << "\n";
  }
  y = (ngp == ngpncd) ? ndbl : 1.0L;
  fout << xmax << " " << y << endl;
}

void statistclass::plot_vrd(ostream& fout, option& opt, ngptype ngp, const string& itemj, stats& statsj) const
// Ngraph用仮想的相対度数分布データを出力する
{
  if (opt.sep()) fout << "%仮想的相対度数分布 集計項目:" << itemj << "\n";
  double range = statsj.max() - statsj.min();
  double xmin = statsj.min() - range * 0.2, xmax = statsj.max() + range * 0.2;
  int n = statsj.n(), ni = 0, n0 = 0, npre = 0;
  int di = (n > (int)RDdef) ? n / (int)RDdef : 0;
  double mean = statsj.mean(), sd = statsj.sd();
  double dxmin = range / (double)n; // dxの最小値
  double cf = dxmin / (double)n;
  double x0 = xmin - dxmin, xpre = xmin - dxmin;
  for (double fdata : statsj.fdlst())
  {
    ++ni;
    double xi = fdata;
    if (di == 0 || ni % di == 0) // プロット出力
    {
      double dxi = xi - xpre, dyi = (double)(ni - npre);
      if (dxi >= dxmin)
      { n0 = npre;  x0 = xpre; }
      else
      {
        dxi = xi - x0;  dyi = (double)(ni - n0);
      }
      double yi = dyi / dxi * cf;
      fout << xi << " " << yi << "\n";
      npre = ni;  xpre = xi;
    }
  }
}

void statistclass::plot_vnrd(ostream& fout, option& opt, ngptype ngp, const string& itemj, stats& statsj) const
// Ngraph用仮想的理論相対度数分布データを出力する
{
  if (opt.sep()) fout << "%仮想的理論相対度数分布 集計項目:" << itemj << "\n";
  double mean = statsj.mean(), sd = statsj.sd();
  double range = sd * 8.0;
  double xmin = mean - sd * 4.0, xmax = mean + sd * 4.0;
  double dx = range / (double)RDdef, x = xmin;
  double h = (statsj.max() - statsj.min()) / (double)statsj.n();
  for (uint i = 0; i <= RDdef; ++i, x += dx)
  {
    ldbl y = sndis((ldbl)((x - mean) / sd)) - sndis((ldbl)((x - h - mean) / sd));
    fout << x << " " << y << "\n";
  }
}

int statistclass::get(const char* c_fname, option& opt)
// データを読み込む
{
  if (c_fname == nullptr) fname.erase();
  else fname = c_fname;
  int nc = fgetdata(fname, id, datalst);
  if (nc < 0) throw runtime_error("実行時エラー！");
  if (opt.check() != 0) throw option_error("実行時オプションエラー！");
  return nc;
}

void statistclass::put(option& opt) // 基礎統計量を出力する
{
  constexpr int LENline = 76;   // 区切り線の長さ
  if (opt.res("res")) put_resid(cout, opt); // 分析結果出力
  ofstream foutdns;
  if (opt.dnsname() != "")      // DANSデータ出力
  {
    if (ofopen(foutdns, (opt.dnsname()).c_str()) == false) return; // ファイルオープン失敗
    put_dnsid(foutdns);
  }
  ofstream foutngp;
  if (opt.ngpname() != "")      // Ngraphデータ出力
  {
    if (ofopen(foutngp, (opt.ngpname()).c_str()) == false) return; // ファイルオープン失敗
  }
  stats statsj;                 // 基礎統計量+k統計量+四分位値
  size_t jstart = 0, jend = id.ni();
  if (opt.sep() == false) { jstart = 0;  jend = 1; } // 全項目を一度に集計
  for (size_t j = jstart; j < jend; ++j)
  {
    statsj.clear();
    if (opt.sep())
    {
      for (auto& cdatai : datalst) statsj.sigma(cdatai.fdata(j));
    }
    else
    {
      for (auto& cdatai : datalst)
      {
        for (size_t i = 0; i < id.ni(); ++i) statsj.sigma(cdatai.fdata(i));
      }
    }
    if (statsj.n() > 0)
    {
      statsj.sort();
      if (opt.res("res"))       // 計算結果出力
      {
        if (opt.sep()) cout << "集計項目:" << id.item(j) << "\n";
        bool putres = false;
        if (opt.res("stat"))
        {
          putline(cout, LENline);
          put_stat(cout, opt, statsj);
          put_mci(cout, opt, id.item(j), statsj);
          putres = true;
        }
        if (opt.res("norm"))
        {
          putline(cout, LENline);
          put_norm(cout, opt, id.item(j), statsj);
          putres = true;
        }
        if (putres) putline(cout, LENline);
        cout << "\n";
      }
      if (opt.dnsname() != "") put_dns(foutdns, opt, id.item(j), statsj);
      if (opt.ngpname() != "") plot_ngp(foutngp, opt, id.item(j), statsj, (j > jstart) ? true : false);
    }
  }
  if (opt.res("res")) cout << endl;
}

int main(int argc, char* argv[]) // メインルーチン
{
  statistclass statist;         // メインクラス変数
  option opt(argv[0]);          // オプション
  try
  {
    opt.setdef();               // オプションの初期設定
    bool fok = true;            // ファイル処理フラッグ
    for (int i = 1; i <= argc; ++i)
    {
      if (i >= argc)
      {
        if (fok) statist.get(nullptr, opt); // 標準入力から読み込む
      }
      else if (*(argv[i]) == CHRopt)
      { opt.set(argv[i] + 1);  fok = true; }
      else
      { statist.get(argv[i], opt);  fok = false; }
    }
    statist.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;
}
