// anovai.C Version 7.45  Copyleft (c) 杉本典夫
// Time-stamp: <2026-04-14 08:16:46 tac>
//
// 対応の無い多群間で一元配置分散分析を行い、結果をDANSデータとして出力する。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -W anovai anovai.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 =            // プログラム使用法表示用
{
  "anovai  Version 7.45         Copyleft(Vo\\oV) 2026 杉本典夫",
  "書式:anovai  [<オプション> …]  [<データファイル> …]",
  "機能:最初の項目を群項目(要因A)とし、以後の項目を集計項目とした一元配置分散",
  "     分析を行い、結果を連続出力します。",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定します。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -f=<ファイル名>:指定したファイルからオプションを読み込みます。",
  "  -dnsname[=<出力ファイル名>]:DANSデータを出力するファイル名を指定します。",
  "    オプション省略時:DANSデータ出力せず  -dnsname:-dnsname=anovai_dns.d",
  "  -dns={rchart|stat}:DANSデータの出力指定です。",
  "    rchart:rchart(レーダーチャート作成コマンド)用データ",
  "    stat:群別統計量(オプション省略時)",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示します。(-h:-h=man)",
  "  -res=[-][{def|all|at|ft|m|r|stat|tt}][,…]:解析結果の出力項目を指定します。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    at:分散分析  ft:等分散性の両側F検定  m:t検定の多重比較  r:相関比",
  "    stat:基礎統計量  tt:2標本t検定"
};

const string Title = "=== 多群の平均値の比較 ==="; // 表題
const string DNSNAMEdef = "anovai_dns.d"; // 省略時DANSファイル名

enum dnstype                    // DANSデータの種類
{
  dnsnon,                       // 出力なし
  dnsrchart,                    // rchart用データ
  dnsstat,                      // 統計量
  dnsetc                        // その他
};

using vecstatv = vector<statv<double>>; // 統計量ベクトル配列の別名定義

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  string v_dnsname;             // DANSファイル名
  dnstype v_dns;                // DANSデータの種類
  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; }
  dnstype dns() const           // DANSデータの種類を返す
  { return v_dns; }
  const string& dnsname() const // DANSデータファイル名を返す
  { return v_dnsname; }
  bool res(const string& code); // 分析結果出力フラッグを返す
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check(size_t ni);         // オプションのチェックとセット
  int set_res();                // -resオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_cc(CCdef), v_dnsname(), v_dns(dnsstat), 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); // プログラム名削除
}

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] == "dns")     // DANSデータの種類指定
  {
    if (opt[1].empty() || opt[1] == "stat") v_dns = dnsstat;
    else if (opt[1] == "rchart") v_dns = dnsrchart;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "dnsname") // DANSファイル名指定
  {
    if (opt[1].empty()) v_dnsname = DNSNAMEdef;
    else v_dnsname = opt[1];
    if (v_dnsname.size() < 3) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "f")       // オプション指定ファイル
  { err = fgetopt(*this, opt[1], true); }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  else if (opt[0] == "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);
}

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

