// wmra.C Version 7.45  Copyleft (c) 杉本典夫
// Time-stamp: <2026-04-27 15:44:08 tac>
//
// 最後の項目を重みw、次の項目を目的変数Y、以後の項目を説明変数Xとした
// 重み付け重回帰分析
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -W wmra wmra.C
//               -I/usr/local/include -L/usr/local/lib -ldans -lfun

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

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

const vecstr USAGE =            // プログラム使用法表示用
{
  "wmra  Version 7.45          Copyleft(Vo\\oV) 2026 杉本典夫",
  "書式:wmra [<オプション> …]  [<データファイル> …]",
  "機能:最初の項目を重みw、次の項目を目的変数Y、以後の項目を説明変数Xとした",
  "     重み付け重回帰分析を行います。",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定します。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -dnsname[=<出力ファイル名>]:DANSデータを出力するファイル名を指定する。",
  "    オプション省略時:DANSデータ出力せず  -dnsname:-dnsname=wmra_dns.d",
  "  -f=<ファイル>:オプションを指定したファイルから読み込みます。",
  "  -fin=<取り込み基準>:変数の取り込み基準Finを指定します。(Fin≧0)",
  "  -fout=<追い出し基準>:変数の追い出し基準Foutを指定します。(0≦Fout≦Fin)",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示します。(-h:-h=man)",
  "  -ix=<説明変数の開始番号>[-<終了番号>]x<説明変数の開始番号>[-<終了番号>][,…]",
  "    交互作用項目を作成する2つの説明変数の番号を指定します。",
  "  -ixc[={-|+}]:交互作用項目の説明変数を中心化するかどうか指定します。",
  "    -:中心化しない(オプション省略時)  +:中心化する  -ixc:ixc=+",
  "  -ngp={est|line|pl}[,…]:Ngraph用データを指定します。",
  "    est:Y推定値/Y補正値  line:回帰直線/偏回帰直線  pl:予測限界",
  "    オプション省略時:全て指定",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データを出力するファイル名を指定する。",
  "    オプション省略時:Ngraph用データ出力せず  -ngpname:-ngpname=wmra_ngp.dat",
  "  -pin=<取り込み基準>:変数の取り込み基準Pinを指定します。(0<Pin≦1)",
  "  -pout=<追い出し基準>:変数の追い出し基準Poutを指定します。(1≧Pout≧Pin)",
  "  -res=[-][{def|all|allx|atab|selx|stat}][,…]:",
  "      解析結果の出力項目を指定します。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    allx:全変数選択結果  atab:分散分析表  selx:変数選択結果  stat:基礎統計量",
  "  -tstat[={c|f}]:検定統計量の種類を指定します。(-tstat:tstat=f)",
  "  -x<説明変数番号>=<説明変数の値>:X軸にしない説明変数の値を指定します。",
  "  -xno=<説明変数の番号>:X軸にする説明変数の番号を指定します。",
  "  -xrange=<最小値>-<最小値>:X軸の最小値と最大値を指定します。"
};

const string Title = "=== 重み付け重回帰分析(weighted multiple regression analysis) ==="; // 表題
const string DNSNAMEdef = "wmra_dns.d"; // 省略時DANSファイル名
const string NGPNAMEdef = "wmra_ngp.dat"; // 省略時Ngraphファイル名
constexpr double FIOdef = 2.0;  // 省略時F値の取り込み・追い出し基準
constexpr double PIOdef = 0.1;  // 省略時p値の取り込み・追い出し基準

enum tstattype                  // 検定統計量の種類
{
  tstatnon,                     // 指定なし
  tstatc,                       // χ2乗値
  tstatf,                       // F値
  tstatetc                      // その他
};

class option                    // オプションクラス
{
private:
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  string v_dnsname;             // DANSファイル名
  pairdbl flevel;               // F値の取り込み基準と追い出し基準
  pairdbl plevel;               // P値の取り込み基準と追い出し基準
  lstr ixstrlst;                // 交互作用項目の説明変数リスト
  lino v_ixlst;                 // 交互作用項目の説明変数ペアリスト
  bool v_ixc;                   // 交互作用項目の中心化フラッグ
  string v_ngpname;             // Ngraphファイル名
  lstr ngplst;                  // Ngraph出力リスト
  mapsb ngpmap;                 // Ngraph出力の有無マップ
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
  tstattype v_tstat;            // 検定統計量の種類
  lstr xlst;                    // 説明変数リスト
  vecdbl xvec;                  // 説明変数ベクトル
  int v_xno;                    // X軸にする説明変数の番号
  pairdbl xrange;               // X軸の最小値と最大値
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 string& dnsname() const // DANSデータファイル名を返す
  { return v_dnsname; }
  double fin() const            // F値の取り込み基準を返す
  { return flevel.first; }
  double fout() const           // F値の追い出し基準を返す
  { return flevel.second; }
  const lino& ixlst()    // 交互作用項目の説明変数番号ペアリストを返す
  { return v_ixlst; }
  bool ixc()                   // 交互作用項目の中心化フラッグを返す
  { return v_ixc; }
  double pin() const            // P値の取り込み基準を返す
  { return plevel.first; }
  double pout() const           // P値の追い出し基準を返す
  { return plevel.second; }
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  bool ngp(const string& code); // Ngraph出力フラッグを返す
  bool res(const string& code); // 分析結果出力フラッグを返す
  tstattype tstat() const       // 検定統計量の種類を返す
  { return v_tstat; }
  double x(size_t xno) const    // 説明変数の値を返す
  { return (0 < xno && xno <= xvec.size()) ? xvec[xno - 1] : BLANKf; }
  int xno() const               // X軸にする説明変数の番号を返す
  { return v_xno; }
  double xmin() const           // X軸の最小値を返す
  { return xrange.first; }
  double xmax() const           // X軸の最大値を返す
  { return xrange.second; }
  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 ni);           // -ixオプションのセット
  int set_x(int nx);            // -xオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_cc(CCdef), v_dnsname(), flevel(FIOdef, FIOdef), ixstrlst(), v_ixlst(), v_ixc(false), plevel(1.0, 1.0), v_ngpname(), ngplst(), ngpmap(), reslst(), resmap(), v_tstat(tstatf), xlst(), xvec(), v_xno(0), xrange(BLANKf, BLANKf)
