// ttest.C Version 7.44  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-29 15:08:55 tac>
//
// 対応のあるt検定(1標本t検定)を多時期で行う。
// 3時期以上については、本来は二元配置分散分析を行うべきだが、
// 諸般の事情により多時期についても計算している。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法  :g++ -O -s -W ttest ttest.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 =            // プログラム使用法表示用
{
  "ttest  Version 7.44          Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:ttest [<オプション> …]  [<データファイル> …]",
  "機能:基準項目と各項目との対応のあるt検定をする。",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定する。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -comp[={-|+}]:欠損のないケースだけ集計指定。(-comp:-comp=+)",
  "  -f=<ファイル名>:指定されたファイルからオプションを読み込む。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示する。(-h:-h=man)",
  "  -key=<基準項目番号>:基準項目の番号を指定する。(オプション省略時:-key=1)",
  "  -ngp={dbar|dci|ddata|dmsd|dmse|vbar|vci|vdata|vmsd|vmse}[,…]",
  "      :Ngraph用データの指定。",
  "    dbar:変化量バーチャート  dci:変化量平均値±CI  ddata:変化量生データ",
  "    dmsd:変化量平均値±SD  dmse:変化量平均値±SE",
  "    vbar:実測値バーチャート  vci:実測値平均値±CI  vdata:実測値生データ",
  "    vmsd:実測値平均値±SD  vmse:実測値平均値±SE",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データを出力するファイル名を指定する。",
  "    オプション省略時:データ出力せず  -ngpname:-ngpname=ttest_ngp.dat",
  "  -res=[-][{def|all|astat|m|pstat|tt}][,…]:解析結果の出力項目を指定する。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    astat:全データの基礎統計量  m:多重検定",
  "    pstat:対応のあるデータの基礎統計量  tt:対応のあるt検定",
  "    オプション省略時:データに合わせて自動判定",
  "  -t=<時間1>,<時間2>[,<時間3>,…]:Ngraph用に項目の時間を指定する。"
};

const string Title = "=== 2時期の平均値の比較 ==="; // 表題
const string NGPNAMEdef = "ttest_ngp.dat"; // 省略時Ngraphファイル名

enum ngptype                    // Ngraph用出力データの種類
{
  ngpnon,                       // 出力なし
  ngpdbar,                      // 変化量バーチャート
  ngpdci,                       // 変化量平均値±信頼区間
  ngpddata,                     // 変化量の生データ
  ngpdmsd,                      // 変化量平均値±SD
  ngpdmse,                      // 変化量平均値±SE
  ngpvbar,                      // 実測値バーチャート
  ngpvci,                       // 実測値平均値±信頼区間
  ngpvdata,                     // 実測値の生データ
  ngpvmsd,                      // 実測値平均値±SD
  ngpvmse,                      // 実測値平均値±SE
  ngpetc                        // その他
};

class option                    // オプション集合体
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  bool v_comp;                  // 欠損のないケースだけ集計フラッグ
  int v_key;                    // 基準項目番号
  ngptype v_ngp;                // Ngraph用出力データの種類
  string v_ngpname;             // Ngraphファイル名
  lstr ngplst;                  // Ngraph出力リスト
  mapsb ngpmap;                 // Ngraph出力の有無マップ
  lstr reslst;                  // 結果出力リスト
  mapsb resmap;                 // 結果出力の有無マップ
  lstr tlst;                    // 項目時間リスト
  vecdbl tvec;                  // 項目時間ベクトル
  testtype v_tt;                // t検定出力フラッグ
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; }
  bool comp() const             // 欠損のないケースだけ集計フラッグを返す
  { return v_comp; }
  int key() const               // 基準項目番号を返す
  { return v_key; }
  ngptype ngp() const           // Ngraph用出力データの種類を返す
  { return v_ngp; }
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  bool ngp(const string& code); // Ngraph出力フラッグを返す
  bool res(const string& code); // 分析結果出力フラッグを返す
  double t(size_t idx) const    // 項目時間を返す
  { return (idx < tvec.size()) ? tvec[idx] : BLANKf; }
  testtype tt() const           // t検定出力フラッグを返す
  { return v_tt; }
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int search_g(const string& gcode) const; // 群コードを検索する
  int check(size_t ni);         // オプションのチェックとセット
  int set_ngp();                // -ngpオプションのセット
  int set_res();                // -resオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_cc(CCdef), v_comp(false), v_key(1), v_ngp(ngpnon), v_ngpname(), ngplst(), ngpmap(), reslst(), resmap(), tlst(), tvec(), v_tt(testaut)