int option::check(size_t ni)    // オプションのセットとチェック
{
  if (ni < 2) return disperr(pname(), "項目が2つ以上必要です！");
  if (set_res() < 0) return -2; // -resオプションのセット
  return 0;                     // エラーなし
}

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["ft"] = resmap["stat"] = resmap["tt"] = true;
  resmap["at"] = resmap["m"] = resmap["r"] = 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["at"] = resmap["ft"] = resmap["m"] = resmap["r"] = resmap["stat"] = resmap["tt"] = true;
      }
      else if (respos == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適当！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

class varclass                  // 変数クラス
{
private:
  set<string> gset;             // 群コードセット
  vector<string> gcodevec;      // 群コードベクトル
  mat<double> memat;            // 項目別・群コード別平均値行列
  double v_memin;               // 最小平均値
  double v_memax;               // 最大平均値
  vecstatv gstatvec;            // 群別・全体基礎統計量ベクトル
  double v_sa;                  // 要因A平方和
  int v_dfa;                    // 要因A自由度
  double v_sr;                  // 残差平方和
  int v_dfr;                    // 残差自由度
  double v_st;                  // 全体平方和
  int v_dft;                    // 全体自由度
  double v_v1;                  // 群1の分散
  int v_df1;                    // 群1の自由度
  double v_v2;                  // 群2の分散
  int v_df2;                    // 群2の自由度
  double v_dm;                  // 2群の平均値の差
  double v_vdmt;                // 平均値の差の不偏分散(t検定)
  double v_vdmw;                // 平均値の差の不偏分散(Welch検定)
  int v_dfw;                    // Welch検定の自由度
public:
  varclass() : gset(), gcodevec(), memat(), v_memin(FDmax), v_memax(-FDmax), gstatvec(), v_sa(0.0), v_dfa(0), v_sr(0.0), v_dfr(0), v_st(0.0), v_dft(0), v_v1(0.0), v_df1(0), v_v2(0.0), v_df2(0), v_dm(0.0), v_vdmt(0.0), v_vdmw(0.0), v_dfw(0)
  // デフォールトコンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  size_t g() const              // 群数を返す
  { return gset.size(); }
  const string& gcode(size_t idx) const // 群コードを返す
  { return (idx < gcodevec.size()) ? gcodevec[idx] : BLANKs; }
  double me(size_t i, size_t j) const // 項目別・群別平均値を返す
  { return (i < memat.r() && j < memat.c()) ? memat(i, j) : 0.0; }
  double memin() const         // 最小平均値を返す
  { return v_memin; }
  double memax() const         // 最大平均値を返す
  { return v_memax; }
  int n(size_t idx) const       // 群の例数を返す
  { return (idx <= g()) ? gstatvec[idx].n() : 0; }
  double mean(size_t idx) const // 群の平均値を返す
  { return (idx <= g()) ? gstatvec[idx].mean() : 0.0; }
  double ss(size_t idx) const   // 群の平方和を返す
  { return (idx <= g()) ? gstatvec[idx].ss() : 0.0; }
  double vr(size_t idx) const   // 群の不偏分散を返す
  { return (idx <= g()) ? gstatvec[idx].vr() : 0.0; }
  double sa() const             // 要因A平方和を返す
  { return v_sa; }
  double va() const             // 要因A分散を返す
  { return (dfa() > 0) ? sa() / (double)(dfa()) : 0.0; }
  int dfa() const               // 要因A自由度を返す
  { return v_dfa; }
  double sr() const             // 残差平方和を返す
  { return v_sr; }
  double vr() const             // 残差分散を返す
  { return (dfr() > 0) ? sr() / (double)(dfr()) : 0.0; }
  int dfr() const               // 残差自由度を返す
  { return v_dfr; }
  double st() const             // 全体平方和を返す
  { return v_st; }
  int dft() const               // 全体自由度を返す
  { return v_dft; }
  double v1() const             // 群1の分散を返す
  { return v_v1; }
  int df1() const               // 群1の自由度を返す
  { return v_df1; }
  double v2() const             // 群2の分散を返す
  { return v_v2; }
  int df2() const               // 群2の自由度を返す
  { return v_df2; }
  double dm() const             // 2群の平均値の差を返す
  { return v_dm; }
  double vdmt() const           // 平均値の差の不偏分散を返す(t検定)
  { return v_vdmt; }
  double vdmw() const           // 平均値の差の不偏分散を返す(Welch検定)
  { return v_vdmw; }
  int dfw() const               // Welch検定の自由度を返す
  { return v_dfw; }
  void setg(size_t ni, ldata& datalst); // 群コードを分類する
  void cala1(ldata& datalst, int jidx); // 分散分析用計算をする
  double caltt();               // 2標本t検定とWelchの検定用計算を行う
};

void varclass::setg(size_t ni, ldata& datalst)
// 群コードを分類する
{
  for (auto& cdatai : datalst)
  {
    if (cdatai.sdata(0) != BLANKs) gset.insert(cdatai.sdata(0));
  }
  gcodevec.resize(g());  gstatvec.resize(g() + 1);
  size_t i = 0;
  for (auto& gcode : gset) gcodevec[i++] = gcode;
  memat.resize(g(), ni - 1);
  memat = 0.0;
}

void varclass::cala1(ldata& datalst, int jidx)
// 分散分析用計算をする
{
  for (size_t i = 0; i <= g(); ++i) gstatvec[i].clear();
  for (auto& cdatai : datalst)
  {
    if (cdatai.sdata(0) != BLANKs && cdatai.fdata(jidx) != BLANKf)
    {
      int i = search_data(gset, cdatai.sdata(0));
      gstatvec[i].sigma(cdatai.fdata(jidx));
    }
  }
  v_sa = 0.0;
  int ng = 0;
  for (size_t i = 0; i < g(); ++i)
  {
    int n = gstatvec[i].n();
    if (n > 0)
    {
      double m = gstatvec[i].mean();
      v_sa += (double)n * m * m;
      gstatvec[g()] += gstatvec[i];
      ++ng;
      memat.assign(i, jidx - 1, gstatvec[i].mean());
      if (m < v_memin) v_memin = m; // 最小平均値のサーチ
      if (m > v_memax) v_memax = m; // 最大平均値のサーチ
    }
    else
    { memat.assign(i, jidx - 1, BLANKf); }
  }
  v_dfa = ng - 1;  v_dft = gstatvec[g()].n() - 1;
  v_dfr = v_dft - v_dfa;
  v_sa -= (double)gstatvec[g()].n() * gstatvec[g()].mean() * gstatvec[g()].mean();
  v_st = gstatvec[g()].ss();
  if (v_dfa < 1 || v_sa < FDeps) v_sa = 0.0;
  v_sr = v_st - v_sa;
  if (v_dft < 1 || v_sr < FDeps) v_sr = 0.0;
}

double varclass::caltt()        // 2標本t検定とWelchの検定用計算を行う
{
  int idx1 = -1, idx2 = -1;
  bool sch = true;
  for (size_t i = 0; sch && i < g(); ++i)
  {
    if (gstatvec[i].n() > 0)
    {
      if (idx1 < 0) idx1 = i;      // 群1のインデックス
      else if (idx2 < 0) idx2 = i; // 群2のインデックス
      else sch = false;
    }
  }
  int n1 = gstatvec[idx1].n(), n2 = gstatvec[idx2].n();
  double me1 = gstatvec[idx1].mean(), me2 = gstatvec[idx2].mean();
  double ss1 = gstatvec[idx1].ss(), ss2 = gstatvec[idx2].ss();
  v_v1 = gstatvec[idx1].vr();  v_v2 = gstatvec[idx2].vr();
  double f = (sgn(v_v1) > 0) ? v_v2 / v_v1 : 0.0;
  v_df1 = n1 - 1;  v_df2 = n2 - 1;
  if (v_v2 < v_v1)
  {
    f = (sgn(v_v2) > 0) ? v_v1 / v_v2 : 0.0;
    v_df1 = n2 - 1;  v_df2 = n1 - 1;
  }
  v_dm = me1 - me2;
  v_dft = n1 + n2 - 2;
  v_vdmt = (ss1 + ss2) / (double)v_dft * (double)(n1 + n2) / (double)n1 / (double)n2;
  double vn1 = v_v1 / (double)n1, vn2 = v_v2 / (double)n2;
  v_vdmw = vn1 + vn2;
  v_dfw = (int)(v_vdmw * v_vdmw / (vn1 * vn1 / (double)v_df1 + vn2 * vn2 / (double)v_df2));
  return f;                     // 分散比Fを返す
}

class anovaiclass               // anovaiメインクラス
{
private:
  option opt;                   // オプション
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout, int iidx);
  // 結果のタイトルとデータ識別情報を出力する
  void put_anova1(ostream& fout, varclass& var);
  // 分散分析の結果を分散分析表形式で出力する
  void put_stat(ostream& fout, const varclass& var) const;
  // 群別基礎統計量を出力する
  void put_rratio(ostream& fout, const varclass& var);
  // 相関比を出力する
  void put_ttest(ostream& fout, ldbl prf, const varclass& var);
  // 2標本t検定結果とWelchの検定結果を出力する
  void put_ftest(ostream& fout, double f, int df1, int df2, ldbl prf);
  // 分散比Fの両側検定と推定結果を出力する
  void put_dnsid(ofstream& fout, option& opt, const varclass& var);
  // DANSデータのタイトルとデータ識別情報を出力する
  void put_dnsstat(ofstream& fout, const varclass& var, int iidx);
  // 分散分析の統計量をケースデータ形式で出力する
  void put_dnsrchart(ofstream& fout, option& opt, const varclass& var);
  // rchart用データをケースデータ形式で出力する
public:
  anovaiclass(const char* c_pname = nullptr) : opt(c_pname), fname(), id(), datalst()
  // コンストラクタ
  {}
  ~anovaiclass()                // デストラクタ
  {}
  const string& pname() const   // プログラム名を返す
  { return opt.pname(); }
  void init()      // 初期設定を行う(コンストラクタ中ではできない処理)
  { opt.setdef(); }
  void setopt(const string& optstr, char devc = OPTdev) // オプションを設定する
  { opt.set(optstr, devc); }
  int get(const char* c_fname); // データを読み込み群コードを集計する
  void put(int nc);             // 集計結果を出力する
};

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

