// anovam.C Version 7.45  Copyleft (c) 杉本典夫
// Time-stamp: <2026-04-08 13:50:01 tac>
//
// 例数、平均値、SD、SEデータから対応の無い多群間で一元配置分散分析を行う。
// 2群の場合は2標本t検定も出力する。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法  :g++ -O -s -W anovam anovam.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 =            // プログラム使用法表示用
{
  "anovam  Version 7.45         Copyleft(Vo\\oV) 2026 杉本典夫",
  "書式:anovam  [<オプション> …]  [<データファイル> …]",
  "機能:例数、平均、SD又はSE又は分散で一元配置分散分析をするんじゃ。",
  "     2標本t検定も含むぞい。",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定するんじゃ。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -f=<ファイル名>:指定したファイルからオプションを読み込むんじゃ。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示するんじゃ。(-h:-h=man)",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データを出力するファイル名を指定する。",
  "    オプション省略時:Ngraph用データ出力せず  -ngpname:-ngpname=anovam_ngp.dat",
  "  -res=[-][{def|all|at|ft|m[{b|d|f|s|t}]|r|stat|tt|vt}][,…]",
  "      :解析結果の出力項目を指定するんじゃ。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    at:分散分析  ft:等分散性の両側F検定  m[{b|d|f|s|t}]:多重比較",
  "      b:Bonferroni型  d:Dunnett型  f:Fisher型  s:Scheffe型  t:Tukey型",
  "    r:相関比  stat:基礎統計量  tt:2標本t検定  vt:等分散性の検定",
  "    オプション省略時:データに合わせて自動判定"
};

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

class option                    // オプション集合体
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  mctype v_m;                   // 多重比較法
  string v_ngpname;             // 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; }
  mctype m() const              // 多重比較法を返す
  { return v_m; }
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  bool res(const string& code); // 分析結果出力フラッグを返す
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check();                  // オプションのチェックとセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_cc(CCdef), v_m(mctuk), v_ngpname(), 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] == "f")       // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  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                          // 間違ったオプション指定
  { err = nonopt(pname(), optstr, USAGE); }
  if (err > 0) throw help_except("ヘルプ表示");
  if (err < 0) throw option_error(string("オプション指定違反→") + optstr);
}