// コンストラクタ
{
  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] == "comp")   // 欠損のないケースだけ集計指定
  {
    if (opt[1].empty() || opt[1] == "+") v_comp = true;
    else if (opt[1] == "-") v_comp = false;
    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] == "key")     // 基準項目番号指定
  {
    char* ptr = nullptr;
    v_key = strtol(opt[1].c_str(), &ptr, 0);
    if (v_key < 1) err = badopt(pname(), optstr, USAGE, 1.0, (double)INT_MAX);
  }
  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] == "res")     // 分析結果の出力指定
  {
    if (strsplit(reslst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "t")       // 項目時間の指定
  {
    if (strsplit(tlst, 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 (v_key > (int)ni) return disperr(pname(), "基準項目番号が変やんか！");
  tvec.resize(ni);
  for (size_t j = 0; j < ni; ++j) tvec[j] = (double)(j + 1);
  if (settvec(tvec, tlst, ni) < 0) 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["dbar"] = ngpmap["dci"] = ngpmap["ddata"] = ngpmap["dmsd"] = ngpmap["dmse"] = ngpmap["vbar"] = ngpmap["vci"] = ngpmap["vdata"] = ngpmap["vmsd"] = ngpmap["vmse"] = true;
  if (ngplst.size() > 0)
  {
    if (v_ngpname == "") return disperr(pname(), "-ngpnameオプションの指定が必要！");
    ngpmap["dbar"] = ngpmap["dci"] = ngpmap["ddata"] = ngpmap["dmsd"] = ngpmap["dmse"] = ngpmap["vbar"] = ngpmap["vci"] = ngpmap["vdata"] = ngpmap["vmsd"] = ngpmap["vmse"] = 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["astat"] = resmap["pstat"] = resmap["tt"] = true;
  resmap["m"] = false;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["astat"] = resmap["pstat"] = resmap["tt"] = 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") // 省略時出力
      {
        resmap["def"] = resmap["astat"] = resmap["pstat"] = resmap["tt"] = true;
      }
      else if (resstr == "all") // 全て出力
      {
        resmap["def"] = resmap["astat"] = resmap["cc"] = resmap["m"] = resmap["pstat"] = resmap["tt"] = true;
      }
      else if (respos == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適当！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

class varclass                  // 変数クラス
{
private:
  int v_nc;                       // ケース数
  vector<statv<double>> astatvec; // 全データ用基礎統計量ベクトル
  vector<statv<double>> pstatvec; // 対データ用基礎統計量ベクトル
  vector<statv<double>> kstatvec; // 基準項目用基礎統計量ベクトル
  vector<statv<double>> dstatvec; // 変化量データ用基礎統計量ベクトル
public:
  varclass(size_t ni) : v_nc(0), astatvec(ni), pstatvec(ni), kstatvec(ni), dstatvec(ni)
  // デフォールトコンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  int cal(option& opt, ldata& datalst); // データを集計して変数に保存する
  int nc() const                // ケース数を返す
  { return v_nc; }
  int an(size_t idx) const      // 全データの例数を返す
  { return (idx < astatvec.size()) ? astatvec[idx].n() : 0; }
  double amean(size_t idx) const // 全データの平均値を返す
  { return (idx < astatvec.size()) ? astatvec[idx].mean() : 0.0; }
  double avr(size_t idx) const  // 全データの不偏分散を返す
  { return (idx < astatvec.size()) ? astatvec[idx].vr() : 0.0; }
  double asd(size_t idx) const  // 全データの標準偏差を返す
  { return (idx < astatvec.size()) ? astatvec[idx].sd() : 0.0; }
  double ase(size_t idx) const  // 全データの標準誤差を返す
  { return (idx < astatvec.size()) ? astatvec[idx].se() : 0.0; }
  int pn(size_t idx) const      // 対データの例数を返す
  { return (idx < pstatvec.size()) ? pstatvec[idx].n() : 0; }
  double pmean(size_t idx) const // 対データの平均値を返す
  { return (idx < pstatvec.size()) ? pstatvec[idx].mean() : 0.0; }
  double pvr(size_t idx) const  // 対データの不偏分散を返す
  { return (idx < pstatvec.size()) ? pstatvec[idx].vr() : 0.0; }
  double psd(size_t idx) const  // 対データの標準偏差を返す
  { return (idx < pstatvec.size()) ? pstatvec[idx].sd() : 0.0; }
  double pse(size_t idx) const  // 対データの標準誤差を返す
  { return (idx < pstatvec.size()) ? pstatvec[idx].se() : 0.0; }
  int kn(size_t idx) const      // 規準項目の例数を返す
  { return (idx < kstatvec.size()) ? kstatvec[idx].n() : 0; }
  double kmean(size_t idx) const // 規準項目の平均値を返す
  { return (idx < kstatvec.size()) ? kstatvec[idx].mean() : 0.0; }
  double kvr(size_t idx) const  // 規準項目の不偏分散を返す
  { return (idx < kstatvec.size()) ? kstatvec[idx].vr() : 0.0; }
  double ksd(size_t idx) const  // 規準項目の標準偏差を返す
  { return (idx < kstatvec.size()) ? kstatvec[idx].sd() : 0.0; }
  double kse(size_t idx) const  // 規準項目の標準誤差を返す
  { return (idx < kstatvec.size()) ? kstatvec[idx].se() : 0.0; }
  int dn(size_t idx) const      // 変化量の例数を返す
  { return (idx < dstatvec.size()) ? dstatvec[idx].n() : 0; }
  double dmean(size_t idx) const // 変化量の平均値を返す
  { return (idx < dstatvec.size()) ? dstatvec[idx].mean() : 0.0; }
  double dvr(size_t idx) const  // 変化量の不偏分散を返す
  { return (idx < dstatvec.size()) ? dstatvec[idx].vr() : 0.0; }
  double dsd(size_t idx) const  // 変化量の標準偏差を返す
  { return (idx < dstatvec.size()) ? dstatvec[idx].sd() : 0.0; }
  double dse(size_t idx) const  // 変化量の標準誤差を返す
  { return (idx < dstatvec.size()) ? dstatvec[idx].se() : 0.0; }
};

int varclass::cal(option& opt, ldata& datalst)
// データを集計して変数に保存する
{
  v_nc = 0;
  bool ok = true;
  for (auto& cdatai : datalst)
  {
    if (opt.comp())             // 欠損のないケースだけ集計
    {
      ok = true;
      for (size_t i = 0; ok && i < cdatai.size(); ++i)
      {
        if (cdatai.fdata(i) == BLANKf) ok = false;
      }
    }
    if (ok)
    {
      ++v_nc;
      double xkey = cdatai.fdata(opt.key() - 1);
      for (size_t i = 0; i < cdatai.size(); ++i)
      {
        double x = cdatai.fdata(i);
        if (x != BLANKf)
        {
          astatvec[i].sigma(x);
          if (xkey != BLANKf)
          {
            pstatvec[i].sigma(x);  kstatvec[i].sigma(xkey);
            dstatvec[i].sigma(x - xkey);
          }
        }
      }
    }
  }
  return v_nc;                  // 有効なケース数を返す
}

class ttestclass                // ttestメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout);
  // 解析結果のタイトルとデータ識別情報を出力する
  void put_astat(ostream& fout, option& opt, varclass& var) const;
  // 全データの基礎統計量を出力する
  void put_ttest(ostream& fout, option& opt, varclass& var);
  // 対応のあるデータの基礎統計量とt検定の結果を出力する
  void plot_ngraph(option& opt, varclass& var); // Ngraph用データを出力する
  void plot_ngpid(ofstream& fout, option& opt);
  // Ngraph用データのタイトルとデータ識別情報を出力する
  void plot_bar(ofstream& fout, option& opt, bool diff, varclass& var) const;
  // Ngraph用にバーチャートデータを出力する
  void plot_ci(ofstream& fout, option& opt, bool diff, varclass& var) const;
  // Ngraph用に平均値±信頼区間を出力する
  void plot_data(ofstream& fout, option& opt, bool diff);
  // Ngraph用に生データを出力する
  void plot_msde(ofstream& fout, option& opt, bool sd, bool diff, varclass& var) const;
  // Ngraph用に平均値±SD/SEを出力する
public:
  ttestclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~ttestclass()                 // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み集計する
  void put(option& opt);        // 対応のあるt検定結果を出力する
};

void ttestclass::put_resid(ostream& fout)
// 集計結果のタイトルとデータ識別情報を出力する
{
  puttitle(fout, Title);
  fout << "\n";
  id.put(fout);
}

void ttestclass::put_astat(ostream& fout, option& opt, varclass& var) const
// 全データの基礎統計量を出力する
{
  constexpr int LENline = 77; // 区切り線の長さ
  fout << "\n○全データの基礎統計量";
  if (opt.comp()) fout << "(欠損がないケースだけ集計)";
  fout << "\n";
  putline(fout, LENline);
  for (size_t i = 0; i < id.ni(); ++i)
  {
    fout << "項目" << setw(4) << i + 1 << ":";
    putmsd(fout, var.an(i), var.amean(i), var.avr(i));
  }
  putline(fout, LENline);
}

void ttestclass::put_ttest(ostream& fout, option& opt, varclass& var)
// 対応のあるデータの基礎統計量とt検定の結果を出力する
{
  constexpr int LENline = 77; // 区切り線の長さ
  fout << "\n○対応のあるデータと対応のあるt検定(1標本t検定、one sample t-test)\n";
  for (size_t i = 0; i < id.ni(); ++i)
  {
    int ino = i + 1;
    if (var.pn(i) > 0 && ino != opt.key())
    {
      fout << "\n項目" << setw(4) << ino << ":" << id.item(i) << "\n";
      if (opt.res("pstat"))
      {
        putline(fout, LENline);
        fout << "項目" << setw(4) << ino << ":";
        putmsd(fout, var.pn(i), var.pmean(i), var.pvr(i));
        fout << "項目" << setw(4) << opt.key() << ":";
        putmsd(fout, var.kn(i), var.kmean(i), var.kvr(i));
      }
      if (opt.res("tt") || opt.res("m"))
      {
        putline(fout, LENline);
        fout << "変化量  :";
        putmsd(fout, var.dn(i), var.dmean(i), var.dvr(i));
      }
      int dfd = var.dn(i) - 1;
      if (dfd > 0)
      {
        double md = var.dmean(i), sed = var.dse(i);
        double t = (sed > 0.0) ? md / sed : 0.0;
        ldbl pr = ttop(t, dfd);
        if (opt.res("tt"))
        {
          fout << "  t=" << t << "  自由度=" << dfd << "  有意確率:p=" << pr << pind(pr) << "\n";
        }
        double ci = calmcci(mcfis, opt.cc(), 2, 1, dfd, 1, sed * sed);
        fout << "  変化量平均値の" << opt.cc() * 100.0 << "%信頼区間=";
        fout << md << "±" << ci << "(" << md - ci << "-" << md + ci << ")\n";
        if ((id.ni() > 2 && opt.res("def")) || opt.res("m"))
        {
          int dc = id.ni() - 1;
          pr = pr * (ldbl)dc;
          if (pr > 1.0L) pr = 1.0L;
          fout << "  Bonferroni型多重検定の有意確率:p=" << pr << pind(pr) << "\n";
          ci = calmcci(mcbon, opt.cc(), 2, 1, dfd, dc, sed * sed);
          fout << "  変化量平均値の" << opt.cc() * 100.0 << "%同時信頼区間=";
          fout << md << "±" << ci << "(" << md - ci << "-" << md + ci << ")\n";
        }
        putline(fout, LENline);
      }
    }
  }
}

void ttestclass::plot_ngraph(option& opt, varclass& var)
// Ngraph用データを出力する
{
  ofstream fout;
  if (ofopen(fout, (opt.ngpname()).c_str()) == false) return; // ファイルオープン失敗
  bool delml = false;
  if (opt.ngp("ddata"))         // 変化量生データ
  {
    fout << "-s1 -zx0," << opt.t(id.ni() - 1) + 1.0 << ",1 -dline -w40\n";
    plot_ngpid(fout, opt);
    plot_data(fout, opt, true);
    delml =true;
  }
  if (opt.ngp("dbar"))          // 変化量バーチャート
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -zx0," << opt.t(id.ni() - 1) + 1.0 << ",1 -drectangle_solid_fill -x1 -y3\n";
      delml =true;
    }
    plot_ngpid(fout, opt);
    plot_bar(fout, opt, true, var);
  }
  if (opt.ngp("dci"))           // 変化量平均値±CI
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -zx0," << opt.t(id.ni() - 1) + 1.0 << ",1 -dline -w80\n";
      delml =true;
    }
    plot_ngpid(fout, opt);
    plot_ci(fout, opt, true, var);
  }
  if (opt.ngp("dmsd"))          // 変化量平均値±SD
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -zx0," << opt.t(id.ni() - 1) + 1.0 << ",1 -dline -w80\n";
      delml =true;
    }
    plot_ngpid(fout, opt);
    plot_msde(fout, opt, true, true, var);
  }
  if (opt.ngp("dmse"))          // 変化量平均値±SE
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -zx0," << opt.t(id.ni() - 1) + 1.0 << ",1 -dline -w80\n";
      delml =true;
    }
    plot_ngpid(fout, opt);
    plot_msde(fout, opt, false, true, var);
  }
  if (opt.ngp("vdata"))         // 実測値生データ
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -zx0," << opt.t(id.ni() - 1) + 1.0 << ",1 -dline -w40\n";
      delml =true;
    }
    plot_ngpid(fout, opt);
    plot_data(fout, opt, false);
  }
  if (opt.ngp("vbar"))          // 実測値バーチャート
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -zx0," << opt.t(id.ni() - 1) + 1.0 << ",1 -drectangle_solid_fill -x1 -y3\n";
      delml =true;
    }
    plot_ngpid(fout, opt);
    plot_bar(fout, opt, false, var);
  }
  if (opt.ngp("vci"))           // 実測値平均値±CI
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -zx0," << opt.t(id.ni() - 1) + 1.0 << ",1 -dline -w80\n";
      delml =true;
    }
    plot_ngpid(fout, opt);
    plot_ci(fout, opt, false, var);
  }
  if (opt.ngp("vmsd"))          // 実測値平均値±SD
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -zx0," << opt.t(id.ni() - 1) + 1.0 << ",1 -dline -w80\n";
      delml =true;
    }
    plot_ngpid(fout, opt);
    plot_msde(fout, opt, true, false, var);
  }
  if (opt.ngp("vmse"))          // 実測値平均値±SE
  {
    if (delml) fout << NGPdelm << "\n";
    else fout << "-s1 -zx0," << opt.t(id.ni() - 1) + 1.0 << ",1 -dline -w80\n";
    plot_ngpid(fout, opt);
    plot_msde(fout, opt, false, false, var);
  }
}

