// lma.C Version 7.45  Copyleft (c) 杉本典夫
// Time-stamp: <2026-04-27 17:15:17 tac>
//
// 最初の項目を目的変数Y、それ以後の項目を説明変数Xとした
// 線形モデルによる分散分析・共分散分析
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -W lma lma.C
//               -I/usr/lcoal/include -L/usr/local/lib -ldans -lfun

#include <dans.H>               // dans用ヘッダファイル(fun.H,mat.H含む)

using namespace dans;           // dans名前空間使用(名前空間matvec含む)

const vecstr USAGE =            // プログラム使用法表示用
{
  "lma  Version 7.45          Copyleft(Vo\\oV) 2026 杉本典夫",
  "書式:lma [<オプション> …]  [<データファイル> …]",
  "機能:最初の項目を目的変数Yとし、以後の項目を説明変数Xとした",
  "     線形モデルによる分散分析または共分散分析を行います。",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定します。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -f=<ファイル>:オプションを指定したファイルから読み込みます。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示します。(-h:-h=man)",
  "  -ino=<説明変数の開始番号>[-<終了番号>][,…]",
  "    交互作用の要因にする説明変数の番号を指定します。",
  "  -ix=<説明変数の開始番号>[-<終了番号>]x<説明変数の開始番号>[-<終了番号>][,…]",
  "    交互作用項目を作成する2つの説明変数の番号を指定します。",
  "  -ixc[={-|+}]:交互作用項目の説明変数を中心化するかどうか指定します。",
  "    -:中心化しない(オプション省略時)  +:中心化する  -ixc:ixc=+",
  "  -mf=<説明変数の個数>:主効果扱いする説明変数の個数を指定します。",
  "  -mno=<説明変数の開始番号>[-<終了番号>][,…]",
  "    主効果の要因にする説明変数の番号を指定します。",
  "  -res=[-][{def|all|atab|reg|stat}][,…]:解析結果の出力項目を指定します。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    atab:分散分析表  reg:重回帰分析結果  stat:基礎統計量"
};

const string Title = "=== 線形モデルによる分散分析 ==="; // 表題
const string NGPNAMEdef = "lma_ngp.dat"; // 省略時Ngraphファイル名

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  lino v_inolst;                // 交互作用の説明変数の番号リスト
  lstr ixstrlst;                // 交互作用項目の説明変数リスト
  lino v_ixlst;                 // 交互作用項目の説明変数ペアリスト
  bool v_ixc;                   // 交互作用項目の中心化フラッグ
  int v_mf;                     // 主効果の説明変数の個数
  lino v_mnolst;                // 主効果の説明変数の番号リスト
  lstr reslst;                  // 解析結果出力リスト
  mapsb resmap;                 // 解析結果出力の有無マップ
