// anova1.C Version 7.45  Copyleft (c) 杉本典夫
// Time-stamp: <2026-04-08 13:47:51 tac>
//
// 対応の無い多群間で一元配置分散分析を行う。
// 2群の場合は2標本t検定も出力する。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法  :g++ -O -s -W anova1 anova1.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 =            // プログラム使用法表示用
{
  "anova1  Version 7.45         Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:anova1  [<オプション> …]  [<データファイル> …]",
  "機能:最初の項目を群項目(要因A)とし、次の項目を集計項目とした一元配置分散",
  "     分析(2標本t検定含む)を行います。",
  "<オプション>",
  "  -al=<有意水準>:検出力分析の有意水準αを指定します。(0<α<1,省略時:0.05)",
  "  -cc[=<信頼係数>]:信頼係数を指定します。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -dm=<検出差>:検出力分析の検出差を指定します。(省略時:標本平均値の差)",
  "  -f=<ファイル名>:指定したファイルからオプションを読み込みます。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示します。(-h:-h=man)",
  "  -ngp={bar|box|ci|data|msd|mse|sd[<係数>]|sr[<係数>]}[,…]",
  "      :Ngraph用データを指定します。",
  "    bar:バーチャート  box:箱ヒゲプロット  ci:平均値±CI  data:個々の生データ",
  "    msd:平均値±SD  mse:平均値±SE  sd[<係数>]:平均値±<係数>×標準偏差の折れ線",
  "    sr[<係数>]:平均値±<係数>×残差  オプション省略時:全て指定",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データを出力するファイル名を指定する。",
  "    オプション省略時:データ出力せず  -ngpname:-ngpname=anova1_ngp.dat",
  "  -res=[-][{def|all|at|ft|m[{a|d}][{b|d|f|s|t}]|mci|pa|r|stat",
  "      |tt|vt}][,…]:解析結果の出力項目を指定します。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    at:分散分析  ft:等分散性の両側F検定",
  "    m[{a|d}][{b|d|f|s|t}]:多重比較(ma:全ての2群比較",
  "    md:群1と他の群の比較)  mci:平均値の信頼区間  pa:検出力分析  r:相関比",
  "    stat:基礎統計量  tt:2標本t検定  vt:等分散性の検定",
  "    オプション省略時:データに合わせて自動判定",
  "  -po=<検出力>:検出力分析の検出力(1-β)を指定します。(0＜1-β＜1,省略時:0.8)"
};

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

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

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_al;                  // 有意水準
  double v_cc;                  // 信頼係数
  double v_dm;                  // 検出差
  mctype v_ma;                  // 全ての2群ごとの多重比較法
  mctype v_md;                  // 対照群と他の群の多重比較法
  double v_po;                  // 検出力
  double v_tv;                  // -ngp=s[dr]<係数>
  string v_ngpname;             // Ngraphファイル名
  lstr ngplst;                  // Ngraph出力リスト
  mapsb ngpmap;                 // Ngraph出力の有無マップ
  lstr reslst;                  // 結果出力リスト
  mapsb resmap;                 // 結果出力の有無マップ
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  double al() const             // 有意水準を返す
  { return v_al; }
  double cc() const             // 信頼係数を返す
  { return v_cc; }
  double dm() const             // 検出差を返す
  { return v_dm; }
  mctype ma() const             // 全ての2群ごとの多重比較法を返す
  { return v_ma; }
  mctype md() const             // 対照群と他の群の多重比較法を返す
  { return v_md; }
  double po() const             // 検出力を返す
  { return v_po; }
  double tv() const             // SDの係数を返す
  { return v_tv; }
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  bool ngp(const string& code); // Ngraph出力フラッグを返す
  bool res(const string& code); // 解析結果出力フラッグを返す
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check(size_t ni);         // オプションのチェックとセット
  int set_ngp();                // -ngpオプションのセット
  int set_res();                // -resオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_al(ALdef), v_cc(CCdef), v_dm(BLANKf), v_ma(mctuk), v_md(mcdun), v_po(POdef), v_tv(BLANKf), v_ngpname(), ngplst(), ngpmap(), reslst(), resmap()
// コンストラクタ
{
  if (c_pname == nullptr) return; // 初期設定せず
  v_pname = getname(c_pname);     // パス名なしのプログラム名取得
  v_path = c_pname;
  strsize_t idx = v_path.find(v_pname);
  if (idx != string::npos) v_path.erase(idx); // プログラム名削除
}

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

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

void option::set(const string& optstr, char devc)
// オプション指定によってオプションを設定する
{
  string opt[2];                // オプションのopt[0]:項目,opt[1]:内容
  if (getopt(opt, optstr, devc) < 1) return; // オプション指定なし
  int err = 0;
  if (opt[0] == "al")           // 有意水準の指定
  {
    v_al = opt[1].empty() ? ALdef : strtodf(opt[1]);
    if (v_al <= 0.0 || 1.0 <= v_al) err = badopt(pname(), optstr, USAGE, 0.0, 1.0);
  }
  else 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] == "dm")      // 検出差の指定
  {
    v_dm = opt[1].empty() ? BLANKf : strtodf(opt[1]);
    if (v_dm < 0.0) err = badopt(pname(), optstr, USAGE, 0.0, FDmax);
  }
  else if (opt[0] == "f")       // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else if (opt[0] == "ngp")     // Ngraph出力項目指定
  {
    if (strsplit(ngplst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "ngpname")     // Ngraphファイル名指定
  {
    if (opt[1].empty()) v_ngpname = NGPNAMEdef;
    else v_ngpname = opt[1];
    if (v_ngpname.size() < 3) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "po")      // 検出力の指定
  {
    v_po = opt[1].empty() ? POdef : strtodf(opt[1]);
    if (v_po <= 0.0 || 1.0 <= v_po) err = badopt(pname(), optstr, USAGE, 0.0, 1.0);
  }
  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_ngp() < 0) return -2; // -ngpオプションのセット
  if (set_res() < 0) return -3; // -resオプションのセット
  return 0;                     // エラーなし
}