void ttestclass::plot_ngpid(ofstream& fout, option& opt)
// Ngraph用データのタイトルとデータ識別情報を出力する
{
  fout << "%";  puttitle(fout, Title);
  fout << "%\n";
  id.put(fout, formngp);
  if (opt.comp()) fout << "%欠損がないケースだけ集計\n";
}

void ttestclass::plot_bar(ofstream& fout, option& opt, bool diff, varclass& var) const
// Ngraph用にバーチャートデータを出力する
{
  double t0 = opt.t(0);
  if (diff == false)            // 実測値プロット
  {
    fout << "%全体の実測値平均値のバーチャート\n";
    for (size_t i = 0; i < id.ni(); ++i)
    {
      double ti = opt.t(i);
      fout << ti - 0.2 << " 0 " << ti + 0.2 << " " << var.amean(i) << "\n";
    }
    fout << NGPdelm2 << "\n";
  }
  if (diff) fout << "%対応のある変化量";
  else fout << "%対応のある実測値";
  fout << "平均値のバーチャート\n";
  for (size_t i = 1; i < id.ni(); ++i)
  {
    double ti = opt.t(i), m0 = var.kmean(i), mi = var.pmean(i);
    if (diff) { m0 = 0.0;  mi = var.dmean(i); }
    if (i > 1) fout << NGPdelm2 << "\n";
    fout << "%項目1-" << i + 1 << "\n";
    fout << t0 - 0.2 << " 0 " << t0 + 0.2 << " " << m0 << "\n";
    fout << ti - 0.2 << " 0 " << ti + 0.2 << " " << mi << "\n";
  }
}