void anovaiclass::put_anova1(ostream& fout, varclass& var)
// 分散分析の結果を分散分析表形式で出力する
{
  if (opt.res("stat"))
  {
    fout << "\n";
    put_stat(fout, var);
  }
  if (var.dfa() < 1)
  { fout << "※分散分析不可能！\n"; }
  else if (var.dfa() == 1)      // 2群の場合
  {
    if (opt.res("r"))
    {
      fout << "\n";
      put_rratio(fout, var);
    }
    if (opt.res("at"))
    {
      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"))
    {
      double f = var.caltt();
      ldbl prf = ftop(f, var.df1(), var.df2());
      prf = (prf < 0.5L) ? 2.0L * prf : 1.0L; // 両側p値にする
      if (opt.res("ft"))
      {
        fout << "\n";
        put_ftest(fout, f, var.df1(), var.df2(), prf);
      }
      fout << "\n";
      put_ttest(fout, prf, var);
    }
  }
  else                          // 3群以上の場合
  {
    if (opt.res("def") || opt.res("r"))
    {
      fout << "\n";
      put_rratio(fout, var);
    }
    if (opt.res("def") || opt.res("at"))
    {
      fout << "\n・一元配置分散分析(one-way layout analysis of variance)\n";
      putatab(fout, "群(要因A)", var.sa(), var.dfa(), var.sr(), var.dfr(), var.st(), var.dft());
    }
  }
  fout << "\n" << endl;
}

void anovaiclass::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.g(); ++i)
  {
    if (var.n(i) > 0)
    {
      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.g()), var.mean(var.g()), var.vr(var.g()));
}