int option::set_ngp()           // -ngpオプションのセット
{
  ngpmap["bar"] = ngpmap["box"] = ngpmap["ci"] = ngpmap["data"] = ngpmap["msd"] = ngpmap["mse"] = ngpmap["sd"] = ngpmap["sr"] = true;
  if (ngplst.size() > 0)
  {
    if (v_ngpname == "") return disperr(pname(), "-ngpnameオプションの指定が必要！");
    ngpmap["bar"] = ngpmap["box"] = ngpmap["ci"] = ngpmap["data"] = ngpmap["msd"] = ngpmap["mse"] = ngpmap["sd"] = ngpmap["sr"] = false;
    for (const string& ngpstr : ngplst)
    {
      mapsb_citr ngppos = ngpmap.find(ngpstr);
      if (ngpstr[0] == 's')
      {
        switch (ngpstr[1])
        {
        case 'd' : ngpmap["sd"] = true;
          break;
        case 'r' : ngpmap["sr"] = true;
          break;
        default : return disperr(pname(), "-ngp=" + ngpstr + "オプション指定が不適切です！");
        }
        v_tv = (ngpstr.size() > 2) ? strtodf(ngpstr.substr(2)) : BLANKf;
        if (v_tv < 0.0) return disperr(pname(), "-ngp" + ngpstr + "オプション指定が不適切です！");
      }
      else 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["ft"] = resmap["stat"] = resmap["tt"] = true;
  resmap["at"] = resmap["ma"] = resmap["md"] = resmap["mci"] = resmap["pa"] = 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["ma"] = resmap["md"] = resmap["mci"] = resmap["pa"] = resmap["r"] = resmap["stat"] = resmap["tt"] = resmap["vt"] = true;
      }
      else if (resstr == "m")      // -res=mと指定された時
      { resmap["ma"] = resmap["md"] = true; }
      else if (resstr[0] == 'm' && resstr[1] == 'a' && resstr.size() > 2)
      {
        switch (resstr[2])
        {
        case 'b' : v_ma = mcbon;
          break;
        case 'f' : v_ma = mcfis;
          break;
        case 's' : v_ma = mcsch;
          break;
        case 't' : v_ma = mctuk;
          break;
        default : return disperr(pname(), "-res=" + resstr + "オプション指定が不適切です！");
        }
        resmap["ma"] = true;
      }
      else if (resstr[0] == 'm' && resstr[1] == 'd' && resstr.size() > 2)
      {
        switch (resstr[2])
        {
        case 'b' : v_md = mcbon;
          break;
        case 'd' : v_md = mcdun;
          break;
        case 'f' : v_md = mcfis;
          break;
        case 's' : v_md = mcsch;
          break;
        default : return disperr(pname(), "-res=" + resstr + "オプション指定が不適切です！");
        }
        resmap["md"] = true;
      }
      else if (respos == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切です！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

class varclass                  // 変数クラス
{
private:
  map<string, listdbl> gfdmap;    // 群データマップ
  listdbl dumfdlst;               // ダミーの群データリスト
  vector<string> gcodevec;      // 群コードベクトル
  vector<statv<double>> gstatvec; // 基礎統計量ベクトル
  double v_sa;                  // 要因A平方和
  double v_sr;                  // 残差平方和
public:
  varclass() : gfdmap(), dumfdlst(), gcodevec(), gstatvec(), v_sa(0.0), v_sr(0.0)
  // デフォールトコンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  void cal(ldata& datalst);     // データを集計して変数をセットする
  size_t ng() const             // 群数を返す
  { return gfdmap.size(); }
  const string& gcode(size_t idx) const // 群コードを返す
  { return (idx < ng()) ? gcodevec[idx] : BLANKs; }
  listdbl& fdlst(size_t idx);   // 群データリストを返す
  int n(size_t idx) const       // 群の例数を返す
  { return (idx <= ng()) ? gstatvec[idx].n() : gstatvec[ng()].n(); }
  double mean(size_t idx) const // 群の平均値を返す
  { return (idx <= ng()) ? gstatvec[idx].mean() : gstatvec[ng()].mean(); }
  double ss(size_t idx) const   // 群の平方和を返す
  { return (idx <= ng()) ? gstatvec[idx].ss() : gstatvec[ng()].ss(); }
  double vr(size_t idx) const   // 群の不偏分散を返す
  { return (idx <= ng()) ? gstatvec[idx].vr() : gstatvec[ng()].vr(); }
  double sd(size_t idx) const   // 群の標準偏差を返す
  { return (idx <= ng()) ? gstatvec[idx].sd() : gstatvec[ng()].sd(); }
  double se(size_t idx) const   // 群の標準誤差を返す
  { return (idx <= ng()) ? gstatvec[idx].se() : gstatvec[ng()].se(); }
  int n() const                 // 全体の例数を返す
  { return gstatvec[ng()].n(); }
  double mean() const           // 全体の平均値を返す
  { return gstatvec[ng()].mean(); }
  double ss() const             // 全体の平方和を返す
  { return gstatvec[ng()].ss(); }
  double vr() const             // 全体の不偏分散を返す
  { return gstatvec[ng()].vr(); }
  double sd() const             // 全体の標準偏差を返す
  { return gstatvec[ng()].sd(); }
  double se() const             // 全体の標準誤差を返す
  { return gstatvec[ng()].se(); }
  double sa() const             // 要因A平方和を返す
  { return v_sa; }
  int dfa() const               // 要因A自由度を返す
  { return ng() - 1; }
  double va() const             // 要因A分散を返す
  { return (dfa() > 0) ? sa() / (double)(dfa()) : 0.0; }
  double sr() const             // 残差平方和を返す
  { return v_sr; }
  int dfr() const               // 残差自由度を返す
  { return dft() - dfa(); }
  double vrr() const            // 残差分散を返す
  { return (dfr() > 0) ? sr() / (double)(dfr()) : 0.0; }
  double st() const             // 全体平方和を返す
  { return gstatvec[ng()].ss(); }
  int dft() const               // 全体自由度を返す
  { return gstatvec[ng()].n() - 1; }
};

void varclass::cal(ldata& datalst) // データを集計して変数をセットする
{
  for (auto& cdatai : datalst)
  {
    if (cdatai.sdata(0) != BLANKs && cdatai.fdata(1) != BLANKf)
    {
      map<string, listdbl>::iterator gitr = gfdmap.find(cdatai.sdata(0));
      if (gitr == gfdmap.end()) // 新規群コード
      {
        listdbl tfdlst(1, cdatai.fdata(1));
        gfdmap.insert(map<string, listdbl>::value_type(cdatai.sdata(0), tfdlst));
      }
      else
      { gitr->second.push_back(cdatai.fdata(1)); }
    }
  }
  size_t ng = gfdmap.size();    // 群数
  if (ng < 1) return;           // データなし
  gcodevec.resize(ng);  gstatvec.resize(ng + 1);
  size_t i = 0;
  for (auto& gfdi : gfdmap)
  {
    gcodevec[i] = gfdi.first;  // 群コードをべクルトに格納する
    for (double fdata : gfdi.second) (gstatvec[i]).sigma(fdata); // データ累積
    ++i;
  }
  v_sa = 0.0;
  for (size_t i = 0; i < ng; ++i)
  {
    double m = gstatvec[i].mean();
    v_sa += (double)(gstatvec[i].n()) * m * m;
    gstatvec[ng] += gstatvec[i];
  }
  double m = gstatvec[ng].mean();
  v_sa -= (double)(gstatvec[ng].n()) * m * m;
  if (dfa() < 1 || v_sa < FDeps) v_sa = 0.0;
  v_sr = st() - v_sa;
  if (dfr() < 1 || v_sr < FDeps) v_sr = 0.0;
}

listdbl& varclass::fdlst(size_t idx) // 群のデータリストを返す
{
  if (idx >= gcodevec.size()) return dumfdlst; // ダミー群データを返す
  map<string, listdbl>::iterator gpos = gfdmap.find(gcodevec[idx]);
  if (gpos == gfdmap.end()) return dumfdlst; // ダミー群データを返す
  return gfdmap[gcodevec[idx]];              // 群データリストを返す
}

class anova1class               // anova1メインクラス
{
private:
  string fname;                 // ファイル名
  void put_resid(ostream& fout, idd& id) const;
  // 分析結果のタイトルとデータ識別情報を出力する
  void put_res2g(ostream& fout, option& opt, const varclass& var) const;
  // 2群比較の解析結果を出力する
  void put_resmg(ostream& fout, option& opt, const varclass& var) const;
  // 3群以上の時の解析結果を出力する
  void put_stat(ostream& fout, const varclass& var) const;
  // 基礎統計量を出力する
  void put_mci(ostream& fout, option& opt, const varclass& var) const;
  // 群別平均値の信頼区間出力
  void put_rratio(ostream& fout, option& opt, const varclass& var) const;
  // 相関比の出力
  void put_ttest(ostream& fout, option& opt, const varclass& var) const;
  // 2標本t検定とWelchの検定
  void put_ftest(ostream& fout, option& opt, double f, int df1, int df2, ldbl pr) const;
  // 分散比Fの両側検定と推定結果を出力して両側p値を返す
  void put_pa(ostream& fout, option& opt, const varclass& var) const;
  // 検出力分析の検定
  void put_vtest(ostream& fout, const varclass& var) const;
  // 分散の均一性の検定結果出力
  void put_mcop(ostream& fout, const varclass& var, mctype mc, bool ma) const;
  // 2群間の多重比較出力
  void put_mdci(ostream& fout, option& opt, const varclass& var, mctype mc, bool ma) const;
  // 2群間の平均値の差の信頼区間出力
  void put_mpa(ostream& fout, option& opt, const varclass& var, mctype mc, bool ma) const;
  // 2群間の平均値の差の検出力分析結果出力
  void put_dns(ostream& fout, option& opt, const varclass& var) const;
  // 分散分析の結果をデータファイル形式で出力する
  void plot_ngp(option& opt, idd& id, const ldata& datalst, varclass& var) const;
  // Ngraph用データ出力
  void plot_ngpid(ofstream& fout, idd& id) const;
  // Ngraph用データのタイトルとデータ識別情報を出力する
  void plot_bar(ofstream& fout, const varclass& var) const;
  // Ngraph用にバーチャートデータを出力する
  void plot_box(ofstream& fout, varclass& var) const;
  // Ngraph用に箱ヒゲプロットデータを出力する
  void plot_ci(ofstream& fout, option& opt, const varclass& var) const;
  // Ngraph用に平均±信頼区間データを出力する
  void plot_data(ofstream& fout, varclass& var) const;
  // Ngraph用に個々のデータを出力する
  void plot_msd(ofstream& fout, const varclass& var) const;
  // Ngraph用に平均±SDを出力する
  void plot_mse(ofstream& fout, const varclass& var) const;
  // Ngraph用に平均±SEを出力する
  void plot_sdt(ofstream& fout, option& opt, const varclass& var) const;
  // Ngraph用に平均±係数×SDを出力する
  void plot_srt(ofstream& fout, option& opt, const varclass& var) const;
  // Ngraph用に平均±係数×残差SDを出力する
public:
  anova1class() : fname()       // コンストラクタ
  {}
  ~anova1class()                // デストラクタ
  {}
  int get(const char* c_fname, option& opt, idd& id, ldata& datalst);
  // データを読み込む
  void put(option& opt, idd& id, ldata& datalst) const; // 解析結果を出力する
};

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

void anova1class::put_res2g(ostream& fout, option& opt, const varclass& var) const
// 2群比較の解析結果を出力する
{
  if (opt.res("r"))             // 相関比出力
  {
    fout << "\n";
    put_rratio(fout, opt, var);
  }
  if (opt.res("at"))           // 2群の時は-res=atの時だけ分散分析出力
  {
    if (opt.res("vt"))          // 等分散性の検定出力
    {
      fout << "\n";
      put_vtest(fout, var);
    }
    fout << "\n・一元配置分散分析(one-way layout analysis of variance)\n";
    putatab(fout, "群(要因A)", var.sa(), var.dfa(), var.sr(), var.dfr(), var.st(), var.dft());
  }
  if (opt.res("tt"))            // 2標本t検定出力
  {
    fout << "\n";
    put_ttest(fout, opt, var);
  }
  if (opt.res("pa"))            // 検出力分析出力
  {
    fout << "\n";
    put_pa(fout, opt, var);
  }
}

void anova1class::put_resmg(ostream& fout, option& opt, const varclass& var) const
// 3群以上の時の解析結果を出力する
{
  if (opt.res("def") || opt.res("r")) // 省略時または-res=rの時に相関比出力
  {
    fout << "\n";
    put_rratio(fout, opt, var);
  }
  if (opt.res("vt"))            // 等分散性の検定出力
  {
    fout << "\n";
    put_vtest(fout, var);
  }
  if (opt.res("def") || opt.res("at")) // 省略時または-res=atの時に分散分析出力
  {
    fout << "\n";
    fout << "・一元配置分散分析(one-way layout analysis of variance)\n";
    putatab(fout, "群(要因A)", var.sa(), var.dfa(), var.sr(), var.dfr(), var.st(), var.dft());
  }
  if (opt.res("ma"))            // 全ての2群の多重比較出力
  {
    fout << "\n";
    put_mcop(fout, var, opt.ma(), true);
    fout << "\n";
    put_mdci(fout, opt, var, opt.ma(), true);
    if (opt.res("pa"))          // 検出力分析出力
    {
      fout << "\n";
      put_mpa(fout, opt, var, opt.ma(), true);
    }
  }
  if (opt.res("md"))            // 群1と他群の多重比較出力
  {
    fout << "\n";
    put_mcop(fout, var, opt.md(), false);
    fout << "\n";
    put_mdci(fout, opt, var, opt.md(), false);
    if (opt.res("pa"))          // 検出力分析出力
    {
      fout << "\n";
      put_mpa(fout, opt, var, opt.md(), false);
    }
  }
}

void anova1class::put_stat(ostream& fout, const varclass& var) const
// 基礎統計量を出力する
{
  constexpr int LENline = 77;   // 区切り線の長さ
  constexpr int LENhed = 9;     // ヘッダの長さ
  fout << "群       :群別基礎統計量\n";
  putline(fout, LENline);
  for (size_t i = 0; i < var.ng(); ++i)
  {
    string gcodei(var.gcode(i));
    int npad = gcodei.size() - strclen(gcodei, CHARwide);
    fout << left << setw(LENhed + npad) << gcodei << right << ":";
    putmsd(fout, var.n(i), var.mean(i), var.vr(i));
  }
  putline(fout, LENline);
  string tstr("全体");
  int npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(LENhed + npad) << tstr << right << ":";
  putmsd(fout, var.n(), var.mean(), var.vr());
}

void anova1class::put_mci(ostream& fout, option& opt, const varclass& var) const
// 群別平均値の信頼区間を出力する
{
  constexpr int LENline = 57;   // 区切り線の長さ
  constexpr int LENhed = 9;     // ヘッダの長さ
  fout << "・群別平均値の" << opt.cc() * 100.0 << "%信頼区間(confidence interval)\n";
  fout << "群             平均値    絶対精度        下限        上限\n";
  putline(fout, LENline);
  for (size_t i = 0; i < var.ng(); ++i)
  {
    string gcodei(var.gcode(i));
    int npad = gcodei.size() - strclen(gcodei, CHARwide);
    int df = var.n(i) - 1;
    double me = var.mean(i);
    double t = (df > 0) ? (double)(ptot((ldbl)(1.0 - opt.cc()), df)) : 0.0;
    double ci = t * var.se(i);
    fout << left << setw(LENhed + npad) << gcodei << right;
    fout  << setw(12) << me << setw(12) << ci << setw(12) << me - ci << setw(12) << me + ci << "\n";
  }
  putline(fout, LENline);
  string tstr("全体");
  int npad = tstr.size() - strclen(tstr, CHARwide);
  int df = var.n() - 1;
  double me = var.mean();
  double t = (df > 0) ? (double)(ptot((ldbl)(1.0 - opt.cc()), df)) : 0.0;
  double ci = t * var.se();
  fout << left << setw(LENhed + npad) << tstr << right;
  fout  << setw(12) << me << setw(12) << ci << setw(12) << me - ci << setw(12) << me + ci << "\n";
}

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

void anova1class::put_ttest(ostream& fout, option& opt, const varclass& var) const
// 2標本t検定とWelchの検定結果を出力する
{
  int n1 = var.n(0), n2 = var.n(1);
  double ss1 = var.ss(0), ss2 = var.ss(1);
  double v1 = var.vr(0), v2 = var.vr(1);
  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 = 1.0L;
  if (df1 > 0 && df2 > 0 && sgn(ss1 * ss2) > 0) pr = ftop((ldbl)f, df1, df2);
  pr = (pr < 0.5L) ? 2.0L * pr : 1.0L; // 両側有意確率にする
  if (opt.res("ft"))                   // 等分散性の両側F検定出力
  {
    put_ftest(fout, opt, f, df1, df2, pr);
    fout << "\n";
  }
  fout << "・対応のないt検定(2標本t検定、「→」の付いた方を採用して下さい)\n";
  string sel1 = "  →", sel2 = "    ";
  if (pr <= 0.05L) { sel1 = "    ";  sel2 = "  →"; }
  double d = var.mean(1) - var.mean(0);
  int df = n1 + n2 - 2;
  double vd = (ss1 + ss2) / (double)df * (double)(n1 + n2) / (double)n1 / (double)n2;
  double t = (vd > 0.0) ? d / sqrt(vd) : 0.0;
  pr = ttop((ldbl)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 = 0.0;
  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(ss1 * ss2) <= 0) return;
  v1 /= n1;  v2 /= n2;
  vd = v1 + v2;
  df = (int)(vd * vd / (v1 * v1 / (n1 - 1) + v2 * v2 / (n2 -1)));
  t = (vd > 0.0) ? d / sqrt(vd) : 0.0;
  pr = ttop((ldbl)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";
}

void anova1class::put_ftest(ostream& fout, option& opt, double f, int df1, int df2, ldbl pr) const
// 分散比Fの両側検定と推定結果を出力する
{
  fout << "・等分散性の両側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";
}

void anova1class::put_pa(ostream& fout, option& opt, const varclass& var) const
// 検出力分析結果を出力する
{
  int n1 = var.n(0), n2 = var.n(1);
  double ss1 = var.ss(0), ss2 = var.ss(1);
  if (n1 < 2 || n2 < 2 || (ss1 <= 0.0 && ss2 <= 0.0)) return; // 出力なし
  int nt = n1 + n2;
  int df = nt - 2;
  double sd = sqrt((ss1 + ss2) / (double)df);
  double dm = (opt.dm() != BLANKf) ? opt.dm() : fabs(var.mean(1) - var.mean(0));
  if (dm < FDeps) dm = sd;      // 平均値の差が0の時はSDにする
  double d = dm / sd;
  double lamda = d / sqrt(1.0 / (double)n1 + 1.0 / (double)n2);
  double ta = (double)(ptot((ldbl)(opt.al()), df));
  double t2b = lamda - ta;
  double b = (t2b > 0.0) ? (double)(ttop((ldbl)t2b, df)) * 0.5 : 1.0 - (double)(ttop((ldbl)(-t2b), df)) * 0.5;
  double b0 = (double)(ncttolp((ldbl)ta, df, (ldbl)lamda));
  fout << "・検出力分析:有意水準=" << opt.al() << "  平均値の差の標準偏差=" << sd << "\n";
  fout << "  検出差" << dm << "(効果量=" << d << ")の時の検出力=" << 1.0 - b0 << "(β=" << b0 << ")\n";
  t2b = (double)(ptot((ldbl)(2.0 * (1.0 - opt.po())), df));
  double s = (double)n2 / (double)n1;
  double nc = ta * ta / 4.0 / (1.0 - 1.0 / (double)nt);
  if (n1 <= nc) nc = 0.0;
  d = (ta + t2b) * sqrt((1.0 + s) / s / ((double)n1 - nc));
  dm = d * sd;
  fout << "  検出力=" << opt.po() << "(β=" << 1.0 - opt.po() << ")の時の検出差=" << dm << "(効果量=" << d << ")\n";
  ta = (double)(ptoz((ldbl)(opt.al())));
  t2b = (double)(ptoz((ldbl)(2.0 * (1.0 - opt.po()))));
  if (opt.po() < 0.5) t2b = (double)(-ptoz((ldbl)(2.0 * opt.po())));
  dm = (opt.dm() != BLANKf) ? opt.dm() : fabs(var.mean(1) - var.mean(0));
  if (dm < FDeps) dm = sd;      // 平均値の差が0の時はSDにする
  d = dm / sd;
  double nsqrt = (ta + t2b) / d;
  double n0 = 2.0 * nsqrt * nsqrt;
  n1 = (int)(n0 + 0.99);
  n1 = (n1 < 2) ? n1 + 1 : (int)(n0 + ta * ta / 4.0 / (1.0 - 1.0 / (double)n1) + 0.99);
  fout << "  検出力=" << opt.po() << "・検出差=" << dm << "の時の1群必要例数=" << n1 << "(2群合計必要例数=" << n1 * 2 << ")" << "\n";
}

void anova1class::put_vtest(ostream& fout, const varclass& var) const
// 分散の均一性の検定結果を出力する
{
  int ng = var.ng();
  int d = -ng;
  double w = 0.0, ss = 0.0, lnss = 0.0;
  for (int i = 0; i < ng; ++i)
  {
    int df = var.n(i) - 1;
    d += var.n(i);
    w += 1.0 / (double)df;
    ss += (double)df * var.vr(i);
    lnss += (double)df * log(var.vr(i));
  }
  int df = ng - 1;
  double m = (d > 0 && ss > 0.0) ? d * log(ss / 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((ldbl)csqr, df);
  fout << "・等分散性の検定(Bartlettの検定):χ2乗=" << csqr;
  fout << "  自由度=" << df << "  p=" << pr << pind(pr) << "\n";
}

void anova1class::put_mcop(ostream& fout, const varclass& var, mctype mc, bool ma) const
// 2群間の多重比較結果を出力する
{
  constexpr int LENline = 75;   // 区切り線の長さ
  constexpr int LENhed = 9;     // ヘッダの長さ
  double vr = var.vrr();        // 残差分散
  if (vr <= 0.0) return;        // 出力なし
  int ng = var.ng();
  int df1 = 1, df2 = var.dfr(), dc = 1;
  string method = "Tukey";
  switch (mc)
  {
  case mcbon : method = "Bonferroni";
    dc = ma ? ng * (ng - 1) / 2 : ng - 1;
    break;
  case mcdun : method = "Dunnett";
    break;
  case mcfis : method = "Fisher";
    break;
  case mcsch : method = "Scheffe";
    df1 = ng - 1;
  }
  fout << "・群(要因A)の" << method << "型多重比較(" << method << " type multiple comparison)\n";
  if (mc == mctuk)
  {
    fout << "群       - 群                q値        群数      自由度   有意確率p値\n";
  }
  else if (mc == mcdun)
  {
    fout << "群       - 群                d値        群数      自由度   有意確率p値\n";
  }
  else
  {
    fout << "群       - 群                F値   第1自由度   第2自由度   有意確率p値\n";
  }
  putline(fout, LENline);
  int iend = ma ? ng - 1 : 1;
  for (int i = 0; i < iend; ++i)
  {
    string gcodei(var.gcode(i));
    int npadi = gcodei.size() - strclen(gcodei, CHARwide);
    double ni = (double)(var.n(i));
    double mi = var.mean(i);
    for (int j = i + 1; j < ng; ++j)
    {
      double nj = (double)(var.n(j));
      double dm = var.mean(j) - 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));
      ldbl pr = calmcpr(mc, f, ng, df1, df2, dc, r);
      string gcodej(var.gcode(j));
      int npad = gcodej.size() - strclen(gcodej, CHARwide);
      fout << left << setw(LENhed + npad) << gcodej;
      fout << "- ";
      fout << setw(LENhed + npadi) << gcodei << right << setw(12) << f;
      if (mc == mctuk || mc == mcdun) fout << setw(12) << ng;
      else fout << setw(12) << df1;
      fout << setw(12) << df2 << setw(14) << pr << pind(pr) << "\n";
    }
  }
  putline(fout, LENline);
}

void anova1class::put_mdci(ostream& fout, option& opt, const varclass& var, mctype mc, bool ma) const
// 2群間の平均値の差の信頼区間を出力する
{
  constexpr int LENline = 69;   // 区切り線の長さ
  constexpr int LENhed = 9;     // ヘッダの長さ
  double vr = var.vrr();        // 残差分散
  if (vr <= 0.0) return;        // 出力なし
  int ng = var.ng();
  int df1 = 1, df2 = var.dfr(), dc = 1;
  string method = "Tukey";
  switch (mc)
  {
  case mcbon : method = "Bonferroni";
    dc = ma ? ng * (ng - 1) / 2 : ng - 1;
    break;
  case mcdun : method = "Dunnett";
    break;
  case mcfis : method = "Fisher";
    break;
  case mcsch : method = "Scheffe";
    df1 = ng - 1;
  }
  fout << "・" << method << "型" << opt.cc() * 100.0 << "%同時信頼区間(simultaneous confidence interval)\n";
  fout << "群       - 群         平均値の差    絶対精度        下限        上限\n";
  putline(fout, LENline);
  int iend = ma ? ng - 1 : 1;
  for (int i = 0; i < iend; ++i)
  {
    string gcodei(var.gcode(i));
    int npadi = gcodei.size() - strclen(gcodei, CHARwide);
    int ni = var.n(i);
    double mi = var.mean(i);
    for (int j = i + 1; j < ng; ++j)
    {
      int nj = var.n(j);
      double dm = var.mean(j) - 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(mc, opt.cc(), ng, df1, df2, dc, vd, r);
      string gcodej(var.gcode(j));
      int npad = gcodej.size() - strclen(gcodej, CHARwide);
      fout << left << setw(LENhed + npad) << gcodej;
      fout << "- ";
      fout << setw(LENhed + npadi) << gcodei << right << setw(12) << dm;
      fout << setw(12) << ci << setw(12) << dm - ci << setw(12) << dm + ci;
      fout << "\n";
    }
  }
  putline(fout, LENline);
}

void anova1class::put_mpa(ostream& fout, option& opt, const varclass& var, mctype mc, bool ma) const
// 2群間の平均値の差の検出力分析結果出力
{
  constexpr int LENline = 69;   // 区切り線の長さ
  constexpr int LENhed = 9;     // ヘッダの長さ
  double vr = var.vrr();        // 残差分散
  if (vr <= 0.0) return;        // 出力なし
  double sd = sqrt(vr);
  int ng = var.ng();
  int df1 = 1, df2 = var.dfr(), dc = 1;
  double poc = opt.po(), ac = opt.al();
  double bc = 1.0 - poc, bcc = 1.0;
  double tac = (double)(ptot((ldbl)ac, df2)), zac = (double)(ptoz((ldbl)ac));
  double t2bc = (poc > 0.5) ? (double)(ptot((ldbl)(2.0 * bc), df2)) : (ldbl)(-ptot((ldbl)(2.0 * poc), df2));
double z2bc = (poc > 0.5) ? (double)(ptoz((ldbl)(2.0 * bc))) : (double)(-ptoz((ldbl)(2.0 * poc)));
  string method = "Fisher";
  switch (mc)
  {
  case mcdun :
  case mctuk :
  case mcbon : method = "Bonferroni";
    dc = ma ? ng * (ng - 1) / 2 : ng - 1;
    ac = opt.al() / dc;
    tac = (double)(ptot((ldbl)ac, df2));
    zac = (double)(ptoz((ldbl)ac));
    break;
  case mcsch : method = "Scheffe";
    df1 = ng - 1;
    tac = (double)(sqrt(ptof((ldbl)ac, df1, df2) * (ldbl)df1));
    zac = (double)(sqrt(ptoc((ldbl)ac, df1)));
  }
  fout << "・検出力分析:設定有意水準=" << opt.al() << "  設定検出力=" << opt.po() << "  " << method << "型多重比較\n";
  fout << "群       - 群         設定検出差  調整検出力  調整検出差 1群必要例数\n";
  putline(fout, LENline);
  int iend = ma ? ng - 1 : 1;
  for (int i = 0; i < iend; ++i)
  {
    string gcodei(var.gcode(i));
    int npadi = gcodei.size() - strclen(gcodei, CHARwide);
    double ni = (double)(var.n(i));
    double mi = var.mean(i);
    for (int j = i + 1; j < ng; ++j)
    {
      double nj = (double)(var.n(j));
      double dm = (opt.dm() != BLANKf) ? opt.dm() : fabs(var.mean(j) - mi);
      if (dm < FDeps) dm = sd; // 平均値の差が0の時はSDにする
      double lamda = dm / sqrt((1.0 / ni + 1.0 / nj) * vr);
      double b = ncttolp(tac, df2, lamda);
      double po = 1.0 - b;
      bcc *= b;
      double s = nj / ni;
      double nc = tac * tac / 4.0 / (1.0 - 1.0 / (ni + nj));
      if (ni <= nc) nc = 0.0;
      double d = (tac + t2bc) * sqrt((1.0 + s) / s / (ni - nc));
      double dmc = d * sd;
      double nsqrt = (zac + z2bc) * sd / dm;
      double n0 = 2.0 * nsqrt * nsqrt;
      int n = (int)(n0 + 0.99);
      n = (n < 2) ? n + 1 : (int)(n0 + zac * zac / 4.0 / (1.0 - 1.0 / (double)n) + 0.99);
      string gcodej(var.gcode(j));
      int npad = gcodej.size() - strclen(gcodej, CHARwide);
      fout << left << setw(LENhed + npad) << gcodej;
      fout << "- ";
      fout << setw(LENhed + npadi) << gcodei << right << setw(12) << dm;
      fout << setw(12) << po << setw(12) << dmc << setw(12) << n;
      fout << "\n";
    }
  }
  putline(fout, LENline);
  if (mc != mcfis) fout << "全体の調整検出力=" << 1.0 - bcc << "(β=" << bcc << ")" << "\n";
}

void anova1class::plot_ngp(option& opt, idd& id, const ldata& datalst, varclass& var) const
// Ngraph用データ出力
{
  ofstream fout;
  if (ofopen(fout, (opt.ngpname()).c_str()) == false) return; // ファイルオープン失敗
  bool delml = false;
  if (opt.ngp("bar"))           // バーチャート
  {
    fout << "-s1 -zx0," << var.ng() + 1 << ",1 -drectangle_solid_fill -x1 -y3\n";
    plot_ngpid(fout, id);
    plot_bar(fout, var);
    delml = true;
  }
  if (opt.ngp("box"))           // 箱ひげプロット
  {
    if (delml)                   // 区切り行を出力する
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -zx0," << var.ng() + 1 << ",1 -drectangle -x1 -y3 -w80\n";
      delml = true;
    }
    plot_ngpid(fout, id);
    plot_box(fout, var);
  }
  if (opt.ngp("ci"))            // 平均値±CI
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -zx0," << var.ng() + 1 << ",1 -dmark,78 -w80 -o400\n";
      delml = true;
    }
    plot_ngpid(fout, id);
    plot_ci(fout, opt, var);
    delml = true;
  }
  if (opt.ngp("data"))          // 個々のデータ
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -zx0," << var.ng() + 1 << ",1 -dmark,00 -o300\n";
      delml = true;
    }
    plot_ngpid(fout, id);
    plot_data(fout, var);
  }
  if (opt.ngp("msd"))           // 平均値±SD
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -zx0," << var.ng() + 1 << ",1 -dline -w80\n";
      fout << "%エラーバーを描く時:-derrorbar_y -w80 -o500\n";
      delml = true;
    }
    plot_ngpid(fout, id);
    plot_msd(fout, var);
  }
  if (opt.ngp("mse"))           // 平均値±SE
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -zx0," << var.ng() + 1 << ",1 -dline -w80\n";
      fout << "%エラーバーを描く時:-derrorbar_y -w80 -o500\n";
      delml = true;
    }
    plot_ngpid(fout, id);
    plot_mse(fout, var);
  }
  if (opt.ngp("sd"))            // 平均値±係数×SDの折れ線
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -zx0," << var.ng() + 1 << ",1 -dline -w80\n";
      fout << "%エラーバーを描く時:-derrorbar_y -w80 -o500\n";
      delml = true;
    }
    plot_ngpid(fout, id);
    plot_sdt(fout, opt, var);
  }
  if (opt.ngp("sr"))            // 平均値±係数×残差SD
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -zx0," << var.ng() + 1 << ",1 -dline -w80\n";
      fout << "%エラーバーを描く時:-derrorbar_y -w80 -o500\n";
    }
    plot_ngpid(fout, id);
    plot_srt(fout, opt, var);
  }
}