public:
  option(const char* c_pname = nullptr); // コンストラクタ
  ~option()                              // デストラクタ
  {}
  const string& path() const    // プログラムのパスを返す
  { return v_path; }
  const string& pname() const   // プログラム名を返す
  { return v_pname; }
  double cc() const             // 信頼係数を返す
  { return v_cc; }
  const lino& inolst()         // 交互作用の説明変数の番号リストを返す
  { return v_inolst; }
  const lino& ixlst()    // 交互作用項目の説明変数番号ペアリストを返す
  { return v_ixlst; }
  bool ixc()                   // 交互作用項目の中心化フラッグを返す
  { return v_ixc; }
  int mf() const                // 主効果の説明変数の個数を返す
  { return v_mf; }
  const lino& mnolst()          // 主効果の説明変数の番号リストを返す
  { return v_mnolst; }
  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オプションのセット
  int set_ix(int nx);           // -ixオプションのセット
  int set_mx(int nx, int nxix); // -mf,-mnoオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_cc(CCdef), v_inolst(), ixstrlst(), v_ixlst(), v_ixc(false), v_mf(0), v_mnolst(), 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 resitr = resmap.find(code);
  if (resitr == 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] == "ino")     // 交互作用の説明変数の番号指定
  {
    if (getrangelst(v_inolst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "ix")      // 交互作用項目の説明変数ペア指定
  {
    if (strsplit(ixstrlst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "ixc")     // 交互作用項目の説明変数の中心化指定
  {
    if (opt[1].empty() || opt[1] == "+" ) v_ixc = true;
    else if (opt[1] == "-") v_ixc = false;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "mf")      // 主効果の説明変数の個数指定
  {
    char* ptr = nullptr;
    v_mf = strtol(opt[1].c_str(), &ptr, 0);
    if (v_mf < 1) err = badopt(pname(), optstr, USAGE, 1.0, (double)INT_MAX);
  }
  else if (opt[0] == "mno")     // 主効果の説明変数の番号指定
  {
    if (getrangelst(v_mnolst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "res")     // 分析結果の出力指定
  {
    if (strsplit(reslst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else                          // 間違ったオプション指定
  { err = nonopt(pname(), optstr, USAGE); }
  if (err > 0) throw help_except("ヘルプ表示");
  if (err < 0) throw option_error(string("オプション指定違反→") + optstr);
}

int option::check(size_t ni)    // オプションのチェックとセット
{
  if (ni < 2) return disperr(pname(), "項目数が足りません！");
  if (set_res() < 0) return -2; // -resオプションのセット
  int nx = (int)ni - 1;         // 交互作用項目を含まない説明変数の数
  int nxix = set_ix(nx);        // -ixオプションのセット
  if (nxix < 0) return -3;      // 交互作用項目を含んだ説明変数の数
  if (set_mx(nx, nxix) < 0) return -4; // -mf,-mnoオプションのセット
  return 0;                     // エラーなし
}

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["atab"] = resmap["reg"] = resmap["stat"] = true;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["atab"] = resmap["reg"] = resmap["stat"] = false;
    for (const string& resstr : reslst)
    {
      mapsb_citr resitr = resmap.find(resstr);
      if (resstr == "-")        // 全て出力せず
      { resmap["res"] = false; }
      else if (resstr[0] == '-' && resstr.size() > 1) // 特定の項目を出力せず
      {
        resitr = resmap.find(resstr.substr(1));
        if (resitr == resmap.end())
        { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
        else
        { resmap[resstr.substr(1)] = false; }
      }
      else if (resstr == "def" || resstr == "all")
      { resmap["def"] = resmap["atab"] = resmap["reg"] = resmap["stat"] = true; }
      else if (resitr == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適切！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

int option::set_ix(int nx)      // -ixオプションのセット
{
  if (ixstrlst.empty()) return nx;
  int ixno = nx + 1;            // 交互作用項目の最初の番号
  for (const string& ixstr : ixstrlst)
  {
    string opt[2];
    if (getopt(opt, ixstr, 'x') > 0)
    {
      int low0 = 0, up0 = 0, low1 = 0, up1 = 0;
      if (getrange(low0, up0, opt[0]) < 1)
      { return disperr(pname(), "1番目の説明変数の番号が間違っています→", ixstr); }
      else
      {
        if (getrange(low1, up1, opt[1]) < 1)
        { return disperr(pname(), "2番目の説明変数の番号が間違っています→", ixstr); }
        else
        {
          if (low0 < 1 || low0 > up0) return disperr(pname(), "1番目の説明変数の番号が間違っています→", ixstr);
          if (low1 < 1 || low1 > up1) return disperr(pname(), "2番目の説明変数の番号が間違っています→", ixstr);
          if (ixno <= up0 || ixno <= up1) return disperr(pname(), "説明変数の番号が間違っています→", ixstr);
          for (int j = low1; j <= up1; ++j) // 2番目の番号
          {
            for (int i = low0; i <= up0; ++i) // 1番目の番号
            {
              pairint ixpi {i, j};
              v_ixlst.push_back(ixpi);
            }
          }
          ixno += (up0 - low0 + 1) * (up1 - low1 + 1);
        }
      }
    }
  }
  return ixno - 1;              // 交互作用項目を含めた説明変数の数
}

int option::set_mx(int nx, int nxix) // -mf,-mnoオプションのセット
{
  if (v_mf < 1) v_mf = nx; // 交互作用項目以外の説明変数を全て主効果にする
  else if (nx < v_mf) return disperr(pname(), "主効果の説明変数の個数が間違っています！");
  int nt = 0;
  for (auto& mnopi : v_mnolst)
  {
    if (mnopi.first > mnopi.second || mnopi.second < 1 || v_mf < mnopi.first) return disperr(pname(), "主効果の説明変数の番号が間違っています！");
    nt += mnopi.second - mnopi.first + 1;
  }
  if (v_mf < nt) return disperr(pname(), "主効果の説明変数の合計個数が間違っています！");
  if (v_inolst.empty()) return v_mf; // 主効果の個数を返す
  for (auto& inopi : v_inolst)
  {
    if (inopi.first > inopi.second || inopi.second < 1 || nxix < inopi.second) return disperr(pname(), "交互作用の説明変数の番号が間違っています！");
    for (pairint mnopi : v_mnolst)
    {
      if (mnopi.first <= inopi.first && inopi.first <= mnopi.second) return disperr(pname(), "交互作用の説明変数の番号と主効果の説明変数の番号が重複しています！");
      if (mnopi.first <= inopi.second && inopi.second <= mnopi.second) return disperr(pname(), "交互作用の説明変数の番号と主効果の説明変数の番号が重複しています！");
      nt += mnopi.second - mnopi.first + 1;
    }
  }
  return nt;                    // 主効果＋交互作用効果の個数を返す
}

class varclass                  // 変数クラス
{
private:
  mat<double> v_smat;           // 積和行列
  mat<double> xmat;             // 単純積和行列
  mat<double> rmat;             // 相関行列
  vecdbl v_mevec;               // 平均値ベクトル
  vecdbl v_vrvec;               // 分散ベクトル
  vecdbl ssvec;                 // 積和ベクトル
  vecdbl fssvec;                // 線形モデル用積和ベクトル
  vecint fdfvec;                // 線形モデル用自由度ベクトル
  int v_ni;                     // 項目数
  int v_nxix;                   // 交互作用項目を含んだ説明変数の数
  int v_n;                      // 例数
  double v_sb;                  // 回帰平方和
  restype v_resvar;             // 処理結果の種類
  int set_mean(vecdbl& mevec, ldata& datalst);
  // 交互作用項計算用平均値ベクトルをセットする
  void cal_vif();               // VIFを計算する
public:
  varclass(size_t ni, option& opt); // コンストラクタ
  ~varclass()                   // デストラクタ
  {}
  const mat<double>& smat() const // 積和行列を返す
  { return v_smat; }
  const vecdbl& mevec() const   // 平均値ベクトルを返す
  { return v_mevec; }
  const vecdbl& vrvec() const   // 分散ベクトルを返す
  { return v_vrvec; }
  int ni() const                // 項目数を返す
  { return v_ni; }
  int nx() const                // 説明変数の数を返す
  { return ni() - 1; }
  int nxix() const              // 交互作用項目を含んだ説明変数の数を返す
  { return v_nxix; }
  int np() const                // パラメータ数を返す
  { return nxix() + 1; }
  int n() const                 // 例数を返す
  { return v_n; }
  int dft() const               // 全自由度を返す
  { return n() - 1; }
  int dfr() const               // 残差自由度を返す
  { return dft() - dfb(); }
  int dfb() const               // 回帰自由度を返す
  { return nxix(); }
  double syy() const            // 全平方和を返す
  { return ssvec[np()]; }
  double sr() const            // 残差平方和を返す
  { return v_smat(np(), np()); }
  double sb() const             // 回帰平方和を返す
  { return v_sb; }
  double vr() const             // 残差分散を返す
  { return (dfr() > 0) ? sr() / (double)(dfr()) : 0.0; }
  double b(int idx) const       // 偏回帰係数を返す
  { return (0 <= idx && idx <= np()) ? v_smat(idx, np()) : 0.0; }
  double vb(int idx) const      // 説明変数の分散を返す
  { return (0 <= idx && idx <= np()) ? vr() * v_smat(idx, idx) : 0.0; }
  double bst(int idx) const     // 標準偏回帰係数を返す
  { return (0 <= idx && idx <= np() && syy() > 0.0) ? b(idx) * sqrt(ssvec[idx] / syy()) : 0.0; }
  double rp(int idx) const;     // 偏相関係数を返す
  double vif(int idx) const     // VIFを返す
  { return (0 < idx && idx <= dfb()) ? rmat(idx, idx) : 0.0; }
  double r2() const             // 重寄与率を返す
  { return (ssvec[np()] > 0.0) ? sb() / ssvec[np()] : 0.0; }
  double r2a() const            // 自由度調整済み重寄与率を返す
  { return (ssvec[np()] > 0.0) ? fabs(1.0 - vr() * (double)(dft()) / ssvec[np()]) : 0.0; }
  double aic() const            // AICを返す
  { return (sr() > 0.0) ? (double)(n()) * log(sr() / (double)(n())) + (double)(2 * (dfb() + 2)) : 0.0; }
  double ss(int idx) const      // 線形モデル用積和を返す
  { return (0 <= idx && idx <= np() + 1) ? fssvec[idx] : 0.0; }
  int df(int idx) const         // 線形モデル用自由度を返す
  { return (0 <= idx && idx <= np() + 1) ? fdfvec[idx] : 0; }
  double me(int idx) const      // 平均値を返す
  { return (0 <= idx && idx <= np()) ? v_mevec[idx] : 0.0; }
  restype resvar() const        // 処理結果の種類を返す
  { return v_resvar; }
  void cal(option& opt, ldata& datalst); // データを集計して各種統計量を求める
  void callma(option& opt);              // 線形モデルによる分析
};

varclass::varclass(size_t ni, option& opt) : v_smat(), xmat(), rmat(), v_mevec(), v_vrvec(), ssvec(), fssvec(), fdfvec(), v_ni(ni), v_nxix(0), v_n(0), v_sb(0.0), v_resvar(resnon)
// コンストラクタ
{
  v_nxix = ni - 1 + (opt.ixlst()).size();
  size_t dim = v_nxix + 2;
  v_smat.resize(dim, dim);  xmat.resize(dim, dim);
  v_mevec.resize(dim, 0.0);  v_vrvec.resize(dim, 0.0);  ssvec.resize(dim, 0.0);
  fssvec.resize(dim + 1, 0.0);  fdfvec.resize(dim + 1, 0);
}

double varclass::rp(int idx) const // 偏相関係数を返す
{
  if (idx < 0 || idx > np()) return 0.0;
  double bi = b(idx);
  double v = sr() * v_smat(idx, idx) + bi * bi;
  return (v > 0.0) ? bi / sqrt(v) : 0.0;
}

void varclass::cal(option& opt, ldata& datalst)
// データを集計して各種統計量を求める
{
  vecdbl mevec(nx() + 1);       // 交互作用項計算用平均値ベクトル
  v_n = set_mean(mevec, datalst);
  if (v_n < np())               // 重回帰計算不可能
  {
    v_resvar = reserr;
    return;
  }
  int dim = np() + 1;
  vecdbl fdatavec(dim, 0.0);
  for (auto& cdatai : datalst)
  {
    if (cdatai.fdata(0) != BLANKf)
    {
      setixdata(fdatavec, opt.ixlst(), opt.ixc(), mevec, cdatai);
      for (int i = 0; i < dim; ++i)
      {
        double x = fdatavec[i];
        for (int j = i; j < dim; ++j) v_smat.add(i, j, x * fdatavec[j]);
      }
    }
  }
  for (int i = 0; i < np(); ++i)
  {
    xmat.assign(i, i, v_smat(i, i));
    for (int j = i + 1; j <= np(); ++j)
    {
      v_smat.assign(j, i, v_smat(i, j));
      if (j < np())
      { xmat.assign(i, j, v_smat(i, j));  xmat.assign(j, i, v_smat(i, j)); }
    }
  }
  v_smat.sweep(0);
  for (int i = 0; i < dim; ++i)
  {
    v_mevec[i] = v_smat(0, i);
    v_vrvec[i] = (dft() > 0) ? v_smat(i, i) / (double)(dft()) : 0.0;
    ssvec[i] = v_smat(i, i);
  }
  v_resvar = resok;
}

int varclass::set_mean(vecdbl& mevec, ldata& datalst)
// 交互作用項計算用平均値ベクトルをセットする
{
  int ni = nx() + 1;            // 項目数
  mevec.assign(ni, 0.0);
  vecdbl fdatavec(ni, 0.0);
  int n = 0;
  for (auto& cdatai : datalst)
  {
    bool blank = false;
    int j = 0;
    while (blank == false && j < ni)
    {
      fdatavec[j] = cdatai.fdata(j);
      if (fdatavec[j] == BLANKf) blank = true;
      else ++j;
    }
    if (blank)
    { cdatai.assign(0, BLANKs); } // 後の計算のために目的変数を欠損にする
    else
    {
      ++n;
      for (j = 0; j < ni; ++j) mevec[j] += fdatavec[j];
    }
  }
  for (int j = 0; j < ni; ++j) mevec[j] /= (double)n;
  return n;                     // 有効な例数を返す
}

void varclass::callma(option& opt)
// 線形モデルによる分析
{
  int iend = opt.mf() + 1;
  for (int i = 1; i < iend; ++i) v_smat.sweep(i);
  int idx = 0;
  if ((opt.mnolst()).empty() == false)
  {
    mat<double> tsmat(v_smat);
    double sr = tsmat(np(), np());
    for (auto& mnoi : opt.mnolst())
    {
      for (int i = mnoi.first; i <= mnoi.second; ++i) tsmat.sweep(i);
      fssvec[idx] = tsmat(np(), np()) - sr;
      fdfvec[idx] = mnoi.second - mnoi.first + 1;
      for (int i = mnoi.first; i <= mnoi.second; ++i) tsmat.sweep(i);
      ++idx;
    }
  }
  for (int i = iend; i < np(); ++i) v_smat.sweep(i);
  if ((opt.inolst()).empty() == false)
  {
    mat<double> tsmat(v_smat);
    double sr = tsmat(np(), np());
    for (const auto& inoi : opt.inolst())
    {
      for (int i = inoi.first; i <= inoi.second; ++i) tsmat.sweep(i);
      fssvec[idx] = tsmat(np(), np()) - sr;
      fdfvec[idx] = inoi.second - inoi.first + 1;
      for (int i = inoi.first; i <= inoi.second; ++i) tsmat.sweep(i);
      ++idx;
    }
  }
  cal_vif();
  if (v_smat(np(), np()) < FDeps) v_smat.assign(np(), np(), 0.0);
  v_sb = syy() - v_smat(np(), np());
  if (v_sb < FDeps) v_sb = 0.0;
}

void varclass::cal_vif()        // VIFを計算する
{
  rmat.resize(np(), np());
  rmat.assign(0, 0, xmat(0, 0));
  for (int j = 1; j < np(); ++j)
  { rmat.assign(0, j, xmat(0, j));  rmat.assign(j, 0, xmat(0, j)); }
  int k = 1;
  for (int i = 1; i < np(); ++i)
  {
    rmat.assign(i, i, xmat(i, i));
    for (int j = 1; j < np(); ++j)
    { rmat.assign(i, j, xmat(i, j));  rmat.assign(j, i, xmat(i, j)); }
  }
  rmat.sweep(0);                // 積和行列にする
  smat2rmat(rmat);              // 相関行列にする
  mat<double> tmat(rmat.sub(1, 1, nxix(), nxix()));
  tmat.inv();                   // 相関行列部分を逆行列にする
  for (int i = 0; i < nxix(); ++i)
  { for (int j = 0; j < nxix(); ++j) rmat.assign(1 + i, 1 + j, tmat(i, j)); }
}

class lmaclass                  // lmaメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout, option& opt);
  // 解析結果のタイトルとデータ識別情報を出力する
  int set_itemix(idd& tid, option& opt, int ni) const;
  // 交互作用項目名をセットする
  void put_lma(ostream& fout, option& opt, const varclass& var) const;
  // 線形モデルによる分析結果を出力する
  void put_lmatab(ostream& fout, option& opt, const varclass& var) const;
  // 線形モデルによる分析分析表を出力する
  void put_bci(ostream& fout, option& opt, const varclass& var) const;
  // 偏回帰係数の信頼区間を出力する
public:
  lmaclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~lmaclass()                   // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み累積する
  void put(option& opt);                     // 集計結果を出力する
};

void lmaclass::put_resid(ostream& fout, option& opt)
// 解析結果のタイトルとデータ識別情報を出力する
{
  puttitle(fout, Title);
  fout << "\n";
  id.put_title(fout);
  id.put_stra(fout);
  idd tid(id);
  tid.set(id.title(), id.ni() + (opt.ixlst()).size()); // 項目数を変更する
  set_itemix(tid, opt, id.ni());
  fout << "目的変数y   :" << tid.item(0) << "\n";
  for (size_t i = 1; i < tid.ni(); ++i)
  { fout << "説明変数x" << setw(3) << i << ":" << tid.item(i) << "\n"; }
}

int lmaclass::set_itemix(idd& tid, option& opt, int ni) const
// 交互作用項目名をセットする
{
  ostringstream ostr;
  int j = ni;
  for (const auto& ixi : opt.ixlst())
  {
    ostr.str("");
    if (opt.ixc())              // 偏差変換
    {
      if (ixi.first < ni) ostr << tid.item(ixi.first) << "偏差 x ";
      else ostr << tid.item(ixi.first) << " x ";
      if (ixi.second < ni) ostr << tid.item(ixi.second) << "偏差";
      else ostr << tid.item(ixi.second);
    }
    else
    { ostr << tid.item(ixi.first) << " x " << tid.item(ixi.second); }
    tid.assign_item(j, ostr.str());
    ++j;
  }
  return ni - j;                // セットした交互作用項目名数を返す
}

void lmaclass::put_lma(ostream& fout, option& opt, const varclass& var) const
// 線形モデルによる分析結果を出力する
{
  fout << "\n・重回帰分析結果\n";
  fout << "                     標準                                          有意確率\n";
  fout << "変数  偏回帰係数  偏回帰係数  偏相関係数        VIF       偏F値       p値\n";
  constexpr int LENline = 78;   // 区切り線の長さ
  putline(fout, LENline);
  double bi = var.b(0), vbi = var.vb(0);
  double fi = (vbi > 0.0) ? bi * bi / vbi : 0.0;
  ldbl pr = ftop(fi, 1, var.dfr());
  fout << "切片" << setw(12) << bi << setw(35) << "";
  fout << setw(12) << fi << setw(12) << pr << pind(pr) << "\n";
  for (int i = 1; i < var.np(); ++i)
  {
    bi = var.b(i);  vbi = var.vb(i);
    fi = (vbi > 0.0) ? bi * bi / vbi : 0.0;
    pr = ftop(fi, 1, var.dfr());
    fout << "x" << setw(3) << i << setw(12) << bi << setw(12) << var.bst(i);
    fout << setw(12) << var.rp(i) << setw(11) << var.vif(i) << setw(12) << fi;
    fout << setw(12) << pr << pind(pr) << "\n";
  }
  putline(fout, LENline);
  put_bci(fout, opt, var);
  fout << "重寄与率　:R2乗=" << var.r2();
  fout << "  自由度調整済重寄与率　:R'2乗=" << var.r2a() << "\n";
  fout << "重相関係数:R   =" << sqrt(var.r2()) << "  自由度調整済重相関係数:R'   =" << sqrt(var.r2a()) << "\n";
  fout << "AIC(赤池の情報量基準)=" << var.aic() << endl;
}

void lmaclass::put_bci(ostream& fout, option& opt, const varclass& var) const
// 偏回帰係数の信頼区間を出力する
{
  fout << "\n変数    偏回帰係数      標準誤差" << setw(5) << opt.cc() * 100.0;
  fout << "%絶対精度  信頼区間下限  信頼区間上限\n";
  constexpr int LENline = 74;   // 区切り線の長さ
  putline(fout, LENline);
  double bi = var.b(0), vbi = var.vb(0);
  double ci = calmcci(mcfis, opt.cc(), 1, 1, var.dfr(), 1, vbi);
  fout << "切片" << setw(14) << bi << setw(14) << sqrt(vbi);
  fout << setw(14) << ci << setw(14) << bi - ci;
  fout << setw(14) << bi + ci << "\n";
  for (int i = 1; i < var.np(); ++i)
  {
    bi = var.b(i);  vbi = var.vb(i);
    ci = calmcci(mcfis, opt.cc(), 1, 1, var.dfr(), 1, vbi);
    fout << "x" << setw(3) << i << setw(14) << bi;
    fout << setw(14) << sqrt(vbi) << setw(14) << ci;
    fout << setw(14) << bi - ci << setw(14) << bi + ci << "\n";
  }
  putline(fout, LENline);
}

void lmaclass::put_lmatab(ostream& fout, option& opt, const varclass& var) const
// 線形モデルによる分析分析表を出力する
{
  constexpr int LENline = 78;   // 区切り線の長さ
  string tstr("分散分析表(ANOVA table)");
  int npad = tstr.size() - strclen(tstr, CHARwide);
  fout << "\n";
  fout << setw(50 + npad) << tstr;
  fout << "\n";
  fout << "要因              平方和    自由度    平均平方和         F値   有意確率p値\n";
  putline(fout, LENline);
  double v = 0.0, f = 0.0, pr = 0.0;
  double sb = var.sb();
  int dfb = var.dfb();
  if ((opt.mnolst()).empty() && (opt.inolst()).empty()) // 説明変数の番号指定なし
  {
    calvfp(v, f, pr, sb, dfb, var.sr(), var.dfr());
    tstr = "回帰";
    npad = tstr.size() - strclen(tstr, CHARwide);
    fout << left << setw(10 + npad) << tstr << right << setw(14) << sb;
    fout << setw(10) << dfb << setw(14) << v << setw(12) << f;
    fout << setw(14) << pr << pind(pr) << "\n";
  }
  else
  {
    int idx = 0;
    for (const auto& mnoi : opt.mnolst())
    {
      double ss = var.ss(idx);
      int df = var.df(idx);
      calvfp(v, f, pr, ss, df, var.sr(), var.dfr());
      ostringstream ostr;
      ostr << "M:x" << mnoi.first << "-" << mnoi.second;
      fout << left << setw(10) << ostr.str() << right << setw(14) << ss;
      fout << setw(10) << df << setw(14) << v << setw(12) << f;
      fout << setw(14) << pr << pind(pr) << "\n";
      sb -= ss;  dfb -= df;  ++idx;
    }
    for (const auto& inoi : opt.inolst())
    {
      double ss = var.ss(idx);
      int df = var.df(idx);
      calvfp(v, f, pr, ss, df, var.sr(), var.dfr());
      ostringstream ostr;
      ostr << "I:x" << inoi.first << "-" << inoi.second;
      fout << left << setw(10) << ostr.str() << right << setw(14) << ss;
      fout << setw(10) << df << setw(14) << v << setw(12) << f;
      fout << setw(14) << pr << pind(pr) << "\n";
      sb -= ss;  dfb -= df;  ++idx;
    }
    if (dfb > 0)
    {
      calvfp(v, f, pr, sb, dfb, var.sr(), var.dfr());
      tstr = "その他";
      npad = tstr.size() - strclen(tstr, CHARwide);
      fout << left << setw(10 + npad) << tstr << right << setw(14) << sb;
      fout << setw(10) << dfb << setw(14) << v << setw(12) << f;
      fout << setw(14) << pr << pind(pr) << "\n";
    }
  }
  tstr = "残差";
  npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << var.sr();
  fout << setw(10) << var.dfr() << setw(14) << var.vr() << "\n";
  putline(fout, LENline);
  tstr = "全体";
  npad = tstr.size() - strclen(tstr, CHARwide);
  fout << left << setw(10 + npad) << tstr << right << setw(14) << var.syy();
  fout << setw(10) << var.dft() << endl;
}

int lmaclass::get(const char* c_fname, option& opt)
// データを読み込み保存・累積する
{
  if (c_fname == nullptr) fname.erase();
  else fname = c_fname;
  int nc = fgetdata(fname, id, datalst);
  if (nc < 0) throw runtime_error("実行時エラー！");
  if (opt.check(id.ni()) != 0) throw option_error("実行時オプションエラー！");
  return nc;
}

void lmaclass::put(option& opt) // 集計結果を出力する
{
  varclass var(id.ni(), opt);   // 変数オブジェクト
  var.cal(opt, datalst);        // 各種の計算実行
  if (opt.res("res"))
  {
    put_resid(cout, opt);
    if (opt.res("stat"))
    {
      putmsdvec(cout, var.n(), var.mevec(), var.vrvec(), var.nxix(), 1);
      putmat(cout, "・相関行列(correlation matrix)", matr, var.n(), var.smat(), var.nxix(), 1);
    }
    if (var.resvar() == reserr)
    {
      cout << "※重回帰式計算不可能！\n\n" << endl;
      return;
    }
  }
  if (var.resvar() == reserr) return;
  var.callma(opt);
  if (opt.res("res"))
  {
    if (opt.res("reg")) put_lma(cout, opt, var);
    if (opt.res("atab")) put_lmatab(cout, opt, var);
    cout << "\n" << endl;
  }
}

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