void ttestclass::plot_ci(ofstream& fout, option& opt, bool diff, varclass& var) const
// Ngraph用に平均値±信頼区間を出力する
{
  if (diff == false)            // 実測値プロット
  {
    fout << "%全体の実測値平均値±" << opt.cc() * 100.0 << "%CI(-zx0," << opt.t(id.ni() - 1) + 1.0 << ",1 -derrorbar_y -w0 -o500)\n";
    for (size_t i = 0; i < id.ni(); ++i)
    {
      double t = ptot(1.0 - opt.cc(), var.an(i) - 1);
      double ci = t * var.ase(i);
      fout << opt.t(i) << " " << var.amean(i) << " " << ci << " " << -ci << "\n";
    }
    fout << NGPdelm2 << "\n";
  }
  if (diff) fout << "%対応のある変化量";
  else fout << "%対応のある実測値";
  fout << "平均値±" << opt.cc() * 100.0 << "%CI\n";
  for (size_t i = 1; i < id.ni(); ++i)
  {
    double m0 = var.kmean(i), se0 = var.kse(i);
    double mi = var.pmean(i), sei = var.pse(i);
    int df = var.kn(i) - 1;
    if (diff)
    {
      m0 = se0 = 0.0;
      mi = var.dmean(i);  sei = var.dse(i);
    }
    if (i > 1) fout << NGPdelm2 << "\n";
    fout << "%項目1-" << i + 1 << "\n";
    double t = (0.0 < opt.cc() && opt.cc() < 1.0) ? ptot(1.0 - opt.cc(), df) : 0.0;
    double ci = t * se0;
    fout << opt.t(0) << " " << m0 << " " << ci << " " << -ci << "\n";
    ci = t * sei;
    fout << opt.t(i) << " " << mi << " " << ci << " " << -ci << "\n";
  }
}

