// sstat.C Version 7.45  Copyleft (c) 杉本典夫
// Time-stamp: <2026-01-03 16:46:18 tac>
//
// 実数をデータを層別集計し、基礎統計量を出力する。複数項目を一度に集計可能。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -w sstat sstat.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 =            // プログラム使用法表示用
{
  "sstat  Version 7.45          Copyleft(Vo\\oV) 2026 杉本典夫",
  "書式:sstat  [<オプション> …]  [<データファイル> …]",
  "機能:各項目の基礎統計量を層別集計して出力する。",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定する。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -dnsname[=<出力ファイル名>]:DANSデータを出力するファイル名を指定する。",
  "    オプション省略時:DANSデータ出力せず  -dnsname:-dnsname=sstat_dns.d",
  "    出力項目:基礎統計量",
  "  -f=<ファイル名>:指定されたファイルからオプションを読み込む。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示する。(-h:-h=man)",
  "  -ol[={-|md|mn|sd|<基準値>}[<係数>]]:外れ値の基準値を指定する。(-ol:-ol=sd2)",
  "    -:外れ値を出力しない(オプション省略時)",
  "  -pt=<％tile>[,…]:パーセンタイルを出力する。(0＜％tile≦100)",
  "  -res=[-][{def|all|norm|si|st|stat}][,…]",
  "      :解析結果の出力項目を指定する。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    norm:正規性の指標  si:層ごとの集計結果を出力する",
  "    st:全体の集計結果を出力する  stat:基礎統計量",
  "    オプション省略時:全て出力"
};

const string Title = "=== 基礎統計量 ==="; // 表題
const string DNSNAMEdef = "sstat_dns.d"; // 省略時DANSファイル名

enum oltype                     // 外れ値の基準値の種類
{
  olnon,                        // 出力なし
  olmn,                         // 平均値
  olmd,                         // 中央値
  olsd,                         // SD
  oletc                         // その他
};