void anovaiclass::put_rratio(ostream& fout, const varclass& var)
// 相関比を出力する
{
  double r2 = var.sa() / var.st();
  fout << "・相関比:η=" << sqrt(r2) << "  寄与率:η2乗=" << r2 << "\n";
  double f = (sgn(var.vr()) > 0) ? var.va() / var.vr() : 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 anovaiclass::put_ttest(ostream& fout, ldbl prf, const varclass& var)
// 2標本t検定結果とWelchの検定結果を出力する
{
  fout << "・対応のないt検定(2標本t検定、「→」の付いた方を採用して下さい)\n";
  double sd = sqrt((var.ss(0) + var.ss(1)) / (double)var.dft());
  double d = 0.0;
  if (sd < FDeps) sd = 0.0;
  else d = var.dm() / sd;
  fout << "  平均値の差=" << var.dm() << "  平均値の差のSD=" << sd << "  効果量=" << d << "\n";
  string sel1 = "  →", sel2 = "    ";
  if (prf <= 0.05L) { sel1 = "    ";  sel2 = "  →"; }
  double t = var.dm() / sqrt(var.vdmt());
  ldbl pr = ttop(t, var.dft());
  fout << sel1 << "等分散t検定(=分散分析)    :t値=" << t;
  fout << "  自由度=" << var.dft() << "  p=" << pr << pind(pr) << "\n";
  if (id.ni() > 2 && opt.res("m"))
  {
    pr = pr * (ldbl)(id.ni() - 1);
    if (pr > 1.0L) pr = 1.0L;
    fout << sel1 << "Bonferroni型多重検定:有意確率p=" << pr << pind(pr) << "\n";
  }
  double ci = calmcci(mcfis, opt.cc(), 2, 1, var.dft(), 1, var.vdmt());
  fout << sel1 << "平均値の差の" << opt.cc() * 100.0 << "%信頼区間=";
  fout << var.dm() << "±" << ci << "(" << var.dm() - ci << "-" << var.dm() + ci << ")\n";
  if (id.ni() > 2 && opt.res("m"))
  {
    ci = calmcci(mcbon, opt.cc(), 2, 1, var.dft(), id.ni() - 1, var.vdmt());
    fout << sel1 << "平均値の差の" << opt.cc() * 100.0 << "%同時信頼区間=";
    fout << var.dm() << "±" << ci << "(" << var.dm() - ci << "-" << var.dm() + ci << ")\n";
  }
  t = var.dm() / sqrt(var.vdmw());
  pr = ttop(t, var.dfw());
  fout << sel2 << "不等分散t検定(Welchの検定):t値=" << t;
  fout << "  自由度=" << var.dfw() << "  p=" << pr << pind(pr) << "\n";
  if (id.ni() > 2 && opt.res("m"))
  {
    pr = pr * (ldbl)(id.ni() - 1);
    if (pr > 1.0L) pr = 1.0L;
    fout << sel2 << "Bonferroni型多重検定:有意確率p=" << pr << pind(pr) << "\n";
  }
  ci = calmcci(mcfis, opt.cc(), 2, 1, var.dfw(), 1, var.vdmw());
  fout << sel2 << "平均値の差の" << opt.cc() * 100.0 << "%信頼区間=";
  fout << var.dm() << "±" << ci << "(" << var.dm() - ci << "-" << var.dm() + ci << ")\n";
  if (id.ni() > 2 && opt.res("m"))
  {
    ci = calmcci(mcbon, opt.cc(), 2, 1, var.dfw(), id.ni() - 1, var.vdmw());
    fout << sel2 << "平均値の差の" << opt.cc() * 100.0 << "%同時信頼区間=";
    fout << var.dm() << "±" << ci << "(" << var.dm() - ci << "-" << var.dm() + ci << ")\n";
  }
}

void anovaiclass::put_ftest(ostream& fout, double f, int df1, int df2, ldbl prf)
// 分散比Fの両側検定と推定結果を出力する
{
  fout << "・等分散性の両側F検定\n";
  fout << "  F=" << f << " 第1自由度=" << df1 << " 第2自由度=" << df2 << "  両側p=" << prf << pind(prf) << "\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 anovaiclass::put_dnsid(ofstream& fout, option& opt, const varclass& var)
// DANSデータのタイトルとデータ識別情報を出力する
{
  int tni = (opt.dns() == dnsrchart) ? id.ni() + 3 : 2 * (var.g() + 1) + 7;
  idd tid(id.title(), tni);
  if (opt.dns() == dnsrchart)
  {
    tid.assign_item(0, id.item(0));
    for (size_t i = 1; i < id.ni(); ++i) tid.assign_item(i, id.item(i) + "の平均値");
    tid.assign_item(id.ni(), "座標軸の最小値");
    tid.assign_item(id.ni() + 1, "座標軸の最大値");
    tid.assign_item(id.ni() + 2, "座標軸の目盛り間隔");
  }
  else
  {
    int j = 1;
    tid.assign_item(0, "項目名");
    for (size_t i = 0; i < var.g(); ++i)
    {
      tid.assign_item(j++, id.item(0) + ":" + var.gcode(i) + "の例数");
      tid.assign_item(j++, id.item(0) + ":" + var.gcode(i) + "の平均値");
    }
    tid.assign_item(j++, "全体の例数");
    tid.assign_item(j++, "全体の平均値");
    tid.assign_item(j++, "要因Aの分散");
    tid.assign_item(j++, "残差分散");
    tid.assign_item(j++, "総分散");
    tid.assign_item(j++, "F値");
    tid.assign_item(j++, "p値");
    tid.assign_item(j++, "相関比");
  }
  tid.put(fout, formtag);
}

void anovaiclass::put_dnsstat(ofstream& fout, const varclass& var, int iidx)
// 分散分析の統計量をケースデータ形式で出力する
{
  int tni = 2 * (var.g() + 1) + 7;
  datac cdata(tni);
  cdata.assign(0, id.item(iidx));
  double fa = (sgn(var.vr()) > 0) ? var.va() / var.vr() : 0.0;
  double r2 = (sgn(var.st()) > 0) ? var.sa() / var.st() : 0.0;
  double pr = ftop(fa, var.dfa(), var.dfr());
  int j = 1;
  for (size_t i = 0; i < var.g(); ++i)
  {
    cdata.assign(j++, (double)(var.n(i)));
    if ((double)(var.n(i)) > 0) cdata.assign(j++, var.mean(i));
    else cdata.assign(j++, BLANKs);
  }
  cdata.assign(j++, (double)(var.n(var.g())));
  cdata.assign(j++, var.mean(var.g()));
  cdata.assign(j++, var.va());
  cdata.assign(j++, var.vr());
  cdata.assign(j++, var.vr(var.g()));
  cdata.assign(j++, fa);
  cdata.assign(j++, pr);
  cdata.assign(j++, sqrt(r2));
  cdata.put(fout);
}

void anovaiclass::put_dnsrchart(ofstream& fout, option& opt, const varclass& var)
// rhcart用データをケースデータ形式で出力する
{
  double arange = var.memax() - var.memin();
  double da = arange * 0.1, ainc = arange * 1.2 / 4.0; // 目盛り間隔は4つ
  double amin = var.memin() - da, amax = var.memax() + da;
  int tni = id.ni() + 3;
  datac cdata(tni);
  for (size_t i = 0; i < var.g(); ++i)
  {
    cdata.assign(0, var.gcode(i));
    for (size_t j = 1; j < id.ni(); ++j) cdata.assign(j, var.me(i, j - 1));
    cdata.assign(id.ni(), amin); // 座標軸の最小値
    cdata.assign(id.ni() + 1, amax); // 座標軸の最大値
    cdata.assign(id.ni() + 2, ainc); // 座標軸の目盛り間隔
    cdata.put(fout);
  }
}

int anovaiclass::get(const char* c_fname) // データを読み込み群コードを集計する
{
  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 anovaiclass::put(int nc)   // 集計結果を出力する
{
  varclass var;                 // 計算用変数クラス
  var.setg(id.ni(), datalst);
  if (var.g() < 1) return;      // 群コードがない
  ofstream foutdns;
  if (opt.dnsname() != "")      // DANSデータ出力
  {
    if (ofopen(foutdns, (opt.dnsname()).c_str()) == false) return; // ファイルオープン失敗
    put_dnsid(foutdns, opt, var);
  }
  for (size_t j = 1; j < id.ni(); ++j)
  {
    var.cala1(datalst, j);
    if (opt.res("res"))         // 集計結果を分散分析表形式で出力する
    {
      put_resid(cout, j);
      put_anova1(cout, var);
    }
    if (opt.dnsname() != "" && opt.dns() == dnsstat) put_dnsstat(foutdns, var, j); // 集計結果をケースデータ形式で出力する
  }
  if (opt.dnsname() != "" && opt.dns() == dnsrchart) put_dnsrchart(foutdns, opt, var);
}

int main(int argc, char* argv[]) // メインルーチン
{
  anovaiclass anovai(argv[0]);  // メインクラス変数
  try
  {
    anovai.init();              // anovaiの初期設定
    bool fok = true;            // ファイル処理フラッグ
    int nc = 0;                 // 読み込んだケース数
    for (int i = 1; i <= argc; ++i)
    {
      if (i >= argc)
      {
        if (fok) nc += anovai.get(nullptr); // 標準入力から読み込む
      }
      else if (*(argv[i]) == CHRopt)
      { anovai.setopt(argv[i] + 1);  fok = true; }
      else
      { nc += anovai.get(argv[i]);  fok = false; }
    }
    anovai.put(nc);
  }
  catch (help_except& helpex)   // ヘルプファイル表示
  { return 0; }
  catch (option_error& opterr)  // オプション指定エラー
  { return 1; }
  catch (runtime_error& rterr)  // 実行時エラー
  {
    disperr(anovai.pname(), rterr.what());
    return 2;
  }
  catch (...)                   // 予期せぬ例外
  {
    disperr(anovai.pname(), "内容不明なエラー発生！");
    return 255;
  }
  return 0;
}