int option::check()             // オプションのセットとチェック
{
  resmap["res"] = resmap["def"] = resmap["ft"] = resmap["stat"] = resmap["tt"] = true;
  resmap["at"] = resmap["m"] = resmap["r"] = resmap["vt"] = false;
  if (reslst.size() > 0)
  {
     resmap["def"] = resmap["ft"] = resmap["stat"] = 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["ft"] = resmap["stat"] = resmap["tt"] = true; }
      else if (resstr == "all")
      {
        resmap["def"] = resmap["at"] = resmap["ft"] = resmap["m"] = resmap["r"] = resmap["stat"] = resmap["tt"] = resmap["vt"] = true;
      }
      else if (resstr[0] == 'm' && resstr.size() > 1)
      {
        switch (resstr[1])
        {
        case 'b' : v_m = mcbon;
          break;
        case 'd' : v_m = mcdun;
          break;
        case 'f' : v_m = mcfis;
          break;
        case 's' : v_m = mcsch;
          break;
        case 't' : v_m = mctuk;
          break;
        default : return disperr(pname(), "-res=" + resstr + "オプション指定が不適当！");
        }
        resmap["m"] = true;
      }
      else if (respos == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適当！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

class gstat                     // 群統計量クラス
{
private:
  string v_code;                // 群コード
  statv<double> v_ystat;        // 基礎統計量
public:
  gstat() : v_code(), v_ystat() // デフォールトコンストラクタ
  {}
  gstat(const string& code0, int n0, double mean0, double sd0);
  // コンストラクタ
  gstat(const gstat& grp) : v_code(grp.v_code), v_ystat(grp.v_ystat)
  // コピーコンストラクタ
  {}
  virtual ~gstat()              // デストラクタ
  {}
  const gstat& operator=(const gstat& gstat0); // コピー代入演算子
  int set(const string& code0, int n0, double mean0, double sd0);
  // データをセットする
  const string& code() const    // 群コードを返す
  { return v_code; }
  const statv<double>& ystat() const // 基礎統計量を返す
  { return v_ystat; }
};

gstat::gstat(const string& code0, int n0, double mean0, double sd0) : v_code(code0), v_ystat()
// コンストラクタ
{
  v_ystat.set(n0, mean0, sd0);
}

const gstat& gstat::operator=(const gstat& gstat0) // コピー代入演算子
{
  if (&gstat0 == this) return *this; // 自分自身への代入を防ぐ
  v_code = gstat0.v_code;  v_ystat = gstat0.v_ystat;
  return *this;
}

int gstat::set(const string& code0, int n0, double mean0, double sd0)
// データをセットする
{
  v_code = code0;
  return v_ystat.set(n0, mean0, sd0);
}

class anovamclass               // anovamメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  list<gstat> gstatlst;         // 群統計量リスト
  void get_stat(gstat& tgstat, const datac& cdata); // 統計量を取得する
  void put(ostream& fout, ofstream& foutngp, option& opt, int nd);
  // 集計結果を出力する
  void put_resid(ostream& fout, int nd);
  // 分析結果のタイトルとデータ識別情報を出力する
  double put_stat(ostream& fout, option& opt, statv<double>& tstat);
  // 群別基礎統計量の出力と累積集計
  void put_rratio(ostream& fout, option& opt, double sa, double sr, double st, int dfa, int dfr);
  // 相関比の出力
  void put_ttest(ostream& fout, option& opt, const statv<double>& tstat);
  // 2標本t検定とWelchの検定
  ldbl put_ftest(ostream& fout, option& opt, double f, int df1, int df2);
  // 分散比Fの検定と推定結果を出力し両側p値を返す
  void put_vtest(ostream& fout, const statv<double>& tstat) const;
  // 分散の均一性の検定結果出力
  void put_mcop(ostream& fout, option& opt, const statv<double>& tstat, double vr) const;
  // 2群間の多重比較出力
  void put_mcopd(ostream& fout, option& opt, const statv<double>& tstat, double vr) const;
  // Dunnett型多重比較出力
  void put_mdci(ostream& fout, option& opt, const statv<double>& tstat, double vr) const;
  // 2群間の平均値の差の信頼区間出力
  void put_mdcid(ostream& fout, option& opt, const statv<double>& tstat, double vr) const;
  // 2群間の平均値の差のDunnett型信頼区間出力
  void plot_ngpid(ofstream& fout, int nd);
  // Ngraph用データのタイトルとデータ識別情報を出力する
  void plot_ngp(ofstream& fout, option& opt); // Ngraph用データ出力
public:
  anovamclass() : fname(), id(), gstatlst() // コンストラクタ
  {}
  ~anovamclass()                // デストラクタ
  {}
  int ng() const                // 群数を返す
  { return (int)gstatlst.size(); }
  void getput(const char* c_fname, option& opt);
  // データを読み込んで集計結果を出力する
};

void anovamclass::get_stat(gstat& tgstat, const datac& cdata)
// 統計量を取得する
{
  if (cdata.sdata(0) == BLANKs) throw runtime_error("群コードの指定が間違っとるぞい！");
  int n = 0;
  double mean = BLANKf, sd = BLANKf, se = BLANKf, v = BLANKf;
  string str[2];                // 統計量指定str[0]:項目,str[1]:データ
  for (size_t j = 1; j < cdata.size(); ++j)
  {
    if (getopt(str, cdata.sdata(j)) > 0)
    {
      strlower(str[0]);
      char* ptr = nullptr;
      if (str[0] == "n") n = strtol((str[1]).c_str(), &ptr, 0);
      else if (str[0] == "m") mean = strtodf(str[1]);
      else if (str[0] == "sd") sd = strtodf(str[1]);
      else if (str[0] == "se") se = strtodf(str[1]);
      else if (str[0] == "v") v = strtodf(str[1]);
    }
  }
  if (n < 1 || mean == BLANKf) throw runtime_error("例数・平均値の指定が間違っとるぞい！");
  if (sd == BLANKf)             // sdなし
  {
    if (v != BLANKf) sd = sqrt(v);
    else if (se != BLANKf) sd = se * sqrt((double)n);
  }
  if (sd == BLANKf) throw runtime_error("バラツキの指定が間違っとるぞい！");
  tgstat.set(cdata.sdata(0), n, mean, sd);
}

void anovamclass::put(ostream& fout, ofstream& foutngp, option& opt, int nd)
// 集計結果を出力する
{
  statv<double> tstat;
  if (opt.res("res"))           // 計算結果出力
  {
    put_resid(fout, nd);
    if (ng() < 1)
    {
      fout << "※有効なデータなし！\n\n" << endl;
      return;
    }
    double sa = put_stat(fout, opt, tstat);
    if (ng() < 2)               // 分散分析不可能
    {
      fout << "※分散分析不可能！\n\n" << endl;
      return;
    }
    double st = tstat.ss();
    double sr = st - sa;
    int dfa = ng() - 1;
    int dft = tstat.n() - 1, dfr = tstat.n() - ng();
    if (sr < FDeps) sr = 0.0;
    if (ng() == 2)
    {
      if (opt.res("r")) put_rratio(fout, opt, sa, sr, st, dfa, dfr);
      if (opt.res("at"))        // 2群の時はatを指定された時だけ出力
      {
        if (opt.res("vt")) put_vtest(fout, tstat);
        fout << "\n・一元配置分散分析(one-way layout analysis of variance)\n";
        putatab(fout, "群(要因A)", sa, dfa, sr, dfr, st, dft);
      }
      if (opt.res("tt")) put_ttest(fout, opt, tstat);
    }
    else
    {
      if (opt.res("def") || opt.res("r")) put_rratio(fout, opt, sa, sr, st, dfa, dfr);
      if (opt.res("vt")) put_vtest(fout, tstat);
      if (opt.res("def") || opt.res("at")) // 3群以上の時は省略時も出力
      {
        fout << "\n・一元配置分散分析(one-way layout analysis of variance)\n";
        putatab(fout, "群(要因A)", sa, dfa, sr, dfr, st, dft);
      }
      if (opt.res("m"))
      {
        double vr = (dfr > 0) ? sr / dfr : 0.0;
        if (opt.m() == mcdun)   // Dunnett型
        {
          put_mcopd(fout, opt, tstat, vr);
          put_mdcid(fout, opt, tstat, vr);
        }
        else
        {
          put_mcop(fout, opt, tstat, vr);
          put_mdci(fout, opt, tstat, vr);
        }
      }
    }
    fout << "\n" << endl;
  }
  if (opt.ngpname() != "")      // Ngraph用データ出力
  {
    plot_ngpid(foutngp, nd);
    plot_ngp(foutngp, opt);
  }
}

void anovamclass::put_resid(ostream& fout, int nd)
// 分析結果のタイトルとデータ識別情報を出力する
{
  if (id.ni() < 2) id.set(fname, 2);
  puttitle(fout, Title);
  fout << "\n";
  id.put_title(fout);
  id.put_stra(fout);
  fout << "群項目(要因A):" << id.item(0) << "\n";
  fout << "集計項目     :" << id.item(1) << "\n";
}

double anovamclass::put_stat(ostream& fout, option& opt, statv<double>& tstat)
// 群別基礎統計量の出力と累積集計
{
  constexpr int LENline = 77;   // 区切り線の長さ
  constexpr int LENhed = 9;     // ヘッダの長さ
  if (opt.res("stat"))
  {
    fout << "\n群       :群別基礎統計量\n";
    putline(fout, LENline);
  }
  double sa = 0.0;
  for (auto& gstati : gstatlst)
  {
    int n = gstati.ystat().n();
    double m = gstati.ystat().mean();
    int npad = gstati.code().size() - strclen(gstati.code(), CHARwide);
    if (opt.res("stat"))
    {
      fout << left << setw(LENhed + npad) << gstati.code() << right << ":";
      putmsd(fout, n, m, gstati.ystat().vr());
    }
    tstat += gstati.ystat();  sa += n * m * m;
  }
  if (opt.res("stat"))
  {
    putline(fout, LENline);
    string tstr("全体");
    int npad = tstr.size() - strclen(tstr, CHARwide);
    fout << left << setw(LENhed + npad) << tstr << right << ":";
    putmsd(fout, tstat.n(), tstat.mean(), tstat.vr());
  }
  sa -= tstat.n() * tstat.mean() * tstat.mean();
  if (sa < FDeps) sa = 0.0;
  return sa;                    // 要因Aの平方和を返す
}

void anovamclass::put_rratio(ostream& fout, option& opt, double sa, double sr, double st, int dfa, int dfr)
// 相関比の出力
{
  double r2 = sa / st;
  fout << "\n・相関比:η=" << sqrt(r2) << "  寄与率:η2乗=" << r2 << "\n";
  double va = (dfa > 0) ? sa / (double)dfa : 0.0;
  double vr = (dfr > 0) ? sr / (double)dfr : 0.0;
  double f = (vr > 0.0) ? va / vr : 0.0;
  double crl = 0.0, cru = 0.0;
  calcrci(crl, cru, opt.cc(), dfa, dfr, f);
  fout << "  相関比の" << opt.cc() * 100.0 << "%信頼区間:下限=" << crl << " 上限=" << cru << "\n";
}

void anovamclass::put_ttest(ostream& fout, option& opt, const statv<double>& tstat)
// 2標本t検定とWelchの検定
{
  list<gstat>::const_iterator gitr = gstatlst.begin();
  const statv<double>* ystat1 = &(gitr->ystat());
  ++gitr;
  const statv<double>* ystat2 = &(gitr->ystat());
  int n1 = ystat1->n(), n2 = ystat2->n();
  double ss1 = ystat1->ss(), ss2 = ystat2->ss();
  double v1 = ystat1->vr(), v2 = ystat2->vr();
  double f = (v2 > 0.0) ? v1 / v2 : 0.0;
  int df1 = n1 - 1, df2 = n2 - 1;
  if (v1 < v2)
  {
    f = (v1 > 0.0) ? v2 / v1 : 0.0;
    df1 = n2 - 1;  df2 = n1 - 1;
  }
  ldbl pr = put_ftest(fout, opt, f, df1, df2);
  fout << "\n・対応のないt検定(2標本t検定、「→」の付いた方を採用して下さい)\n";
  string sel1 = "  →", sel2 = "    ";
  if (pr <= 0.05L) { sel1 = "    ";  sel2 = "  →"; }
  double d = ystat2->mean() - ystat1->mean();
  int df = n1 + n2 - 2;
  double vd = (ss1 + ss2) / df * (double)(n1 + n2) / (double)n1 / (double)n2;
  double t = (sgn(vd) > 0) ? d / sqrt(vd) : 0.0;
  pr = ttop(t, df);
  double dsd = sqrt((ss1 + ss2) / (double)df);
  double dc = 0.0;
  if (dsd < FDeps) dsd = 0.0;
  else dc = d / dsd;
  fout << "  平均値の差=" << d << "  平均値の差のSD=" << dsd << "  効果量=" << dc << "\n";
  fout << sel1 << "等分散  t検定(=分散分析) :t値=" << t;
  fout << "  自由度=" << df << "  p=" << pr << pind(pr) << "\n";
  double ci = calmcci(mcfis, opt.cc(), 2, 1, df, 1, vd);
  fout << sel1 << "平均値の差の" << opt.cc() * 100.0 << "%信頼区間=";
  fout << d << "±" << ci << "(" << d - ci << "-" << d + ci << ")\n";
  if (df1 < 1 || df2 < 1 || sgn(v1 * v2) <= 0) return;
  double vn1 = v1 / (double)n1, vn2 = v2 / (double)n2;
  vd = vn1 + vn2;
  df = (int)(vd * vd / (vn1 * vn1 / (double)(n1 - 1) + vn2 * vn2 / (double)(n2 -1)));
  t = d / sqrt(vd);
  pr = ttop(t, df);
  fout << sel2 << "不等分散t検定(Welchの検定):t値=" << t;
  fout << "  自由度=" << df << "  p=" << pr << pind(pr) << "\n";
  ci = calmcci(mcfis, opt.cc(), 2, 1, df, 1, vd);
  fout << sel2 << "平均値の差の" << opt.cc() * 100.0 << "%信頼区間=";
  fout << d << "±" << ci << "(" << d - ci << "-" << d + ci << ")\n";
}

ldbl anovamclass::put_ftest(ostream& fout, option& opt, double f, int df1, int df2)
// 分散比Fの両側検定と推定結果を出力し両側p値を返す
{
  ldbl pr = ftop(f, df1, df2);
  pr = (pr < 0.5L) ? 2.0L * pr : 1.0L; // 両側有意確率にする
  if (opt.res("ft"))
  {
    fout << "\n・等分散性の両側F検定\n";
    fout << "  F=" << f << " 第1自由度=" << df1 << " 第2自由度=" << df2 << "  両側p=" << pr << pind(pr) << "\n";
    double fl = 0.0, fu = 0.0;
    calfci(fl, fu, opt.cc(), df1, df2, f);
    fout << "  分散比Fの" << opt.cc() * 100.0 << "%信頼区間：下限=" << fl << " 上限=" << fu << "\n";
  }
  return pr;                    // 両側p値を返す
}

void anovamclass::put_vtest(ostream& fout, const statv<double>& tstat) const
// 分散の均一性の検定結果出力
{
  int d = -ng();
  double w = 0.0, ss = 0.0, lnss = 0.0;
  for (auto& gstati : gstatlst)
  {
    int df = gstati.ystat().n() - 1;
    d += gstati.ystat().n();  w += 1.0 / (double)df;
    ss += (double)df * gstati.ystat().vr();
    lnss += (double)df * log(gstati.ystat().vr());
  }
  int df = ng() - 1;
  double m = (d > 0 && sgn(ss) > 0) ? (double)d * log(ss / (double)d) - lnss : 0.0;
  double c = (df > 0 && d > 0) ? 1.0 + 1.0 / 3.0 / (double)df * (w - 1.0 / (double)d) : 0.0;
  double csqr = (sgn(m) >= 0 && sgn(c) > 0) ? m / c : 0.0;
  ldbl pr = ctop(csqr, df);
  fout << "\n・等分散性の検定(Bartlettの検定):χ2乗=" << csqr;
  fout << "  自由度=" << df << "  p=" << pr << pind(pr) << "\n";
}

void anovamclass::put_mcop(ostream& fout, option& opt, const statv<double>& tstat, double vr) const
// 2群間の多重比較出力
{
  if (vr <= 0.0) return;        // 出力なし
  int df1 = 1, df2 = tstat.n() - ng(), dc = 1;
  string method = "Tukey";
  switch (opt.m())
  {
  case mcbon : method = "Bonferroni";  dc = ng() * (ng() - 1) / 2;
    break;
  case mcfis : method = "Fisher";
    break;
  case mcsch : method = "Scheffe";  df1 = ng() - 1;
  }
  fout << "\n・群(要因A)の" << method << "型多重比較(" << method << " type multiple comparison)\n";
  if (opt.m() == mctuk)
  { fout << "群      - 群               q値        群数      自由度   有意確率p値\n"; }
  else
  { fout << "群      - 群               F値   第1自由度   第2自由度   有意確率p値\n"; }
  constexpr int LENline = 72;   // 区切り線の長さ
  putline(fout, LENline);
  list<gstat>::const_iterator gitr1 = gstatlst.begin();
  for (int i = 0; i < ng() - 1; ++i, ++gitr1)
  {
    string codei(gitr1->code());
    int npadi = codei.size() - strclen(codei, CHARwide);
    int ni = (gitr1->ystat()).n();
    double mi = (gitr1->ystat()).mean();
    list<gstat>::const_iterator gitr2 = gitr1;
    ++gitr2;
    while (gitr2 != gstatlst.end())
    {
      int nj = (gitr2->ystat()).n();
      double dm = (gitr2->ystat()).mean() - mi;
      double f = dm * dm / df1 / (1.0 / (double)ni + 1.0 / (double)nj) / vr;
      double pr = calmcpr(opt.m(), f, ng(), df1, df2, dc);
      int npad = gitr2->code().size() - strclen(gitr2->code(), CHARwide);
      fout << left << setw(8 + npad) << gitr2->code();
      fout << "- ";
      fout << setw(8 + npadi) << codei << right << setw(12) << f;
      if (opt.m() == mctuk) fout << setw(12) << ng();
      else fout << setw(12) << df1;
      fout << setw(12) << df2 << setw(14) << pr << pind(pr) << "\n";
      ++gitr2;
    }
  }
  putline(fout, LENline);
}

void anovamclass::put_mcopd(ostream& fout, option& opt, const statv<double>& tstat, double vr) const
// Dunnett型多重比較出力
{
  if (vr <= 0.0) return;        // 出力なし
  int df1 = 1, df2 = tstat.n() - ng(), dc = 1;
  fout << "\n・群(要因A)のDunnett型多重比較(Dunnett type multiple comparison)\n";
  fout << "群      - 群               d値        群数      自由度   有意確率p値\n";
  constexpr int LENline = 72;   // 区切り線の長さ
  putline(fout, LENline);
  list<gstat>::const_iterator gitr = gstatlst.begin();
  string codei(gitr->code());
  int npadi = codei.size() - strclen(codei, CHARwide);
  double ni = (double)((gitr->ystat()).n());
  double mi = (gitr->ystat()).mean();
  ++gitr;
  for (int j = 1; j < ng(); ++j, ++gitr)
  {
    double nj = (double)((gitr->ystat()).n());
    double dm = (gitr->ystat()).mean() - mi;
    double f = dm * dm / df1 / (1.0 / ni + 1.0 / nj) / vr;
    double r = (ni == nj) ? 0.5 : sqrt(ni / (ni + nj) * nj / (ni + nj));
    double pr = calmcpr(opt.m(), f, ng(), df1, df2, dc, r);
    int npad = gitr->code().size() - strclen(gitr->code(), CHARwide);
    fout << left << setw(8 + npad) << gitr->code();
    fout << "- ";
    fout << setw(8 + npadi) << codei << right << setw(12) << f;
    fout << setw(12) << ng() << setw(12) << df2 << setw(14) << pr << pind(pr) << "\n";
  }
  putline(fout, LENline);
}

void anovamclass::put_mdci(ostream& fout, option& opt, const statv<double>& tstat, double vr) const
// 2群間の平均値の差の信頼区間出力
{
  if (vr <= 0.0) return;        // 出力なし
  int df1 = 1, df2 = tstat.n() - ng(), dc = 1;
  string method = "Tukey";
  switch (opt.m())
  {
  case mcbon : method = "Bonferroni";  dc = ng() * (ng() - 1) / 2;
    break;
  case mcfis : method = "Fisher";
    break;
  case mcsch : method = "Scheffe";  df1 = ng() - 1;
  }
  fout << "\n・" << method << "型" << opt.cc() * 100.0 << "%同時信頼区間(simultaneous confidence interval)\n";
  fout << "群      - 群        平均値の差    絶対精度        下限        上限\n";
  constexpr int LENline = 66;   // 区切り線の長さ
  putline(fout, LENline);
  list<gstat>::const_iterator gitr1 = gstatlst.begin();
  for (int i = 0; i < ng() - 1; ++i, ++gitr1)
  {
    string codei(gitr1->code());
    int npadi = codei.size() - strclen(codei, CHARwide);
    int ni = (gitr1->ystat()).n();
    double mi = (gitr1->ystat()).mean();
    list<gstat>::const_iterator gitr2 = gitr1;
    ++gitr2;
    while (gitr2 != gstatlst.end())
    {
      int nj = (gitr2->ystat()).n();
      double dm = (gitr2->ystat()).mean() - mi;
      double vd = (1.0 / (double)ni + 1.0 / (double)nj) * vr;
      double ci = calmcci(opt.m(), opt.cc(), ng(), df1, df2, dc, vd);
      int npad = gitr2->code().size() - strclen(gitr2->code(), CHARwide);
      fout << left << setw(8 + npad) << gitr2->code() << "- ";
      fout << setw(8 + npadi) << codei << right << setw(12) << dm;
      fout << setw(12) << ci << setw(12) << dm - ci << setw(12) << dm + ci;
      fout << "\n";
      ++gitr2;
    }
  }
  putline(fout, LENline);
}

void anovamclass::put_mdcid(ostream& fout, option& opt, const statv<double>& tstat, double vr) const
// 2群間の平均値の差のDunnett型信頼区間出力
{
  if (vr <= 0.0) return;        // 出力なし
  int df1 = 1, df2 = tstat.n() - ng(), dc = 1;
  fout << "\n・Dunnett型" << opt.cc() * 100.0 << "%同時信頼区間(simultaneous confidence interval)\n";
  fout << "群      - 群        平均値の差    絶対精度        下限        上限\n";
  constexpr int LENline = 66;   // 区切り線の長さ
  putline(fout, LENline);
  list<gstat>::const_iterator gitr = gstatlst.begin();
  string codei(gitr->code());
  int npadi = codei.size() - strclen(codei, CHARwide);
  int ni = (gitr->ystat()).n();
  double mi = (gitr->ystat()).mean();
  ++gitr;
  for (int j = 1; j < ng(); ++j, ++gitr)
  {
    int nj = (gitr->ystat()).n();
    double dm =  (gitr->ystat()).mean() - mi;
    double vd = (1.0 / ni + 1.0 / nj) * vr;
    double r = (ni == nj) ? 0.5 : sqrt(ni / (ni + nj) * nj / (ni + nj));
    double ci = calmcci(opt.m(), opt.cc(), ng(), df1, df2, dc, vd, r);
    int npad = gitr->code().size() - strclen(gitr->code(), CHARwide);
    fout << left << setw(8 + npad) << gitr->code() << "- ";
    fout << setw(8 + npadi) << codei << right << setw(12) << dm;
    fout << setw(12) << ci << setw(12) << dm - ci << setw(12) << dm + ci;
    fout << "\n";
  }
  putline(fout, LENline);
}

void anovamclass::plot_ngpid(ofstream& fout, int nd)
// Ngraph用データのタイトルとデータ識別情報を出力する
{
  if (id.ni() < 2) id.set(fname, 2);
  if (nd > 1)
  { fout << NGPdelm << "\n"; }
  else
  { fout << "-s1 -zx0," << gstatlst.size() + 1 << ",1 -dmark,78 -w80 -o400\n"; }
  fout << "%";  puttitle(fout, Title);
  fout << "%\n";
  id.put_title(fout, formngp);
  id.put_stra(fout, formngp);
  fout << "%群項目(要因A):" << id.item(0) << "\n";
  fout << "%集計項目     :" << id.item(1) << "\n";
}

void anovamclass::plot_ngp(ofstream& fout, option& opt) // Ngraph用データ出力
{
  fout << "%平均値±SE(-zx0," << gstatlst.size() + 1 << ",1 -derrorbar_y -w40 -o500)\n";
  int i = 0;
  for (auto& gstati : gstatlst)
  {
    double se = gstati.ystat().se();
    fout << i + 1 << " " << gstati.ystat().mean();
    fout << " " << se << " " << -se << "\n";
    ++i;
  }
  fout << NGPdelm2 << "\n";
  fout << "%平均値±SD\n";
  i = 0;
  for (auto& gstati : gstatlst)
  {
    double sd = gstati.ystat().sd();
    fout << i + 1 << " " << gstati.ystat().mean();
    fout << " " << sd << " " << -sd << endl;
    ++i;
  }
  fout << NGPdelm2 << "\n";
  fout << "%平均値±" << opt.cc() * 100.0 << "%信頼区間\n";
  i = 0;
  for (auto& gstati : gstatlst)
  {
    double t = ptot(1.0 - opt.cc(), gstati.ystat().n() - 1);
    double ci = t * gstati.ystat().se();
    double sd = gstati.ystat().sd();
    fout << i + 1 << " " << gstati.ystat().mean();
    fout << " " << ci << " " << -ci << endl;
    ++i;
  }
}

void anovamclass::getput(const char* c_fname, option& opt)
// データを読み込んで集計結果を出力する
{
  if (opt.check() != 0) throw option_error("実行時オプションエラー！");
  if (c_fname == nullptr) fname.erase();
  else fname = c_fname;
  istream* finp = nullptr;
  if (fname.empty())
  { finp = &cin;  fname = "stdin"; }
  else
  {
    finp = new ifstream(fname.c_str());
    if (!(*finp)) throw runtime_error(string("ファイルがありません→") + fname);
  }
  ofstream foutngp;
  if (opt.ngpname() != "")      // Ngraph用データ出力
  {
    if (ofopen(foutngp, (opt.ngpname()).c_str()) == false) return; // ファイルオープン失敗
  }
  string sline;                 // データ入力用
  datac cdata(5);               // データオブジェクト
  int nd = 1;                   // データの数
  while (getsline(sline, *finp) != EOF)
  {
    strchomp(sline);
    if (sline.empty())          // ブランク行
    {
      if (gstatlst.size() > 0)   // 群統計量がある時
      {
        put(cout, foutngp, opt, nd);
        id.clear();  gstatlst.clear();  ++nd;
      }
    }
    else if (sline[0] == CHRtags) // データ識別情報
    {
      if (gstatlst.size() > 0)   // 群統計量がある時
      {
        put(cout, foutngp, opt, nd);
        id.clear();  gstatlst.clear();  ++nd;
      }
      if (id.append(sline) < 0) throw runtime_error(string("データ識別情報が間違っています→") + sline);
      if (id.ni() > 2) throw runtime_error(string("項目が多すぎます→") + sline);
    }
    else if (cdata.getdata(sline) > 0) // 統計量データ読み込み
    {
      gstat tgstat;
      get_stat(tgstat, cdata);
      gstatlst.push_back(tgstat);
    }
  }
  if (finp != &cin) delete finp;
  put(cout, foutngp, opt, nd);
}

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