class option                    // オプション集合体
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  string v_dnsname;             // DANSファイル名
  oltype v_ol;                  // 外れ値の基準値
  double v_olc;                 // 外れ値の係数
  lstr ptlst;                   // パーセンタイルリスト
  vecdbl ptvec;                 // パーセンタイルベクトル
  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; }
  const string& dnsname() const // DANSデータファイル名を返す
  { return v_dnsname; }
  oltype ol() const             // 外れ値の基準値を返す
  { return v_ol; }
  double olc() const            // 外れ値の係数を返す
  { return v_olc; }
  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); // 分析結果出力フラッグを返す
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check(size_t ni);         // オプションのセットとチェック
  int set_pt();                 // -ptオプションのセット
  int set_res();                // -resオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_cc(CCdef), v_dnsname(), v_ol(olnon), v_olc(0.0), ptlst(), ptvec(), 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::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] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else if (opt[0] == "ol")      // 外れ値の基準値と係数の指定
  {
    if (opt[1].empty())
    { v_ol = olsd; }
    else if (opt[1] == "-")
    { v_ol = olnon; }
    else if (opt[1][0] == 'm')
    {
      if (opt[1][1] == 'd') v_ol = olmd;
      else if (opt[1][1] == 'n') v_ol = olmn;
      else err = badopt(pname(), optstr, USAGE);
    }
    else if (opt[1][0] == 's' && opt[1][1] == 'd')
    { v_ol = olsd; }
    else if (isdigit(opt[1][0]) || (opt[1][0] == '-' && isdigit(opt[1][1])))
    {
      v_ol = oletc;  v_olc = strtodf(opt[1]);
      if (v_olc == BLANKf) err = badopt(pname(), optstr, USAGE);
    }
    else
    { err = badopt(pname(), optstr, USAGE); }
    if (err == 0 && v_ol != oletc)
    {
      v_olc = (opt[1].size() < 3) ? 2.0 : strtodf(opt[1].substr(2));
      if (v_olc < 0.0 || 1000.0 < v_olc) 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                          // 間違ったオプション指定
  { err = nonopt(pname(), optstr, USAGE); }
  if (err > 0) throw help_except("ヘルプ表示");
  if (err < 0) throw option_error(string("オプション指定違反→") + optstr);
}

int option::check(size_t ni)    // オプションのセットとチェック
{
  if (ni < 2) return disperr(pname(), "項目が2つ以上必要です！");
  if (set_pt() < 0) return -2;  // -ptオプションのセット
  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_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["norm"] = resmap["si"] = resmap["st"] = resmap["stat"] = true;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["norm"] = resmap["si"] = resmap["st"] = 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["si"] = resmap["st"] = resmap["stat"] = true;
      }
      else if (respos == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

class stats : public statk<double> // 基礎統計量+k統計量+四分位値クラス
{
private:
  double v_min;                 // 最小値
  double v_q1;                  // 第1四分位値(25パーセンタイル)
  double v_med;                 // 中央値(第2四分位値)
  double v_q3;                  // 第3四分位値(75パーセンタイル)
  double v_max;                 // 最大値
  int v_nol;                    // 外れ値の例数
  list<double> v_fdlst;         // 実数データリスト
public:
  stats() : statk(), v_min(BLANKf), v_q1(BLANKf), v_med(BLANKf), v_q3(BLANKf), v_max(BLANKf), v_nol(0), v_fdlst()
  // デフォールトコンストラクタ
  {}
  stats(int n0, double sig0, double sig20, double sigmax30, double sigmax40, double min0, double q10, double med0, double q30, double max0, double nol0, list<double>& fdls0) : statk(n0, sig0, sig20, sigmax30, sigmax40), v_min(min0), v_q1(q10), v_med(med0), v_q3(q30), v_max(max0), v_nol(nol0), v_fdlst(fdls0)
  // コンストラクタ
  {}
  stats(const stats& stats0) : statk(stats0.vn, stats0.sigmax, stats0.sigmax2, stats0.sigmax3, stats0.sigmax4), 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_nol(stats0.v_nol), v_fdlst(stats0.v_fdlst)
  // コピーコンストラクタ
  {}
  virtual ~stats()              // デストラクタ
  {}
  const stats& operator=(const stats& stats0); // コピー代入演算子
  void clear();                                // データをクリアする
  int sigma(double x);                         // データの累積集計
  int sort(oltype ol, double olc);             // データをソートする
  const list<double>& fdlst() const // 実数データリストを返す
  { return v_fdlst; }
  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; }
  int nol() const               // 外れ値の例数を返す
  { return v_nol; }
  double ptv(double pt);        // パーセンタイル値を返す
};

const stats& stats::operator=(const stats& stats0)
// コピー代入演算子
{
  if (&stats0 == this) return *this; // 自分自身への代入を防ぐ
  vn = stats0.vn;  sigmax = stats0.sigmax;  sigmax2 = stats0.sigmax2;
  sigmax3 = stats0.sigmax3;  sigmax4 = stats0.sigmax4;
  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_nol = stats0.v_nol;
  v_fdlst = stats0.v_fdlst;
  return *this;
}

void stats::clear()             // データをクリアする
{
  vn = 0;  sigmax = sigmax2 = sigmax3 = sigmax4 = 0.0;
  v_min = v_q1 = v_med = v_q3 = v_max = BLANKf;  v_nol = 0;
  v_fdlst.clear();
}

int stats::sigma(double x)      // データの累積集計
{
  if (x == BLANKf) return 0;    // 累積せず
  double x2 = x * x;
  ++vn;  sigmax += x;  sigmax2 += x2;
  sigmax3 += x * x2;  sigmax4 += x2 * x2;
  v_fdlst.push_back(x);
  return 1;                     // 累積
}

int stats::sort(oltype ol, double olc) // データをソートする
{
  if (v_fdlst.empty()) return 0; // データ無し
  if (v_min != BLANKf) return 1; // ソート済み
  v_fdlst.sort();
  list<double>::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;
  int idx = (int)fidx;
  fditr = v_fdlst.begin();
  for (int 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 = (int)fidx;
  fditr = v_fdlst.begin();
  for (int 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 = (int)fidx;
  fditr = v_fdlst.begin();
  for (int i = 0; i < idx; ++i) ++fditr;
  xpre = *fditr;
  ++fditr;
  v_q3 = xpre + (*fditr - xpre) * (fidx - (double)idx);
  double lim = 0.0;
  switch (ol)
  {
  case olmn : lim = olc * mean();
    break;
  case olmd : lim = olc * v_med;
    break;
  case olsd : lim = mean() + olc * sd();
    break;
  default : lim = olc;
  }
  idx = 0;
  for (fditr = v_fdlst.begin(); *fditr <= lim && fditr != v_fdlst.end(); ++fditr) ++idx;
  v_nol = this->n() - idx;       // 外れ値の例数
  return 1;
}

double stats::ptv(double pt)    // パーセンタイル値を返す
{
  if (v_fdlst.empty()) return BLANKf; // データ無し
  if (v_min == BLANKf) sort(olnon, 0.0); // ソートする
  list<double>::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;
  int idx = (int)fidx;
  fditr = v_fdlst.begin();
  for (int i = 0; i < idx; ++i) ++fditr;
  double xpre = *fditr;
  ++fditr;
  return xpre + (*fditr - xpre) * (fidx - (double)idx);
}

using vecstats = vector<stats>; // 基礎統計量ベクトル配列の別名定義
using mapsvecs = map<string, vecstats>; // 連想配列用マップの短名定義
using mapsvecs_itr = map<string, vecstats>::iterator; // 反復子の短名定義
using mapsvecs_citr = map<string, vecstats>::const_iterator; // const反復子の短名定義

class varclass                  // 変数クラス
{
private:
  mapsvecs v_ssvecmap;          // q基礎統計量ベクトルマップ
  vecstats v_ssvectot;          // 全体の基礎統計量ベクトル
public:
  varclass(size_t ni) : v_ssvecmap(), v_ssvectot(ni - 1)
  // デフォールトコンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  size_t nmap() const           // マップ数を返す
  { return v_ssvecmap.size(); }
  mapsvecs& ssvecmap()          // 基礎統計量ベクトルマップを返す
  { return v_ssvecmap; }
  mapsvecs_itr ssvecpos(const string& code) // マップ反復子を返す
  { return v_ssvecmap.find(code); }
  void cal(option& opt, ldata& datalst); // データを集計する
};

void varclass::cal(option& opt, ldata& datalst)
// データを集計する
{
  for (auto& cdatai : datalst)
  {
    if (cdatai.sdata(0) != BLANKs) // 層別コードあり
    {
      mapsvecs_itr spos = v_ssvecmap.find(cdatai.sdata(0));
      if (spos != v_ssvecmap.end()) // 既存コード
      {
        for (size_t j = 1; j < cdatai.size(); ++j)
        {
          double fdata = cdatai.fdata(j);
          if (fdata != BLANKf)
          { spos->second[j - 1].sigma(fdata);  v_ssvectot[j - 1].sigma(fdata); }
        }
      }
      else                      // 新規コード
      {
        vecstats tssvec(cdatai.size() - 1);
        for (size_t j = 1; j < cdatai.size(); ++j)
        {
          double fdata = cdatai.fdata(j);
          if (fdata != BLANKf)
          { tssvec[j - 1].sigma(fdata);  v_ssvectot[j - 1].sigma(fdata); }
        }
        v_ssvecmap.insert(mapsvecs::value_type(cdatai.sdata(0), tssvec));
      }
    }
  }
  v_ssvecmap.insert(mapsvecs::value_type("全体", v_ssvectot));
}

class sstatclass                // sstatメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout, const string& code);
  // 分析結果のタイトルとデータ識別情報を出力する
  void put_tab(ostream& fout, option& opt, vecstats& statsvec);
  // 集計結果を集計表形式で出力する
  void put_stat(ostream& fout, option& opt, stats& statsj);
  // 基礎統計量を出力する
  void put_mci(ostream& fout, option& opt, stats& statsj) const;
  // 各種の信頼区間を出力する
  void put_norm(ostream& fout, option& opt, stats& statsj);
  // 正規性の指標を出力する
  void cal_sk(double& g1, double& v1, string& skew, double& g2, double& v2, string& kur, stats& statj);
  // 歪度と尖度を計算する
  void put_dnsid(ofstream& fout, option& opt);
  // DANSデータのタイトルとデータ識別情報を出力する
  void put_dns(ofstream& fout, option& opt, mapsvecs_itr spos);
  // 集計結果をデータファイル形式で出力する
public:
  sstatclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~sstatclass()                 // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込む
  void put(option& opt);                     // 基礎統計量を出力する
};

void sstatclass::put_resid(ostream& fout, const string& code)
// 分析結果のタイトルとデータ識別情報を出力する
{
  puttitle(fout, Title);
  fout << "\n";
  idd tid(id);
  stra tstra;
  tstra.get('s', 1, code);
  lstra tstrals;
  tstrals.push_back(tstra);
  tid.append_stra(tstrals);     // 層別コードを追加する
  tid.put_title(fout);
  tid.put_stra(fout);
}

void sstatclass::put_tab(ostream& fout, option& opt, vecstats& statsvec)
// 集計結果を集計表形式で出力する
{
  constexpr int LENline = 76;   // 区切り線の長さ
  int nss = (int)(id.ni()) - 1;
  for (int j = 0; j < nss; ++j)
  {
    if (statsvec[j].n() > 0)
    {
      statsvec[j].sort(opt.ol(), opt.olc());
      fout << "集計項目:" << id.item(j + 1) << "\n";
      int n = statsvec[j].n();
      double m = statsvec[j].mean(), sd = statsvec[j].sd();
      bool putres = false;
      if (opt.res("stat"))
      {
        putline(fout, LENline);
        put_stat(fout, opt, statsvec[j]);
        putres = true;
      }
      if (opt.res("norm"))
      {
        putline(fout, LENline);
        put_norm(fout, opt, statsvec[j]);
        putres = true;
      }
      if (putres) putline(fout, LENline);
      fout << "\n";
    }
  }
  fout << endl;
}

void sstatclass::put_stat(ostream& fout, option& opt, stats& statsj)
// 基礎統計量を出力する
{
  constexpr uint Npt = 4;       // 1行あたりのパーセンタイルの個数
  int n = statsj.n();
  double m = statsj.mean(), sd = statsj.sd();
  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 << "\n";
  put_mci(fout, opt, statsj);
/*  double t = ptot(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 = ptoz(1.0 - opt.cc());
  int r = (int)((double)(n + 1) * 0.5 - t * sqrt((double)n) * 0.5);
  int i = 1;
  list<double>::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 << "\n"; */
}

void sstatclass::put_mci(ostream& fout, option& opt, 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;
  list<double>::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 sstatclass::put_norm(ostream& fout, option& opt, stats& statsj)
// 正規性の指標を出力する
{
  int n = statsj.n();
  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;
  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(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)(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";
  double zc = ptoz(1.0 - opt.cc());
  double 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 sstatclass::cal_sk(double& g1, double& v1, string& skew, double& g2, double& v2, string& kur, stats& statj)
// 歪度と尖度を計算する
{
  g1 = v1 = g2 = v2 = 0.0;  skew.erase();  kur.erase();
  int n = statj.n();
  if (n < 4) return;
  double k2 = statj.k2(), k3 = statj.k3(), k4 = statj.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 sstatclass::put_dnsid(ofstream& fout, option& opt)
// データ識別情報をデータファイル形式で出力する
{
  int nitem = 14;               // 出力項目数
  if (opt.ol() != olnon) nitem += 1;
  id.put_title(fout, formtag);
  id.put_stra(fout, formtag);
  idd tid(id.title(), 1 + nitem);
  tid.assign_item(0, id.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, "尖度");
  tid.assign_item(14, "合計");
  if (opt.ol() != olnon)
  {
    stringstream ostr;
    switch (opt.ol())
    {
    case olmn : ostr << opt.olc() << "x平均値";
      break;
    case olmd : ostr << opt.olc() << "x中央値";
      break;
    case olsd : ostr << "平均値+" << opt.olc() << "SD";
      break;
    default : ostr << opt.olc();
      break;
    }
    tid.assign_item(15, ostr.str()+"を越える割合");
  }
  tid.put_item(fout, formtag);
}

void sstatclass::put_dns(ofstream& fout, option& opt, mapsvecs_itr spos)
// 集計結果をデータファイル形式で出力する
{
  int nitem = 14;               // 出力項目数
  if (opt.ol() != olnon) nitem += 1;
  for (size_t i = 0; i < id.ni() - 1; ++i)
  {
    datac cdata(1 + nitem);
    stringstream ostr;
    ostr << spos->first << "-" << id.item(i + 1);
    cdata.assign(0, ostr.str());
    if (spos->second[i].n() < 1)
    { for (int j = 1; j <= nitem; ++j) cdata.assign(j, BLANKs); }
    else
    {
      spos->second[i].sort(opt.ol(), opt.olc());
      cdata.assign(1, (double)(spos->second[i].n()));
      cdata.assign(2, spos->second[i].mean());
      cdata.assign(3, spos->second[i].sd());
      cdata.assign(4, spos->second[i].se());
      cdata.assign(5, spos->second[i].min());
      cdata.assign(6, spos->second[i].q1());
      cdata.assign(7, spos->second[i].med());
      cdata.assign(8, spos->second[i].q3());
      cdata.assign(9, spos->second[i].max());
      cdata.assign(10, spos->second[i].max() - spos->second[i].min());
      double cv = (spos->second[i].mean() == 0.0) ? 0.0 : spos->second[i].sd() / spos->second[i].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, spos->second[i]);
      cdata.assign(12, g1);
      cdata.assign(13, g2);
      cdata.assign(14, spos->second[i].sum());
      if (opt.ol() != olnon)
      {
        double fn = spos->second[i].nol();
        cdata.assign(15, fn / (double)(spos->second[i].n()));
      }
    }
    cdata.put(fout);
  }
}

int sstatclass::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(id.ni()) != 0) throw option_error("実行時オプションエラー！");
  return nc;
}

void sstatclass::put(option& opt) // 基礎統計量を出力する
{
  ofstream foutdns;
  if (opt.dnsname() != "")      // DANSデータ出力
  {
    if (ofopen(foutdns, (opt.dnsname()).c_str()) == false) return; // ファイルオープン失敗
    put_dnsid(foutdns, opt);
  }
  varclass var(id.ni());        // 計算用変数クラス
  var.cal(opt, datalst);
  if (var.nmap() < 2) return;   // データなし
  size_t i = 1;
  for (mapsvecs_itr sitr = (var.ssvecmap()).begin(); sitr != (var.ssvecmap()).end(); ++sitr, ++i)
  {
    if (i < var.nmap())
    {
      if (opt.res("si"))
      {
        put_resid(cout, sitr->first);
        put_tab(cout, opt, sitr->second);
      }
    }
    else if (i > 2)
    {
      if (opt.res("st"))
      {
        put_resid(cout, sitr->first);
        put_tab(cout, opt, sitr->second);
      }
    }
    if (opt.dnsname() != "") put_dns(foutdns, opt, sitr);
  }
}

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