void anova1class::plot_ngpid(ofstream& fout, idd& id) const
// Ngraph用データのタイトルとデータ識別情報を出力する
{
  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 anova1class::plot_bar(ofstream& fout, const varclass& var) const
// Ngraph用にバーチャートデータを出力する
{
  fout << "%バーチャート\n";
  for (size_t i = 0; i < var.ng(); ++i)
  {
    double x = (double)(i + 1);
    fout << x - 0.2 << " 0 " << x + 0.2 << " " << var.mean(i) << "\n";
  }
}

struct pro                      // 最小・最大値と髭の下限・上限と四分位値
{
  double min, wlow, q1, q2, q3, wup, max;
};

void anova1class::plot_box(ofstream& fout, varclass& var) const
// Ngraph用に箱ヒゲプロットデータを出力する
{
  size_t ng = var.ng();
  vector<struct pro> provec(ng);
  fout << "%箱ヒゲプロット：箱\n";
  for (size_t i = 0; i < ng; ++i)
  {
    int n = var.n(i);
    vecdbl fdvec(n);
    (var.fdlst(i)).sort();      // 群のデータリストをソートする
    int j = 0;
    for (listdbl::const_iterator fditr = (var.fdlst(i)).begin(); fditr != (var.fdlst(i)).end(); ++fditr, ++j) fdvec[j] = *fditr;
    provec[i].min = provec[i].wlow = fdvec[0];
    provec[i].max = provec[i].wup = fdvec[n - 1];
    provec[i].q1 = fdvec[0];  provec[i].q2 = fdvec[0];
    provec[i].q3 = fdvec[n - 1];
    if (n > 1)
    {
      double q1fidx = (double)(n - 1) * 0.25;
      double q2fidx = (double)(n - 1) * 0.5;
      double q3fidx = (double)(n - 1) * 0.75;
      int q1idx = (int)q1fidx, q2idx = (int)q2fidx, q3idx = (int)q3fidx;
      provec[i].q1 = fdvec[q1idx] + (fdvec[q1idx + 1] - fdvec[q1idx]) * (q1fidx - (double)q1idx);
      provec[i].q2 = fdvec[q2idx] + (fdvec[q2idx + 1] - fdvec[q2idx]) * (q2fidx - (double)q2idx);
      provec[i].q3 = fdvec[q3idx] + (fdvec[q3idx + 1] - fdvec[q3idx]) * (q3fidx - (double)q3idx);
    }
    double h = provec[i].q3 - provec[i].q1;
    double hlow = provec[i].q1 - h * 1.5, hup = provec[i].q3 + h * 1.5;
    bool schwlow = true, schwup = true;
    for (j = 0; (schwlow || schwup) && j < n; ++j)
    {
      if (schwlow && hlow <= fdvec[j])
      { schwlow = false;  provec[i].wlow = fdvec[j]; }
      if (fdvec[j] <= hup) provec[i].wup = fdvec[j];
      else schwup = false;
    }
    double x = (double)i + 0.9;
    fout << x << " " << provec[i].q1 << " ";
    fout << x + 0.2 << " " << provec[i].q2 << "\n";
    fout << x << " " << provec[i].q2 << " ";
    fout << x + 0.2 << " " << provec[i].q3 << "\n";
  }
  fout << NGPdelm2 << "\n";
  fout << "%箱ヒゲプロット：ヒゲ(-zx0," << ng + 1 << ",1 -derrorbar_y -x1 -y2 -w40 -o500)\n";
  for (size_t i = 0; i < ng; ++i)
  {
    fout << i + 1 << " " << provec[i].q1;
    fout << " " << provec[i].wlow - provec[i].q1 << " 0\n";
    fout << i + 1 << " " << provec[i].q3;
    fout << " 0 " << provec[i].wup - provec[i].q3 << "\n";
  }
  fout << NGPdelm2 << "\n";
  fout << "%箱ヒゲプロット：平均値(-s1 -zx0," << ng + 1 << ",1 -dmark,72 -w80 -o400)\n";
  for (size_t i = 0; i < ng; ++i) fout << i + 1 << " " << var.mean(i) << "\n";
  fout << NGPdelm2 << "\n";
  fout << "%箱ヒゲプロット：外れ値(-s1 -zx0," << ng + 1 << ",1 -dmark,2 -w80 -o400)\n";
  for (size_t i = 0; i < ng; ++i)
  {
    for (double fdata : var.fdlst(i))
    {
      if (fdata < provec[i].wlow) fout << i + 1 << " " << fdata << "\n";
      else if (provec[i].wup < fdata) fout << i + 1 << " " << fdata << "\n";
    }
  }
}

void anova1class::plot_ci(ofstream& fout, option& opt, const varclass& var) const
// Ngraph用に平均±信頼区間データを出力する
{
  size_t ng = var.ng();
  fout << "%" << opt.cc() * 100.0 << "%信頼区間(-zx0," << ng + 1 << ",1 -derrorbar_y -w0 -o500)\n";
  for (size_t i = 0; i < ng; ++i)
  {
    double t = (double)(ptot((ldbl)(1.0 - opt.cc()), var.n(i) - 1));
    double ci = t * var.se(i);
    fout << i + 1 << " " << var.mean(i) << " " << ci << " " << -ci << "\n";
  }
}

void anova1class::plot_data(ofstream& fout, varclass& var) const
// Ngraph用に個々のデータを出力する
{
  fout << "%生データ\n";
  for (size_t i = 0; i < var.ng(); ++i)
  {
    for (double fdata : var.fdlst(i)) fout << i + 1 << " " << fdata << "\n";
  }
}

void anova1class::plot_msd(ofstream& fout, const varclass& var) const
// Ngraph用に平均±SDを出力する
{
  size_t ng = var.ng();
  fout << "%平均値±SD(-zx0," << ng + 1 << ",1 -derrorbar_y -w40 -o500)\n";
  for (size_t i = 0; i < ng; ++i)
  {
    double sd = var.sd(i);
    fout << i + 1 << " " << var.mean(i) << " " << sd << " " << -sd << "\n";
  }
}

void anova1class::plot_mse(ofstream& fout, const varclass& var) const
// Ngraph用に平均±SEを出力する
{
  size_t ng = var.ng();
  fout << "%平均値±SE(-zx0," << ng + 1 << ",1 -derrorbar_y -w40 -o500)\n";
  for (size_t i = 0; i < ng; ++i)
  {
    double se = var.se(i);
    fout << i + 1 << " " << var.mean(i) << " " << se << " " << -se << "\n";
  }
}

void anova1class::plot_sdt(ofstream& fout, option& opt, const varclass& var) const
// Ngraph用に平均±係数×SDを出力する
{
  int dfr = var.dfr();          // 残差自由度
  double t = (opt.tv() != BLANKf) ? opt.tv() : (double)(ptot((ldbl)(1.0 - opt.cc()), dfr));
  fout << "%平均値±" << t << "×SD\n";
  for (size_t i = 0; i < var.ng(); ++i) fout << i + 1 << " " << var.mean(i) << "\n";
  fout << NGPdelm2 << "\n";
  fout << "%下限\n";
  for (size_t i = 0; i < var.ng(); ++i)
  { fout << i + 1 << " " << var.mean(i) - t * var.sd(i) << "\n"; }
  fout << NGPdelm2 << "\n";
  fout << "%上限\n";
  for (size_t i = 0; i < var.ng(); ++i)
  { fout << i + 1 << " " << var.mean(i) + t * var.sd(i) << "\n"; }
}

void anova1class::plot_srt(ofstream& fout, option& opt, const varclass& var) const
// Ngraph用に平均±係数×残差SDを出力する
{
  int dfr = var.dfr();          // 残差自由度
  double t = (opt.tv() != BLANKf) ? opt.tv() : (double)(ptot((ldbl)(1.0 - opt.cc()), dfr));
  double tsd = t * sqrt(var.vrr());
  fout << "%平均値±" << t << "×残差SD=平均値±" << tsd << "\n";
  for (size_t i = 0; i < var.ng(); ++i) fout << i + 1 << " " << var.mean(i) << "\n";
  fout << NGPdelm2 << "\n";
  fout << "%下限\n";
  for (size_t i = 0; i < var.ng(); ++i)
  { fout << i + 1 << " " << var.mean(i) - tsd << "\n"; }
  fout << NGPdelm2 << "\n";
  fout << "%上限\n";
  for (size_t i = 0; i < var.ng(); ++i)
  { fout << i + 1 << " " << var.mean(i) + tsd << "\n"; }
}

int anova1class::get(const char* c_fname, option& opt, idd& id, ldata& datalst)
// データを読み込む
{
  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 anova1class::put(option& opt, idd& id, ldata& datalst) const
// 解析結果を出力する
{
  varclass var;                 // 計算用変数クラス
  var.cal(datalst);
  if (opt.res("res"))
  {
    put_resid(cout, id);
    if (var.ng() < 1)
    { cout << "※有効なデータなし！\n"; }
    else
    {
      if (opt.res("stat"))      // 基礎統計量出力
      {
        cout << "\n";
        put_stat(cout, var);
      }
      if (opt.res("mci"))       // 平均値の信頼区間出力
      {
        cout << "\n";
        put_mci(cout, opt, var);
      }
      if (var.ng() < 2) cout << "※分散分析不可能！\n";
      else if (var.ng() == 2) put_res2g(cout, opt, var);
      else put_resmg(cout, opt, var); // 3群以上の時
    }
    cout << "\n" << endl;
  }
  if (var.ng() > 0 && opt.ngpname() != "") plot_ngp(opt, id, datalst, var); // Ngraph用データ出力
}

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