// コンストラクタ
{
  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 ngpitr = ngpmap.find(code);
  if (ngpitr == ngpmap.end()) return false; // code無し！
  return ngpmap[code];                      // true/falseを返す
}

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] == "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] == "fin")     // F値の取り込み基準指定
  {
    flevel.first = opt[1].empty() ? FIOdef : strtodf(opt[1]);
    if (flevel.first < 0.0) err = badopt(pname(), optstr, USAGE, 0.0, FDmax);
    if (flevel.second > flevel.first) flevel.second = flevel.first;
  }
  else if (opt[0] == "fout")    // F値の追い出し基準指定
  {
    flevel.second = opt[1].empty() ? flevel.first : strtodf(opt[1]);
    if (flevel.second < 0.0 || flevel.first < flevel.second) err = badopt(pname(), optstr, USAGE, 0.0, flevel.first);
  }
  else if (opt[0] == "h")       // オンラインヘルプ表示
  { err = help(pname(), USAGE, opt[1]); }
  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] == "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] == "pin")     // p値の取り込み基準指定
  {
    plevel.first = opt[1].empty() ? PIOdef : strtodf(opt[1]);
    if (plevel.first <= 0.0 || 1.0 < plevel.first) err = badopt(pname(), optstr, USAGE, 0.0, 1.0);
  }
  else if (opt[0] == "pout")    // p値の追い出し基準指定
  {
    plevel.second = opt[1].empty() ? plevel.first : strtodf(opt[1]);
    if (plevel.second < plevel.first || 1.0 < plevel.second) err = badopt(pname(), optstr, USAGE, plevel.first, 1.0);
  }
  else if (opt[0] == "res")     // 分析結果の出力指定
  {
    if (strsplit(reslst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "tstat")   // 検定統計量の種類指定
  {
    if (opt[1].empty() || opt[1] == "f") v_tstat = tstatf;
    else if (opt[1] == "c") v_tstat = tstatc;
    else err = badopt(pname(), optstr, USAGE);
  }
  else if (opt[0] == "xno")     // X軸にする説明変数の番号指定
  {
    char* ptr = nullptr;
    v_xno = strtol(opt[1].c_str(), &ptr, 0);
    if (v_xno < 1) err = badopt(pname(), optstr, USAGE, 1, INT_MAX);
  }
  else if (opt[0] == "xrange")  // X軸の最小値と最大値の指定
  {
    if (getrange(xrange.first, xrange.second, opt[1]) == 0) // 範囲指定無
    { err = badopt(pname(), optstr, USAGE, -FDmax, FDmax); }
    else if (xrange.first >= xrange.second) // 範囲指定不適切
    { err = badopt(pname(), optstr, USAGE, -FDmax, FDmax); }
  }
  else if (opt[0][0] == 'x')    // 説明変数の値指定
  {
    if (strsplit(xlst, optstr, 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 < 3) return disperr(pname(), "項目数3つ以上必要です！");
  int nx = (int)ni - 1;         // 交互作用項目を含まない説明変数の数
  if (nx < v_xno) return disperr(pname(), "X軸にする説明変数の番号が間違っています！");
  if (set_ngp() < 0) return -2; // -ngpオプションのセット
  if (set_res() < 0) return -3; // -resオプションのセット
  int nxix = set_ix(nx);        // -ixオプションのセット
  if (nxix < 0) return -4;      // 交互作用項目を含んだ説明変数の数
  if (set_x(nx) < 0) return -5; // -xオプションのセット
  return 0;                     // エラーなし
}

int option::set_ngp()           // -ngpオプションのセット
{
  ngpmap["est"] = ngpmap["line"] = ngpmap["pl"] = true;
  if (ngplst.size() > 0)
  {
    if (v_ngpname == "") return disperr(pname(), "-ngpnameオプションの指定が必要！");
    ngpmap["est"] = ngpmap["line"] = ngpmap["pl"] = false;
    for (const string& ngpstr : ngplst)
    {
      mapsb_citr ngpitr = ngpmap.find(ngpstr);
      if (ngpitr == ngpmap.end())
      { return disperr(pname(), "-ngp=" + ngpstr + "オプション指定が不適切！"); }
      else
      { ngpmap[ngpstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["allx"] = resmap["atab"] = resmap["selx"] = resmap["stat"] = true;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["allx"] = resmap["atab"] = resmap["selx"] = 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["allx"] = resmap["atab"] = resmap["selx"] = 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_x(int nx)       // -xオプションのセット
{
  if (xlst.empty()) return 0;   // 説明変数指定なし
  if ((int)(xlst.size()) > nx) return disperr(pname(), "説明変数の値が多すぎます！");
  xvec.resize(nx, BLANKf);
  for (const string& xstr : xlst)
  {
    string opt[2];
    if (getopt(opt, xstr, OPTdev) > 0)
    {
      char* ptr = nullptr;
      int xno = strtol(opt[0].substr(1).c_str(), &ptr, 0);
      if (xno < 1 || nx < xno || xno == v_xno) return disperr(pname(), "説明変数の番号が間違っています→", xstr);
      double fdata = strtodf(opt[1]);
      if (fdata == BLANKf) return disperr(pname(), "説明変数の値が間違っています→", xstr);
      xvec[xno - 1] = fdata;
    }
  }
  return (int)(xlst.size());    // 説明変数指定の数を返す
}

class varclass                 // 変数クラス
{
private:
  mat<double> v_smat;           // 積和行列
  mat<double> xmat;             // 単純積和行列
  mat<double> rmat;             // 相関行列
  vector<statv<double>> wyxstatvec; // wとyとxの基礎統計量
  vecdbl v_mevec;               // 平均値ベクトル
  vecdbl v_vrvec;               // 分散ベクトル
  vecdbl ssvec;                 // 積和ベクトル
  vecint checkvec;              // 変数選択用ベクトル
  int v_ni;                     // 項目数
  int v_nxix;                   // 交互作用項目を含んだ説明変数の数
  int v_n;                      // 例数
  restype v_resvar;             // 処理結果の種類
  int set_wyxstat(vector<statv<double>>& wyxstatvec, ldata& datalst) const;
  // wとyとxの統計量ベクトルをセットする
  bool set_fdata(vecdbl& fdatavec, option& opt, const vector<statv<double>>& wyxstatvec, datac& cdata) const;
  // 交互作用項目を含むデータをセットする
  int sel_var(int in, double fin, double fout, double pin, double pout);
  // 変数増減法による変数の選択
  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; }
  double me(int idx) const      // 重み付け平均値を返す
  { return (0 <= idx && idx <= np()) ? v_mevec[idx] : 0.0; }
  double vr(int idx) const      // 重み付け分散を返す
  { return (0 <= idx && idx <= np()) ? v_vrvec[idx] : 0.0; }
  int wyxn(int idx) const       // 重み付け例数を返す
  { return (0 <= idx && idx <= np()) ? wyxstatvec[idx].n() : 0.0; }
  double wyxme(int idx) const   // 重み無し平均値を返す
  { return (0 <= idx && idx <= np()) ? wyxstatvec[idx].mean() : 0.0; }
  double wyxsum(int idx) const  // 重み無し合計を返す
  { return (0 <= idx && idx <= np()) ? wyxstatvec[idx].sum() : 0.0; }
  double wyxvr(int idx) const   // 重み無し分散を返す
  { return (0 <= idx && idx <= np()) ? wyxstatvec[idx].vr() : 0.0; }
  int check(int idx) const      // 変数選択フラッグを返す
  { return (0 <= idx && idx < np()) ? checkvec[idx] : 0; }
  int ni() const                // 項目数を返す
  { return v_ni; }
  int nx() const                // 説明変数の数を返す
  { return ni() - 2; }
  int nxix() const              // 交互作用項目を含んだ説明変数の数を返す
  { return v_nxix; }
  int np() const                // パラメータ数を返す
  { return nxix() + 1; }
  int dim() const               // 行列の次元数を返す
  { return nxix() + 2; }
  int n() const                 // 例数を返す
  { return v_n; }
  int dft() const               // 全自由度を返す
  { return n() - 1; }
  int dfb() const               // 回帰自由度を返す
  { return checkvec[0]; }
  int dfr() const               // 残差自由度を返す
  { return dft() - dfb(); }
  double sxx(int idx) const     // 説明変数の平方和を返す
  { return (0 <= idx && idx <= np()) ? ssvec[idx] : 0.0; }
  double syy() const            // 目的変数の平方和を返す
  { return ssvec[np()]; }
  double sr() const             // 残差平方和を返す
  { return v_smat(np(), np()); }
  double vr() const             // 残差分散を返す
  { return (dfr() > 0) ? sr() / (double)(dfr()) : 0.0; }
  double sb() const             // 回帰平方和を返す
  { return syy() - sr(); }
  double b(int idx) const       // 偏回帰係数を返す
  { return (0 <= idx && idx <= np()) ? v_smat(idx, np()) : 0.0; }
  double vcb(int idx) const     // 偏回帰係数の分散係数を返す
  { return (0 <= idx && idx <= np()) ? v_smat(idx, idx) : 0.0; }
  double bst(int idx) const     // 標準偏回帰係数を返す
  { return (0 <= idx && idx <= np() && syy() > 0.0) ? b(idx) * sqrt(sxx(idx) / syy()) : 0.0; }
  double vif(int idx) const     // VIFを返す
  { return (0 < idx && idx <= np()) ? rmat(idx, idx) : 0.0; }
  double rp(int idx) const;     // 偏相関係数を返す
  restype resvar() const        // 処理結果の種類を返す
  { return v_resvar; }
  void cal(option& opt, ldata& datalst); // データを集計して各種統計量を求める
  void stepsel(double fin, double fout, double pin, double pout);
  // 変数増減法による重回帰分析の計算
  void set_yest(datac& estdata, option& opt, double vr, double tcc, datac& cdata) const;
  // 実際のW,Y,X,Y推定値をセットする
  void set_xyest(datac& estdata, option& opt, double vr, double tcc, datac& mdata, datac& cdata) const;
  // 変化説明変数XとY補正値をセットする
  double cal_yest(option& opt, datac& cdata) const; // Yの推定値を計算して返す
  double cal_pi(datac& cdata, double tcc, double vr, option& opt, double yest, statv<double> wyeststat) const;
  // Yの予測限界幅を計算して返す
};

varclass::varclass(size_t ni, option& opt) : v_smat(), xmat(), rmat(), wyxstatvec(), v_mevec(), v_vrvec(), ssvec(), checkvec(), v_ni(ni), v_nxix(0), v_n(0), v_resvar(resnon)
// コンストラクタ
{
  v_nxix = ni - 2 + (opt.ixlst()).size();
  int dim = v_nxix + 2;
  v_smat.resize(dim, dim);  xmat.resize(dim, dim);
  rmat.resize(dim, dim);  wyxstatvec.resize(dim);
  v_mevec.resize(dim, 0.0);  v_vrvec.resize(dim, 0.0);  ssvec.resize(dim, 0.0);
  checkvec.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() * vcb(idx) + bi * bi;
  return (v > 0.0) ? bi / sqrt(v) : 0.0;
}

void varclass::cal(option& opt, ldata& datalst)
// データを集計して各種統計量を求める
{
  v_n = set_wyxstat(wyxstatvec, datalst);
  if (v_n < np())               // 重回帰式計算不可能
  {
    v_resvar = reserr;
    return;
  }
  vecdbl fdatavec(dim(), 0.0);
  for (auto& cdatai : datalst)
  {
    if (set_fdata(fdatavec, opt, wyxstatvec, cdatai))
    {
      double w = fdatavec[0], y = fdatavec[1];
      for (int j = 2; j < dim(); ++j)
      {
        if (j >= ni()) wyxstatvec[j].sigma(fdatavec[j]);
        fdatavec[j - 1] = w * fdatavec[j];
      }
      fdatavec[np()] = w * y;
      for (int i = 0; i < dim(); ++i)
      {
        double fwi = fdatavec[i];
        for (int j = i; j < dim(); ++j) v_smat.add(i, j, fwi * fdatavec[j]);
      }
    }
  }
  xmat = v_smat;
  for (int i = 0; i < np(); ++i)
  {
    for (int j = i + 1; j <= np(); ++j)
    { v_smat.assign(j, i, v_smat(i, j));  xmat.assign(j, i, xmat(i, j)); }
  }
  v_smat.sweep(0);
  v_mevec[0] = wyxstatvec[0].mean();
  v_vrvec[0] = wyxstatvec[0].vr();  ssvec[0] = wyxstatvec[0].ss();
  for (int i = 1; i <= np(); ++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 = (sgn(v_vrvec[np()]) > 0) ? resok : reserr;
}

int varclass::set_wyxstat(vector<statv<double>>& wyxstatvec, ldata& datalst) const
// wとyとxの統計量ベクトルをセットする
{
  vecdbl fdatavec(ni(), 0.0);
  int n = 0;
  for (auto& cdatai : datalst)
  {
    bool blank = false;
    if (cdatai.fdata(0) == BLANKf || sgn(cdatai.fdata(0)) <= 0) blank = true;
    else fdatavec[0] = cdatai.fdata(0); // 重み＞0
    int j = 1;
    while (blank == false && j < ni())
    {
      fdatavec[j] = cdatai.fdata(j);
      if (fdatavec[j] == BLANKf) blank = true;
      else ++j;
    }
    if (blank)
    { cdatai.assign(0, BLANKs); } // 後の計算のためにwをブランクにする
    else
    {
      ++n;
      for (j = 0; j < ni(); ++j) wyxstatvec[j].sigma(fdatavec[j]);
    }
  }
  return n;                     // 例数を返す
}

bool varclass::set_fdata(vecdbl& fdatavec, option& opt, const vector<statv<double>>& wyxstatvec, datac& cdata) const
// 交互作用項目を含むデータをセットする
{
  if (cdata.sdata(0) == BLANKs) // 重みなし＝除外例
  {
    fdatavec.assign(fdatavec.size(), 0.0);
    return false;               // 全データを0にして返す
  }
  for (int j = 0; j < ni(); ++j) fdatavec[j] = cdata.fdata(j);
  if ((opt.ixlst()).empty()) return true; // 交互作用項目なし
  int j = ni();
  for (const auto& ixi : opt.ixlst()) // 交互作用データのセット
  {
    int xidx1 = ixi.first + 1, xidx2 = ixi.second + 1;
    double x1 = fdatavec[xidx1], x2 = fdatavec[xidx2];
    if (opt.ixc())          // 偏差に変換する
    {
      if (xidx1 < ni()) x1 -= wyxstatvec[xidx1].mean();
      if (xidx2 < ni()) x2 -= wyxstatvec[xidx2].mean();
    }
    fdatavec[j] = x1 * x2;
    ++j;
  }
  return true;                  // 全データをセットした
}

void varclass::stepsel(double fin, double fout, double pin, double pout)
// 変数増減法による重回帰分析の計算
{
  if (fin <= 0.0 && pin >= 1.0) // 全変数選択
  {
    for (int i = 1; i <= nxix(); ++i)
    {
      if (checkvec[i] == 0) { v_smat.sweep(i);  checkvec[i] = 1; }
    }
    checkvec[0] = nxix();
  }
  else
  {
    int in = 1, i = 0;
    do {
      in = 1;
      while ((i = sel_var(in, fin, fout, pin, pout)) > 0)
      {
        v_smat.sweep(i);  checkvec[i] = in ? 1 : 0;
        checkvec[0] += in ? 1 : -1;
        in = !in;
      }
    } while (!in);                // 取り込む変数がなくなるまで続ける
  }
  if (checkvec[0] > 0) cal_vif(); // 変数が選択されている時だけVIFを計算する
  if (v_smat(np(), np()) < FDeps) v_smat.assign(np(), np(), 0.0);
}

int varclass::sel_var(int in, double fin, double fout, double pin, double pout)
// 変数増減法による変数の選択
{
  int i = 0;
  double dsel = in ? FDmin : FDmax;
  bool fsel = true;             // F値による変数選択
  if (pin < 1.0) fsel = false;  // 有意確率pによる変数選択
  for (int j = 1; j < np(); ++j)
  {
    if (checkvec[j] != in && fabs(v_smat(j, j)) > FDeps)
    {
      double dj = v_smat(j, np()) * v_smat(j, np()) / v_smat(j, j);
      if ((in && dj > dsel) || (!in && dj < dsel)) { dsel = dj;  i = j; }
    }
  }
  if (i < 1) return i;          // 選択すべき変数なし
  if (in)                       // 変数取り込み
  {
    int df = n() - checkvec[0] - 2;
    double f = (sgndiff(v_smat(np(), np()), dsel) > 0) ? dsel * (double)df / (v_smat(np(), np()) - dsel) : 0.0;
    if (fsel)
    {
      if (f <= fin) i = 0;
    }
    else
    {
      ldbl pr = ftop(f, 1, df);
      if (pr >= (ldbl)pin) i = 0;
    }
  }
  else                          // 変数追い出し
  {
    int df = n() - checkvec[0] - 1;
    double f = dsel * (double)df / v_smat(np(), np());
    if (fsel)
    {
      if (f > fout) i = 0;
    }
    else
    {
      ldbl pr = ftop(f, 1, df);
      if (pr < (ldbl)pout) i = 0;
    }
  }
  return i;                     // 選択すべき変数のインデックスを返す
}

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

void varclass::set_yest(datac& estdata, option& opt, double vr, double tcc, datac& cdata) const
// 実際のW,Y,X,Y推定値をセットする
{
  vecdbl fdatavec(dim(), 0.0);
  if (set_fdata(fdatavec, opt, wyxstatvec, cdata))
  {
    estdata.assign(0, fdatavec[0]);
    estdata.assign(1, fdatavec[1]);
    for (int i = 2; i < dim(); ++i) estdata.assign(i, fdatavec[i]);
    estdata.assign(dim(), cal_yest(opt, cdata));
  }
  else
  { estdata = BLANKs; }
}

void varclass::set_xyest(datac& estdata, option& opt, double vr, double tcc, datac& mdata, datac& cdata) const
// YとXとY補正値をセットする
{
  vecdbl fdatavec(dim(), 0.0);
  int j = 2;
  if (set_fdata(fdatavec, opt, wyxstatvec, cdata))
  {
    double w = fdatavec[0], y = fdatavec[1];
    estdata.assign(0, w);  estdata.assign(1, y);
    double err = y - cal_yest(opt, cdata);
    size_t xidx = 1 + opt.xno();
    mdata.assign(xidx, cdata.fdata(xidx));
    set_fdata(fdatavec, opt, wyxstatvec, mdata);
    for (int i = 2; i < dim(); ++i) estdata.assign(i, fdatavec[i]);
    estdata.assign(dim(), cal_yest(opt, mdata) + err);
  }
  else
  { estdata = BLANKs; }
}

double varclass::cal_yest(option& opt, datac& cdata) const
// Yの推定値を計算して返す
{
  vecdbl fdatavec(dim(), 0.0);
  if (set_fdata(fdatavec, opt, wyxstatvec, cdata) == false) return BLANKf;
  double yest = b(0);
  for (int i = 1; i < np(); ++i)
  {
    if (check(i)) yest += b(i) * fdatavec[i + 1];
  }
  return yest;                  // 推定値を返す
}

double varclass::cal_pi(datac& cdata, double tcc, double vr, option& opt, double yest, statv<double> wyeststat) const
// Yの予測限界幅を計算して返す
{
  vecdbl fdatavec(dim(), 0.0);
  if (set_fdata(fdatavec, opt, wyxstatvec, cdata) == false) return BLANKf;
  double w = fdatavec[0], vpi = 0.0;
  if (opt.xno() > 0)
  {
    vec<double> xv(dim());
    for (int i = 1; i < np(); ++i)
    {
      if (check(i)) xv.assign(i, w * fdatavec[i + 1]);
      else xv.assign(i, 0.0);
    }
    xv.assign(0, w);  xv.assign(np(), 0.0);
    vec<double> tv(xv.trn() * v_smat);
    vpi = 1.0 + inp(tv, xv);
  }
  else
  {
    double xd = w * (yest - wyxme(1));
    xd = (wyeststat.ss() > 0.0) ? xd * xd / wyeststat.ss() : 0.0;
    vpi = 1.0 + wyxme(0) * wyxme(0) / wyxsum(0) + xd;
  }
  if (vpi < FDeps) vpi = 0.0;
  return tcc * sqrt(vr * vpi) / wyxme(0); // Yの予測限界を返す
}

class wmraclass                // wmraメインクラス
{
private:
  string fname;                 // ファイル名
  idd id;                       // データ識別情報
  ldata datalst;                // データリスト
  void put_resid(ostream& fout, option& opt);
  // 解析結果のタイトルとデータ識別情報を出力する
  void put_stat(ostream& fout, const varclass& var);
  // 変数の基礎統計量を出力する
  void put_regtab(ostream& fout, option& opt, double fin, double flout, double pin, double pout, const varclass& var);
  // 重回帰分析結果を出力する
  void put_bci(ostream& fout, option& opt, const varclass& var) const;
  // 偏回帰係数の信頼区間を出力する
  void put_catab(ostream& fout, double sbw, int dfb, double srw, int dfr, double syyw, int dft) const;
  // χ2乗値による分散分析表を出力する
  void put_dns(option& opt, const varclass& var);
  // Yの推定値をデータファイル形式で出力する
  void put_dnsid(ofstream& fout, option& opt, const varclass& var);
  // データファイル形式のデータ識別情報を出力する
  int set_itemix(idd& tid, option& opt, size_t ni) const;
  // 交互作用項目名をセットする
  void plot_regngp(option& opt, const varclass& var);
  // Ngraph用データ出力
  void plot_ngpid(ofstream& fout, option& opt, const varclass& var);
  // Ngraph用にデータ識別情報を出力する
  void plot_est(ofstream& fout, option& opt, const varclass& var);
  // Yの推定値をNgraph用に出力する
  double cal_estimate(option& opt, const mat<double>& smat, const vector<statv<double>>& wyxstatvec, const vecint& checkvec, datac& cdata) const;
  // Yの推定値を計算して返す
  void plot_line(ofstream& fout, option& opt, const varclass& var);
  // 回帰直線をNgraph用に出力する
  double search_xrange(double& xmin, double& xmax, option& opt, const varclass& var);
  // x軸の最小値と最大値を検索して範囲を返す
  void plot_pl(ofstream& fout, option& opt, const varclass& var);
  // 予測限界をNgraph用に出力する
  void plot_pllu(ofstream& fout, option& opt, datac& cdata, const varclass& var, double t, double vr, double xmin, double xmax, statv<double>& wyeststat, bool cup) const;
  // 下限／上限の予測限界をNgraph用に出力する
public:
  wmraclass() : fname(), id(), datalst() // コンストラクタ
  {}
  ~wmraclass()                  // デストラクタ
  {}
  int get(const char* c_fname, option& opt); // データを読み込み累積する
  void put(option& opt);                     // 集計結果を出力する
};

void wmraclass::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 << "重みw       :" << tid.item(0) << "\n";
  fout << "目的変数y   :" << tid.item(1) << "\n";
  for (size_t i = 2; i < tid.ni(); ++i)
  { fout << "説明変数x" << setw(3) << i - 1 << ":" << tid.item(i) << "\n"; }
}

void wmraclass::put_stat(ostream& fout, const varclass& var)
// 変数の基礎統計量を出力する
{
  constexpr int LENline = 72;
  fout << "\n・各変数の基礎統計量と重み付け基礎統計量\n";
  putline(fout, LENline);
  fout << "w   :";
  putmsd(fout, var.wyxn(0), var.wyxme(0), var.wyxvr(0));
  fout << "y   :";
  putmsd(fout, var.wyxn(1), var.wyxme(1), var.wyxvr(1));
  for (int j = 2; j <= var.np(); ++j)
  {
    fout << "x" << setw(3) << j - 1 << ":";
    putmsd(fout, var.wyxn(j), var.wyxme(j), var.wyxvr(j));
  }
  putline(fout, LENline);
  fout << "yw  :";
  putmsd(fout, var.n(), var.me(var.np()), var.vr(var.np()));
  for (int j = 1; j < var.np(); ++j)
  {
    fout << "xw" << setw(2) << j << ":";
    putmsd(fout, var.n(), var.me(j), var.vr(j));
  }
  putline(fout, LENline);
}

void wmraclass::put_regtab(ostream& fout, option& opt, double fin, double flout, double pin, double pout, const varclass& var)
// 重回帰分析結果を出力する
{
  if (fin <= 0.0 && pin >= 1.0)
  { fout << "\n・全変数を選択した結果\n"; }
  else
  {
    fout << "\n・前進的変数増減法(stepwise forward selection method)による変数選択結果\n";
    if (pin < 1.0) fout << "  取り込み基準:p値＜Pin=" << pin << "  追い出し基準:p値≧Pout=" << pout << "\n";
    else fout << "  取り込み基準:偏F値＞Fin=" << fin << "  追い出し基準:偏F値≦Fout=" << flout << "\n";
  }
  fout << "                     標準                                          有意確率\n";
  fout << "変数  偏回帰係数  偏回帰係数  偏相関係数        VIF";
  if (opt.tstat() == tstatc) fout << "      χ2乗値";
  else fout << "       偏F値";
  fout << "       p値\n";
  constexpr int LENline = 78;   // 区切り線の長さ
  putline(fout, LENline);
  double vr = var.vr();
  if (opt.tstat() == tstatc) vr = 1.0; // 検定統計量としてχ2乗値を使う
  double bi = var.b(0), vbi = vr * var.vcb(0);
  double fi = (vbi > 0.0) ? bi * bi / vbi : 0.0;
  double pr = ftop(fi, 1, var.dfr());
  fout << "切片" << setw(12) << bi << setw(35) << "" << setw(12) << fi;
  fout << setw(12) << pr << pind(pr) << "\n";
  int k = 1;
  for (int i = 1; i < var.np(); ++i)
  {
    if (var.check(i))
    {
      bi = var.b(i);  vbi = vr * var.vcb(i);
      fi = (vbi > 0.0) ? bi * bi / vbi : 0.0;
      pr = (opt.tstat() == tstatc) ? ctop(fi, 1) : 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(k) << setw(12) << fi;
      fout << setw(12) << pr << pind(pr) << "\n";
      ++k;
    }
  }
  putline(fout, LENline);
  put_bci(fout, opt, var);
  double r2 = (var.syy() > 0.0) ? var.sb() / var.syy() : 0.0;
  double r2a = (var.syy() > 0.0) ? fabs(1.0 - vr * (double)(var.dft()) / var.syy()) : 0.0;
  double aic = (var.sr() > 0.0) ? (double)var.n() * log(var.sr() / (double)var.n()) + (double)(2 * (var.dfb() + 2)) : 0.0;
  fout << "重寄与率　:R2乗=" << r2;
  fout << "  自由度調整済重寄与率　:R'2乗=" << r2a << "\n";
  fout << "重相関係数:R   =" << sqrt(r2) << "  自由度調整済重相関係数:R'   =" << sqrt(r2a) << "\n";
  fout << "AIC(赤池の情報量基準)=" << aic << "\n\n";
  if (opt.res("atab"))
  {
    if (opt.tstat() == tstatc) put_catab(fout, var.sb(), var.dfb(), var.sr(), var.dfr(), var.syy(), var.dft());
    else putatab(fout, "回帰", var.sb(), var.dfb(), var.sr(), var.dfr(), var.syy(), var.dft());
  }
}

void wmraclass::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 vr = var.vr();
  if (opt.tstat() == tstatc) vr = 1.0; // 検定統計量としてχ2乗値を使う
  double bi = var.b(0), vbi = vr * var.vcb(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)
  {
    if (var.check(i))
    {
      bi = var.b(i);  vbi = vr * var.vcb(i);
      ci = (opt.tstat() == tstatc) ? calmcci(mcfis, opt.cc(), 1, 1, 0, 1, vbi) :  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 wmraclass::put_catab(ostream& fout, double sbw, int dfb, double srw, int dfr, double syyw, int dft) const
// χ2乗値による分散分析表を出力する
{
  const string Head = "分散分析表(ANOVA table)";
  int w = strclen(Head, CHARwide);
  int npad = Head.size() - w;
  fout << setw((68 + w) / 2 + npad) << Head;
  fout << "\n要因                平方和    自由度        χ2乗値   有意確率p値\n";
  constexpr int LENline = 68;
  putline(fout, LENline);
  ldbl pr = ctop(sbw, dfb);
  fout << "回帰" << setw(22) << sbw << setw(10) << dfb << setw(14) << sbw;
  fout << setw(14) << pr << pind(pr) << "\n";
  pr = ctop(srw, dfr);
  fout << "残差" << setw(22) << srw << setw(10) << dfr << setw(14) << srw;
  fout << setw(14) << pr << pind(pr) << "\n";
  putline(fout, LENline);
  pr = ctop(syyw, dft);
  fout << "全体" << setw(22) << syyw << setw(10) << dft << setw(14) << syyw;
  fout << setw(14) << pr << pind(pr) << endl;
}

void wmraclass::put_dns(option& opt, const varclass& var)
// Y推定値をデータファイル形式で出力する
{
  ofstream fout;
  if (ofopen(fout, (opt.dnsname()).c_str()) == false) return; // ファイルオープン失敗
  put_dnsid(fout, opt, var);
  datac estdata(var.np() + 2); // WとYとXとY推定値
  double vr = var.vr();
  double tcc = ptot(1.0 - opt.cc(), var.dfr());
  if (opt.tstat() == tstatc) // 検定統計量としてχ2乗値を使う
  {
    vr = 1.0;
    tcc = ptoz(1.0 - opt.cc());
  }
  if (opt.xno() > 0)            // X軸にする説明変数が指定されている時
  {
    datac mdata(var.ni());
    mdata.assign(0, var.me(0));  mdata.assign(1, var.me(var.np()));
    for (int i = 2; i < var.ni(); ++i)
    {
      int xno = i - 1;
      mdata.assign(i, (opt.x(xno) == BLANKf) ? var.me(xno) : opt.x(xno));
    }
    for (auto& cdatai : datalst)
    {
      var.set_xyest(estdata, opt, vr, tcc, mdata, cdatai);
      estdata.put(fout);        // WとYとX(固定X含む)とY補正値を出力
    }
  }
  else                          // X軸にする説明変数がない時
  {
    for (auto& cdatai : datalst)
    {
      var.set_yest(estdata, opt, vr, tcc, cdatai);
      estdata.put(fout);        // WとYとXとY推定値を出力
    }
  }
}

void wmraclass::put_dnsid(ofstream& fout, option& opt, const varclass& var)
// データファイル形式のデータ識別情報を出力する
{
  idd tid(id);
  tid.set(id.title(), var.np() + 2); // 項目数だけ変更する
  ostringstream ostr;
  if (opt.pin() < 1.0) ostr << "取り込み基準:p値＜Pin=" << opt.pin() << "  追い出し基準:p値≧Pout=" << opt.pout();
  else ostr << "取り込み基準:偏F値＞Fin=" << opt.fin() << "  追い出し基準:偏F値≦Fout=" << opt.fout();
  tid.append_stra(ostr.str(), attnon, avalnon); // 層別情報として追加する
  tid.assign_item(0, id.item(0));  tid.assign_item(1, id.item(1));
  for (int i = 2; i < var.ni(); ++i)
  {
    int xno = i - 1;
    ostr.str("");
    ostr << id.item(i);
    if (opt.xno() < 1)          // X軸にする説明変数なし
    { ostr << " 重み付け平均値=" << var.me(xno); }
    else
    {
      if (xno == opt.xno()) ostr << " X軸にする説明変数";
      else if (opt.x(xno) == BLANKf) ostr << " 重み付け平均値=" << var.me(xno);
      else ostr << " 指定値=" << opt.x(xno);
    }
    if (var.check(xno)) ostr << "(選択)";
    else ostr << "(非選択)";
    tid.assign_item(i, ostr.str());
  }
  if ((opt.ixlst()).size() > 0) // 交互作用項目あり
  {
    idd txid(id);
    txid.set(id.title(), var.ni() + (opt.ixlst()).size());
    set_itemix(txid, opt, var.ni()); // 交互作用項目を追加する
    for (int i = var.ni(); i < (int)(txid.ni()); ++i)
    {
      int xno = i - 1;
      ostr << txid.item(i);
      if (var.check(xno)) ostr << "(選択)";
      else ostr << "(非選択)";
      tid.assign_item(i, ostr.str());
    }
  }
  ostr.str("");
  ostr << id.item(1);
  if (opt.xno() > 0) ostr << "補正値";
  else ostr << "推定値";
  tid.assign_item(tid.ni() - 1, ostr.str());
  tid.put(fout, formtag);
}

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

void wmraclass::plot_regngp(option& opt, const varclass& var)
// Ngraph用データ出力
{
  ofstream fout;
  if (ofopen(fout, (opt.ngpname()).c_str()) == false) return; // ファイルオープン失敗
  bool delml = false;
  if (opt.ngp("est"))            // Yの推測値/補正値
  {
    fout << "-s1 -dmark,00 -o300\n";
    plot_ngpid(fout, opt, var);
    plot_est(fout, opt, var);
    delml = true;
  }
  if (opt.ngp("line"))          // 回帰直線
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline\n";
      delml = true;
    }
    plot_ngpid(fout, opt, var);
    plot_line(fout, opt, var);
  }
  if (opt.ngp("pl"))            // 回帰直線の予測限界
  {
    if (delml) fout << NGPdelm << "\n";
    else fout << "-s1 -dline\n";
    plot_ngpid(fout, opt, var);
    plot_pl(fout, opt, var);
  }
}

void wmraclass::plot_ngpid(ofstream& fout, option& opt, const varclass& var)
// Ngraph用にデータ識別情報を出力する
{
  id.put_title(fout, formngp);
  id.put_stra(fout, formngp);
  idd tid(id);
  tid.set(id.title(), var.np() + 1); // 項目数だけ変更する
  set_itemix(tid, opt, id.ni()); // 交互作用項目名を追加する
  fout << "%重みw       :" << tid.item(0) << "\n";
  fout << "%目的変数y   :" << tid.item(1) << "\n";
  for (size_t i = 2; i < tid.ni(); ++i) fout << "%説明変数x" << setw(3) << i - 1 << ":" << tid.item(i) << "\n";
  if (opt.pin() < 1.0) fout << "%取り込み基準:p値＜Pin=" << opt.pin() << "  追い出し基準:p値≧Pout=" << opt.pout() << "\n";
  else fout << "%取り込み基準:偏F値＞Fin=" << opt.fin() << "  追い出し基準:偏F値≦Fout=" << opt.fout() << "\n";
  for (size_t i = 2; i < id.ni(); ++i)
  {
    int xno = i - 1;
    if (opt.xno() < 1)
    { fout << "%説明変数:x" << xno << "=" << var.me(xno) << "(平均値)"; }
    else
    {
      if (xno == opt.xno())
      { fout << "%X軸にする説明変数:x" << xno; }
      else
      {
        fout << "%X軸にしない説明変数:x" << xno << "=";
        if (opt.x(xno) == BLANKf) fout << var.me(xno) << "(平均値)";
        else fout << opt.x(xno);
      }
    }
    if (var.check(i)) fout << "(選択)\n";
    else fout << "(非選択)\n";
  }
}

void wmraclass::plot_est(ofstream& fout, option& opt, const varclass& var)
// Yの推定値／補正値をNgraph用に出力する
{
  if (opt.xno() > 0)            // X軸にする説明変数が指定されている時
  {
    datac tdata(var.ni());
    tdata.assign(0, var.me(0));  tdata.assign(1, var.me(var.np()));
    for (int i = 2; i < var.ni(); ++i)
    {
      int xno = i - 1;
      tdata.assign(i, (opt.x(xno) == BLANKf) ? var.me(xno) : opt.x(xno));
    }
    fout << "%X軸にする説明変数とY補正値\n";
    for (auto& cdatai : datalst)
    {
      double yest = var.cal_yest(opt, cdatai);
      if (yest != BLANKf)
      {
        double err = cdatai.fdata(1) - yest;
        size_t xidx = 1 + opt.xno();
        tdata.assign(xidx, cdatai.fdata(xidx));
        yest = var.cal_yest(opt, tdata);
        fout << cdatai.fdata(xidx) << " " << yest + err << "\n";
      }
    }
  }
  else                          // X軸にする説明変数が指定されていない時
  {
    fout << "%Y推定値とY\n";
    for (auto& cdatai : datalst)
    {
      double yest = var.cal_yest(opt, cdatai);
      if (yest != BLANKf) fout << yest << " " << cdatai.fdata(1) << "\n";
    }
  }
}

void wmraclass::plot_line(ofstream& fout, option& opt, const varclass& var)
// 回帰直線をNgraph用に出力する
{
  double xmin = FDmax, xmax = -FDmax;
  double xrange = search_xrange(xmin, xmax, opt, var);
  datac tdata(var.ni());
  tdata.assign(0, var.me(0));  tdata.assign(1, var.me(var.np()));
  for (int i = 2; i < var.ni(); ++i)
  {
    int xno = i - 1;
    tdata.assign(i, (opt.x(xno) == BLANKf) ? var.me(xno) : opt.x(xno));
  }
  double yest = xmin;
  if (opt.xno() > 0)
  {
    tdata.assign(1 + opt.xno(), xmin);
    yest = var.cal_yest(opt, tdata);
    fout << "%偏回帰直線\n";
  }
  else
  { fout << "%回帰直線\n"; }
  fout << xmin << " " << yest << "\n";
  yest = xmax;
  if (opt.xno() > 0)
  {
    tdata.assign(1 + opt.xno(), xmax);
    yest = var.cal_yest(opt, tdata);
  }
  fout << xmax << " " << yest << endl;
}

double wmraclass::search_xrange(double& xmin, double& xmax, option& opt, const varclass& var)
// x軸の最小値と最大値を検索して範囲を返す
{
  if (opt.xmin() != BLANKf && opt.xmax() != BLANKf) // X軸の範囲指定あり
  {
    xmin = opt.xmin();  xmax = opt.xmax();
    return xmax - xmin;         // 範囲を返す
  }
  xmin = FDmax;  xmax = -FDmax;
  for (auto& cdatai : datalst)
  {
    double x = (opt.xno() > 0 && var.check(opt.xno())) ? cdatai.fdata(1 + opt.xno()) : var.cal_yest(opt, cdatai);
    if (x != BLANKf)            // xの最大値と最小値をサーチ
    {
      if (x < xmin) xmin = x;
      if (x > xmax) xmax = x;
    }
  }
  double xrange = xmax - xmin, xmg = 0.5;
  xmin -= xrange * xmg;  xmax += xrange * xmg;
  return xmax - xmin;           // 範囲を返す
}

void wmraclass::plot_pl(ofstream& fout, option& opt, const varclass& var)
// 予測限界をNgraph用に出力する
{
  double xmin = FDmax, xmax = -FDmax;
  double xrange = search_xrange(xmin, xmax, opt, var);
  datac mdata(var.ni());
  mdata.assign(0, var.wyxme(0));  mdata.assign(1, var.me(var.np()));
  for (int i = 2; i < var.ni(); ++i)
  {
    int xno = i - 1;
    mdata.assign(i, (opt.x(xno) == BLANKf) ? var.me(xno) : opt.x(xno));
  }
  double vr = var.vr();
  double tcc = ptot(1.0 - opt.cc(), var.dfr());
  if (opt.tstat() == tstatc) // 検定統計量としてχ2乗値を使う
  {
    vr = 1.0;
    tcc = ptoz(1.0 - opt.cc());
  }
  statv<double> wyeststat;      // W・Y推定値の基礎統計量
  for (auto& cdatai : datalst)
  {
    double yest = var.cal_yest(opt,  cdatai);
    if (yest != BLANKf) wyeststat.sigma(cdatai.fdata(0) * yest);
  }
  fout << "%" << opt.cc() * 100.0 << "%予測限界上限\n";
  plot_pllu(fout, opt, mdata, var, tcc, vr, xmin, xmax, wyeststat, true);
  fout << NGPdelm2 << "\n";
  fout << "%" << opt.cc() * 100.0 << "%予測限界下限\n";
  plot_pllu(fout, opt, mdata, var, tcc, vr, xmin, xmax, wyeststat, false);
}

void wmraclass::plot_pllu(ofstream& fout, option& opt, datac& cdata, const varclass& var, double t, double vr, double xmin, double xmax, statv<double>& wyeststat, bool cup) const
// 下限／上限の予測限界をNgraph用に出力する
{
  constexpr int RDdef = 50;     // プロットの繰り返し数
  double dx = (xmax - xmin) / (double)RDdef;
  double yest = xmin, xi = xmin;
  for (int i = 0; i <= RDdef; ++i, xi += dx)
  {
    if (opt.xno() > 0)
    {
      cdata.assign(1 + opt.xno(), xi);
      yest = var.cal_yest(opt, cdata);
    }
    else
    { yest = xi; }
    double pi = var.cal_pi(cdata, t, vr, opt, yest, wyeststat);
    if (cup) yest += pi;
    else yest -= pi;
    fout << xi << " " << yest << "\n";
  }
}

int wmraclass::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 wmraclass::put(option& opt) // 集計結果を出力する
{
  varclass var(id.ni(), opt);   // 変数オブジェクト
  var.cal(opt, datalst);        // 各種の計算実行
  if (opt.res("res"))
  {
    put_resid(cout, opt);
    if (opt.res("stat"))
    {
      put_stat(cout, var);
      putmat(cout, "・相関行列(correlation matrix)", matr, var.n(), var.smat(), var.nxix(), 1);
    }
    if (var.resvar() == reserr)
    { cout << "※重回帰式計算不可能！\n\n" << endl; }
  }
  if (var.resvar() == reserr) return;
  var.stepsel(opt.fin(), opt.fout(), opt.pin(), opt.pout());
  if (opt.res("res") && opt.res("selx")) put_regtab(cout, opt, opt.fin(), opt.fout(), opt.pin(), opt.pout(), var);
  if (opt.dnsname() != "") put_dns(opt, var);
  if (opt.ngpname() != "") plot_regngp(opt, var);
  if (opt.res("res"))
  {
    if (opt.res("allx") && (opt.res("selx") == false || var.dfb() < var.nxix()))
    {                           // 全変数選択結果出力
      if (var.dfb() < var.nxix()) var.stepsel(0.0, 0.0, 1.0, 1.0);
      put_regtab(cout, opt, 0.0, 0.0, 1.0, 1.0, var);
    }
    cout << "\n" << endl;
  }
}

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