void ttestclass::plot_data(ofstream& fout, option& opt, bool diff)
// Ngraph用に生データを出力する
{
  if (diff) fout << "%変化量生データ\n";
  else fout << "%実測値生データ\n";
  bool putdelml = false;
  for (auto& cdatai : datalst)
  {
    if (putdelml) fout << NGPdelm2 << "\n";
    else putdelml = true;
    double x0 = cdatai.fdata(0);
    if (x0 != BLANKf)
    {
      if (diff) fout << opt.t(0) << " " << 0.0 << "\n";
      else fout << opt.t(0) << " " << x0 << "\n";
    }
    for (size_t i = 1; i < cdatai.size(); ++i)
    {
      double xi = cdatai.fdata(i);
      if (xi != BLANKf)
      {
        if (diff)
        {
          if (x0 != BLANKf) fout << opt.t(i) << " " << xi - x0 << "\n";
        }
        else
        { fout << opt.t(i) << " " << xi << "\n"; }
      }
    }
  }
}

void ttestclass::plot_msde(ofstream& fout, option& opt, bool sd, bool diff, varclass& var) const
// Ngraph用に平均値±SD/SEを出力する
{
  string strsde = sd ? "SD" : "SE";
  if (diff == false)            // 実測値プロット
  {
    fout << "%全体の実測値平均値±" << strsde << "(-zx0," << opt.t(id.ni() - 1) + 1.0 << ",1 -derrorbar_y -w40 -o500)\n";
    for (size_t i = 0; i < id.ni(); ++i)
    {
      double sdei = sd ? var.asd(i) : var.ase(i);
      fout << opt.t(i) << " " << var.amean(i) << " " << sdei << " " << -sdei << "\n";
    }
    fout << NGPdelm2 << "\n";
  }
  if (diff) fout << "%対応のある変化量";
  else fout << "%対応のある実測値";
  fout << "平均値±" << strsde << "\n";
  for (size_t i = 1; i < id.ni(); ++i)
  {
    if (i > 1) fout << NGPdelm2 << "\n";
    fout << "%項目1-" << i + 1 << "\n";
    double m0 = var.kmean(i), mi = var.pmean(i);
    double sde0 = sd ? var.ksd(i) : var.kse(i);
    double sdei = sd ? var.psd(i) : var.pse(i);
    if (diff)
    {
      m0 = sde0 = 0.0;
      mi = var.dmean(i);  sdei = sd ? var.dsd(i) : var.dse(i);
    }
    fout << opt.t(0) << " " << m0 << " " << sde0 << " " << -sde0 << "\n";
    fout << opt.t(i) << " " << mi << " " << sdei << " " << -sdei << "\n";
  }
}

int ttestclass::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 ttestclass::put(option& opt) // 対応のあるt検定結果を出力する
{
  varclass var(id.ni());        // 計算用変数クラス
  var.cal(opt, datalst);
  if (opt.res("res"))           // 計算結果出力
  {
    put_resid(cout);
    if (var.nc() < 1)
    {
      cout << "※有効なデータなし！\n\n" << endl;
      return;
    }
    if (opt.res("astat")) put_astat(cout, opt, var);
    if (opt.res("pstat") || opt.res("tt")) put_ttest(cout, opt, var);
    cout << "\n" << endl;
  }
  if (var.nc() > 0 && opt.ngpname() != "") plot_ngraph(opt, var); // グラフデータ出力
}

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