// coreg.C Version 7.44  Copyleft (c) 杉本典夫
// Time-stamp: <2025-10-28 16:56:18 tac>
//
// 相関係数と回帰直線を計算する。
//
// 使用ライブラリ:libdans.so, libfun.so
// コンパイル法 :g++ -O -s -W coreg coreg.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 =            // プログラム使用法表示用
{
  "coreg  Version 7.44          Copyleft(Vo\\oV) 2025 杉本典夫",
  "書式:coreg  [<オプション> …]  [<データファイル> …]",
  "機能:相関係数と回帰直線を計算したるでな。",
  "<オプション>",
  "  -cc[=<信頼係数>]:信頼係数を指定するんじゃ。(0≦信頼係数<1)",
  "    -cc:-cc=0.95(オプション省略時)",
  "  -f=<ファイル名>:指定されたファイルからオプションを読み込むんじゃ。",
  "  -h[=<コマンド名>]:オンラインマニュアルを表示したる。(-h:-h=man)",
  "  -ngp={ce|cl|data|line|ma|pe|pl|te|tl}[,…]:Ngraph用データの指定。",
  "    ce:信頼楕円  cl:信頼限界  data:個々のデータ  line:回帰直線",
  "    ma:主軸回帰直線  pe:予測楕円  pl:予測限界  te:許容楕円  tl:許容限界",
  "    オプション省略時:全て指定",
  "  -ngpname[=<出力ファイル名>]:Ngraph用データを出力するファイル名を指定する。",
  "    オプション省略時:Ngraph用データ出力せず  -ngpname:-ngpname=coreg_ngp.dat",
  "  -res=[-][{def|all|ce|cl|dw[s]|icc|r|reg|sdr|stat|x[<指定値>]|y[<指定値>]}]",
  "       [,…]:解析結果の出力項目を指定するんじゃ。",
  "    -:-を付けた項目を出力しない(-res=-:全ての項目を出力しない)",
  "    def:オプション省略時と同じ項目を出力する  all:全ての項目を出力する",
  "    ce:信頼楕円  cl:信頼限界  dw[s]:ダービン・ワトソン比(dws:x順にソートする)",
  "    icc:級内相関係数  r:相関係数  reg:回帰直線  sdr:回帰残差  stat:基礎統計量",
  "    x[<指定値>]:x=平均値または<指定値>の時のy推定値",
  "    y[<指定値>]:y=平均値または<指定値>の時のx逆推定値",
  "    オプション省略時:データに合わせて自動判定",
  "  -xrange=<最小値>-<最小値>:X軸の最小値と最大値を指定します。"
};

const string Title = "=== 相関係数と回帰直線 ==="; // 表題
const string NGPNAMEdef = "coreg_ngp.dat"; // 省略時Ngraphファイル名
constexpr uint RDdef = 50;      // グラフの描画範囲分割数省略時値

enum dwtype                     // ダービン・ワトソン比の種類
{
  dwnon,                        // なし
  dwput,                        // 出力
  dwsort,                       // x順にソートした値を出力
  dwetc                         // その他
};

enum limtype                    // 限界の種類
{
  limnon,                       // なし
  limconf,                      // 信頼限界
  limtole,                      // 許容限界
  limpred,                      // 予測限界
  limetc                        // その他
};

class option                    // オプション集合体
{
  string v_path;                // プログラムのパス
  string v_pname;               // プログラム名
  double v_cc;                  // 信頼係数
  dwtype v_dw;                  // ダービン・ワトソン比出力フラッグ
  string v_ngpname;             // Ngraphファイル名
  lstr ngplst;                  // Ngraph出力リスト
  mapsb ngpmap;                 // Ngraph出力の有無マップ
  lstr reslst;                  // 分析結果出力リスト
  mapsb resmap;                 // 分析結果出力の有無マップ
  double v_x;                   // yの値を推定するxの値
  double v_y;                   // xの値を逆推定するyの値
  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; }
  dwtype dw() const          // ダービン・ワトソン比の種類を返す
  { return v_dw; }
  const string& ngpname() const // Ngraphデータファイル名を返す
  { return v_ngpname; }
  bool ngp(const string& code); // Ngraph出力フラッグを返す
  bool res(const string& code); // 分析結果出力フラッグを返す
  double x() const              // yの値を推定するxの値を返す
  { return v_x; }
  double xmin() const           // X軸の最小値を返す
  { return xrange.first; }
  double xmax() const           // X軸の最大値を返す
  { return xrange.second; }
  double y() const              // xの値を逆推定するyの値を返す
  { return v_y; }
  void setdef()                 // 省略時オプション指定を設定する
  { getoptdef(*this); }
  void set(const string& optstr, char devc = OPTdev);
  // オプション指定によってオプションを設定する
  int check(size_t ni);         // オプションのチェックとセット
  int set_ngp();                // -ngpオプションのセット
  int set_res();                // -resオプションのセット
};

option::option(const char* c_pname) : v_path(), v_pname(), v_cc(CCdef), v_dw(dwsort), v_ngpname(), ngplst(), ngpmap(), reslst(), resmap(), v_x(DBL_MAX), xrange(BLANKf, BLANKf), v_y(DBL_MAX)
// コンストラクタ
{
  if (c_pname == nullptr) return; // 初期設定せず
  v_pname = getname(c_pname);     // パス名なしのプログラム名取得
  v_path = c_pname;
  strsize_t idx = v_path.find(v_pname);
  if (idx != string::npos) v_path.erase(idx); // プログラム名削除
}

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

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

void option::set(const string& optstr, char devc)
// オプション指定によってオプションを設定する
{
  string opt[2];                // オプションのopt[0]:項目,opt[1]:内容
  if (getopt(opt, optstr, devc) < 1) return; // オプション指定なし
  int err = 0;
  if (opt[0] == "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] == "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] == "res")     // 分析結果の出力指定
  {
    if (strsplit(reslst, opt[1], CHRdelm) < 1) err = badopt(pname(), optstr, USAGE);
  }
  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                          // 間違ったオプション指定
  { 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) throw runtime_error("項目は2つなんや！");
  if (set_ngp() < 0) return -2; // -ngpオプションのセット
  if (set_res() < 0) return -3; // -resオプションのセット
  return 0;                     // エラーなし
}

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

int option::set_res()           // -resオプションのセット
{
  resmap["res"] = resmap["def"] = resmap["r"] = resmap["reg"] = resmap["stat"] = true;
  resmap["ce"] = resmap["cl"] = resmap["dw"] = resmap["icc"] = resmap["sdr"] = resmap["x"] = resmap["y"] = false;
  if (reslst.size() > 0)
  {
    resmap["def"] = resmap["r"] = resmap["reg"] = resmap["stat"] = 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["r"] = resmap["reg"] = resmap["stat"] = true; }
      else if (resstr == "all") // 全て出力
      {
        resmap["def"] = resmap["ce"] = resmap["cl"] = resmap["dw"] = resmap["icc"] = resmap["r"] = resmap["reg"] = resmap["sdr"] = resmap["stat"] = resmap["x"] = resmap["y"] = true;
      }
      else if (resstr == "dw")     // ソートしない時のダービン・ワトソン比
      { resmap[resstr] = true;  v_dw = dwput; }
      else if (resstr == "dws")    // ソートした時のダービン・ワトソン比
      { resmap["dw"] = true;  v_dw = dwsort; }
      else if (resstr[0] == 'x' && resstr.size() > 1) // x=指定値の時のy推定値を出力
      {
        v_x = strtodf(resstr.substr(1));
        if (v_x == BLANKf) return disperr(pname(), "-res=" + resstr + "オプション指定が不適当！");
        resmap["x"] = true;
      }
      else if (resstr[0] == 'y' && resstr.size() > 1) // y=指定値の時のx逆推定値を出力
      {
        v_y = strtodf(resstr.substr(1));
        if (v_y == BLANKf) return disperr(pname(), "-res=" + resstr + "オプション指定が不適当！");
        resmap["y"] = true;
      }
      else if (respos == resmap.end())
      { return disperr(pname(), "-res=" + resstr + "オプション指定が不適当！"); }
      else
      { resmap[resstr] = true; }
    }
  }
  return 0;                     // エラーなし
}

class data_less                 // ケースデータの比較用関数オブジェクト
{
private:
  int kidx;                     // ソートキー項目のインデックス
public:
  data_less(int kidx0) : kidx(kidx0) // コンストラクタ
  {}
  ~data_less()                  // デストラクタ
  {}
  bool operator()(datac& cdata1, datac& cdata2) const; // ケースデータを比較する
};

bool data_less::operator()(datac& cdata1, datac& cdata2) const
// ケースデータを比較する
{
  int comp = 0;
  double fdata1 = cdata1.fdata(kidx);
  double fdata2 = cdata2.fdata(kidx);
  if (fdata1 < fdata2) comp = -1;
  else if (fdata1 > fdata2) comp = 1;
  return (comp < 0) ? true : false; // 小さいときだけtrueを返す
}

class varclass                  // 変数クラス
{
private:
  pairdbl range;                // 最小値と最大値
  statcor<double> statc;        // 相関基礎統計量
public:
  varclass() : range(FDmax, -FDmax), statc()
  // デフォールトコンストラクタ
  {}
  ~varclass()                   // デストラクタ
  {}
  void cal(ldata& datalst);     // データを集計する
  int n() const                 // 例数を返す
  { return statc.n(); }
  double mx() const             // Xの平均値を返す
  { return statc.mx(); }
  double my() const             // Yの平均値を返す
  { return statc.my(); }
  double sxx() const            // Xの平方和散を返す
  { return statc.sxx(); }
  double syy() const            // Yの平方和散を返す
  { return statc.syy(); }
  double sxy() const            // XとYの積和を返す
  { return statc.sxy(); }
  double vxx() const            // Xの不偏分散を返す
  { return statc.vxx(); }
  double vyy() const            // Yの不偏分散を返す
  { return statc.vyy(); }
  double vxy() const            // XとYの共分散を返す
  { return statc.vxy(); }
  double rxy() const            // XとYの相関係数を返す
  { return statc.rxy(); }
  double xmin() const           // Xの最小値を返す
  { return range.first; }
  double xmax() const           // Xの最大値を返す
  { return range.second; }
};

void varclass::cal(ldata& datalst) // データを集計を有効例数を返す
{
  ldata_itr ditr = datalst.begin();
  while (ditr != datalst.end())
  {
    double x = ditr->fdata(0), y = ditr->fdata(1);
    if (statc.sigma(x, y) < 1) // 欠損データはリストから削除する
    { ditr = datalst.erase(ditr); }
    else
    {
      if (x < range.first) range.first = x;
      if (x > range.second) range.second = x;
      ++ditr;
    }
  }
}

class coregclass                // coregメインクラス
{
private:
  string fname;                 // ファイル名
  void put_resid(ostream& fout, idd& id) const;
  // 分析結果のタイトルとデータ識別情報を出力する
  void put_r(ostream& fout, option& opt, const varclass& var) const;
  // 相関係数を出力する
  void put_icc(ostream& fout, const varclass& var) const;
  // Ebelの級内相関係数を出力する
  void put_reg(ostream& fout, option& opt, const varclass& var) const;
  // 回帰直線を出力する
  void put_esty(ostream& fout, option& opt, int n, double b0, double b1, double mx, double xy, double sxx, double vr, int dfr, double t) const;
  // yの推定値を出力する
  void put_estx(ostream& fout, option& opt, int n, double b0, double b1, double mx, double xy, double sxx, double vr, int dfr, double t) const;
  // xの逆推定値を出力する
  void put_regatab(ostream& fout, const varclass& var) const;
  // 回帰直線の分散分析結果を出力する
  void put_dw(ostream& fout, option& opt, ldata& datalst, const varclass& var) const;
  // ダービン・ワトソン比を出力する
  void plot_ngp(option& opt, idd& id, ldata& datalst, const varclass& var) const;
  // Ngraph用データ出力
  void plot_ngpid(ofstream& fout, idd& id) const;
  // Ngraph用データのタイトルとデータ識別情報を出力する
  void plot_ngpdata(ofstream& fout, ldata& datalst) const;
  // Ngraph用に個々のデータを出力する
  void plot_ngpline(ofstream& fout, option& opt, const varclass& var) const;
  // Ngraph用に回帰直線を出力する
  double set_xrange(double& xmin, double& xmax, option& opt, const varclass& var) const;
  // X軸の最小値と最大値を設定して範囲を返す
  void plot_ngplimit(ofstream& fout, option& opt, const varclass& var, limtype limt) const;
  // Ngraph用に信頼・許容・予測限界のデータを出力する
  void plot_ngpelp(ofstream& fout, option& opt, const varclass& var, limtype limt) const;
  // Ngraph用に信頼・許容・予測楕円のデータを出力する
  void plot_ngpma(ofstream& fout, option& opt, const varclass& var) const;
  // Ngraph用に主軸回帰直線を出力する
public:
  coregclass() : fname()        // コンストラクタ
  {}
  ~coregclass()                 // デストラクタ
  {}
  int get(const char* c_fname, option& opt, idd& id, ldata& datalst);
  // データを読み込み累積集計する
  void put(option& opt, idd& id, ldata& datalst) const;
  // 集計結果を出力する
};

void coregclass::put_resid(ostream& fout, idd& id) const
// 分析結果のタイトルとデータ識別情報を出力する
{
  puttitle(fout, Title);
  fout << "\n";
  id.put_title(fout);
  id.put_stra(fout);
  fout << "説明変数x:" << id.item(0) << "\n";
  fout << "目的変数y:" << id.item(1) << "\n";
}

void coregclass::put_r(ostream& fout, option& opt, const varclass& var) const
// 相関係数を出力する
{
  int n = var.n();
  int dfr = n - 2;
  double r = var.rxy();
  fout << "\n・相関係数(correlation coefficient):r=" << r << "\n";
  fout << "  寄与率(coefficeint of determination):r2乗=" << r * r << "\n";
  if (n > 3)
  {
    double zr = 0.0;
    double ci = calrci(zr, opt.cc(), n, r);
    fout << "  相関係数の" << opt.cc() * 100.0 << "%信頼区間(Fisherのz変換):下限=" << tanh(zr - ci) << " 上限=" << tanh(zr + ci) << "\n";
    zr *= sqrt((double)n - 3.0);
    ldbl pr = ztop((ldbl)zr);
    fout << "  相関係数の検定(Fisherのz変換):z=" << zr << "  有意確率:p=" << pr << pind(pr) << "\n";
  }
  if (n < 3 || opt.res("ce") == false) return;
  double vxx = var.vxx(), vyy = var.vyy(), vxy = var.vxy();
  double cf = 1.0 / (double)n, rf = 1.0 + 1.0 / (double)n;
  double ccpr = opt.cc() * 100.0;
  double a = 0.0, b = 0.0, h = 0.0, deg = 180.0 / M_PI;
  calelp(a, b, h, opt.cc(), 2, dfr, vxx * cf, vyy * cf, vxy * cf);
  double ba = tan(h);
  fout << "  " << "信頼楕円の長軸=主軸回帰(Major axis regression):y=" << var.my() - ba * var.mx() << showpos << ba << noshowpos << "x\n";
  fout << "  " << ccpr << "%信頼楕円:長軸=" << a << "  短軸=" << b << "  長軸の回転角(度)=" << h * deg << "°\n";
  calelp(a, b, h, opt.cc(), 2, dfr, vxx, vyy, vxy);
  fout << "  " << ccpr << "%許容楕円:長軸=" << a << "  短軸=" << b << "  長軸の回転角(度)=" << h * deg << "°\n";
  calelp(a, b, h, opt.cc(), 2, dfr, vxx * rf, vyy * rf, vxy * rf);
  fout << "  " << ccpr << "%予測楕円:長軸=" << a << "  短軸=" << b << "  長軸の回転角(度)=" << h * deg << "°\n";
}

void coregclass::put_icc(ostream& fout, const varclass& var) const
// Ebelの級内相関係数を出力する
{
  int n = var.n();
  int dfr = n - 2;
  double r = var.rxy();
  double sxx = var.sxx(), syy = var.syy(), sxy = var.sxy();
  double r11 = (sxx + syy > 0.0) ? 2.0 * sxy / (sxx + syy) : 0.0;
  double r22 = (sxx + syy + 2.0 * sxy > 0.0) ? 4.0 * sxy / (sxx + syy + 2.0 * sxy) : 0.0;
  double rv = (sxx + syy > 0.0) ? 2.0 * sqrt(sxx * syy) / (sxx + syy) : 0.0;
  fout << "\n・Ebelの級内相関係数(Ebel's intraclass correlation coefficient)\n";
  fout << "  r11=" << r11 << "  r22=" << r22 << "  分散一致係数:rv=" << rv << endl;
}

void coregclass::put_reg(ostream& fout, option& opt, const varclass& var) const
// 回帰直線を出力する
{
  int n = var.n();
  int dfr = n - 2;
  double mx = var.mx(), my = var.my();
  double sxx = var.sxx(), syy = var.syy(), sxy = var.sxy();
  double b1 = (sxx > 0.0) ? sxy / sxx : 0.0;
  double b0 = my - b1 * mx;
  double sr = syy - sxy * b1;
  if (sr < FDeps) sr = 0.0;
  double vr = (dfr > 0) ? sr / (double)dfr : 0.0;
  double vb = (sxx > 0.0) ? vr / sxx : 0.0;
  double se0 = (sxx > 0.0) ? sqrt(vr * (1.0 / n + mx * mx / sxx)) : 0.0;
  double seb = (sxx > 0.0) ? sqrt(vr / sxx) : 0.0;
  fout << "\n・回帰直線(regression line):y=" << b0 << showpos << b1 << noshowpos << "x\n";
  fout << "  切片の標準誤差=" << se0 << "  回帰係数の標準誤差=" << seb << "\n";
  if (opt.res("sdr")) fout << "  回帰残差=" << vr << "  回帰残差の標準偏差=" << sqrt(vr) << "\n";
  if (n < 3) return;
  double cf = 1.0 / (double)n, rf = 1.0 + 1.0 / (double)n;
  double t = (dfr > 0) ? (double)(ptot((ldbl)(1.0 - opt.cc()), dfr)) : 0.0;
  double ccpr = opt.cc() * 100.0;
  double ci = calmcci(mcfis, opt.cc(), 1, 1, dfr, 1, vb);
  fout << "  回帰係数の" << ccpr << "%信頼区間=";
  fout << b1 << "±" << ci << "(" << b1 - ci << "-" << b1 + ci << ")\n";
  if (opt.res("cl"))
  {
    fout << "  回帰直線の" << ccpr << "%信頼限界(平均値の信頼区間):\n";
    fout << "    y=" << b0 << showpos << b1 << noshowpos << "x±";
    fout << t << "√{" << vr * cf << showpos << vb << noshowpos << "(x" << showpos << -mx << noshowpos << ")^2}\n";
  }
  if (opt.res("x")) put_esty(fout, opt, n, b0, b1, mx, my, sxx, vr, dfr, t);
  if (opt.res("y")) put_estx(fout, opt, n, b0, b1, mx, my, sxx, vr, dfr, t);
  if (opt.res("cl"))
  {
    fout << "  回帰直線の" << ccpr << "%許容限界(データの信頼区間):\n";
    fout << "    y=" << b0 << showpos << b1 << noshowpos << "x±";
    fout << t << "√(" << vr << ")\n";
    fout << "  回帰直線の" << ccpr << "%予測限界(平均値+データの信頼区間):\n";
    fout << "    y=" << b0 << showpos << b1 << noshowpos << "x±";
    fout << t << "√{" << vr * rf << showpos << vb << noshowpos << "(x" << showpos << -mx << noshowpos << ")^2}\n";
  }
}

void coregclass::put_esty(ostream& fout, option& opt, int n, double b0, double b1, double mx, double my, double sxx, double vr, int dfr, double t) const
// yの推定値を出力する
{
  double x = (opt.x() == DBL_MAX) ? mx : opt.x();
  double y = b0 + b1 * x, dx = x - mx;
  double sey = sqrt(vr * (1.0 / (double)n + dx * dx / sxx));
  double ci = t * sey;
  fout << "    x=" << x << "の時:y=" << y << "±" << ci << "(" << y - ci << "-" << y + ci << ")\n";
}

void coregclass::put_estx(ostream& fout, option& opt, int n, double b0, double b1, double mx, double my, double sxx, double vr, int dfr, double t) const
// xの逆推定値を出力する
{
  if (sgn(b1) == 0) return;     // 逆推定できず
  double y = (opt.y() == DBL_MAX) ? my : opt.y();
  double cf = (y - my) / b1;
  double x = mx + cf;
  double g = t / b1 * t / b1 * vr / sxx;
  if (g >= 1.0) g = 0.0;        // g>1だと計算できないため
  double gr = 1.0 - g;
  double v = vr * (gr / (double)n + cf * cf / sxx);
  double ci = t / b1 * sqrt(v);
  double xl = (x - g * mx - ci) / gr, xu = (x - g * mx + ci) / gr;
  fout << "    y=" << y << "の時:x=" << x << "  下限=" << xl << "  上限=" << xu << "\n";
}

void coregclass::put_regatab(ostream& fout, const varclass& var) const
// 回帰直線の分散分析結果を出力する
{
  int n = var.n();
  int dfr = n - 2;
  double mx = var.mx(), my = var.my();
  double sxx = var.sxx(), syy = var.syy(), sxy = var.sxy();
  double b1 = (sxx > 0.0) ? sxy / sxx : 0.0;
  double b0 = my - b1 * mx;
  double sr = syy - sxy * b1;
  if (sr < FDeps) sr = 0.0;
  fout << "\n・回帰直線の分散分析(=回帰係数の検定=相関係数の検定)\n";
  putatab(fout, "回帰", sxy * b1, 1, sr, dfr, syy, n - 1);
}

void coregclass::put_dw(ostream& fout, option& opt, ldata& datalst, const varclass& var) const
// ダービン・ワトソン比を出力する
{
  int n = var.n();
  double mx = var.mx(), my = var.my();
  double sxx = var.sxx(), syy = var.syy(), sxy = var.sxy();
  double b1 = (sxx > 0.0) ? sxy / sxx : 0.0;
  double b0 = my - b1 * mx;
  double sr = syy - sxy * b1;
//  if (sr < FDeps) return;                             // 計算できず
  if (opt.dw() == dwsort) datalst.sort(data_less(0)); // x順にソートする
  ldata_itr ditr = datalst.begin();
  double e1 = ditr->fdata(1) - (b0 + b1 * ditr->fdata(0));
  double sume2 = e1 * e1, sumed = 0.0, epre = e1;
  ++ditr;
  while (ditr != datalst.end())
  {
    double ei = ditr->fdata(1) - (b0 + b1 * ditr->fdata(0));
    double ed = ei - epre;
    sume2 += ei * ei;  sumed += ed * ed;
    epre = ei;  ++ditr;
  }
  double dw = (sgn(sume2) > 0) ? sumed / sume2 : 0.0;
  fout << "\n・ダービン・ワトソン比";
  if (opt.dw() == dwsort) fout << "(x順にソート後)";
  fout << ":DW=" << dw << "\n";
}

void coregclass::plot_ngp(option& opt, idd& id, ldata& datalst, const varclass& var) const
// Ngraph用データ出力
{
  ofstream fout;
  if (ofopen(fout, (opt.ngpname()).c_str()) == false) return; // ファイルオープン失敗
  bool delml = false;
  if (opt.ngp("data"))          // データ
  {
    fout << "-s1 -dmark,00\n";
    plot_ngpid(fout, id);
    plot_ngpdata(fout, datalst);
    delml = true;
  }
  if (opt.ngp("line"))          // 回帰直線
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline\n";
      delml = true;
    }
    plot_ngpid(fout, id);
    plot_ngpline(fout, opt, var);
  }
  if (opt.ngp("cl"))            // 信頼限界
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline\n";
      delml = true;
    }
    plot_ngpid(fout, id);
    plot_ngplimit(fout, opt, var, limconf);
  }
  if (opt.ngp("tl"))            // 許容限界
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline\n";
      delml = true;
    }
    plot_ngpid(fout, id);
    plot_ngplimit(fout, opt, var, limtole);
  }
  if (opt.ngp("pl"))            // 予測限界
  {
    if (delml)
    { fout << NGPdelm << "\n"; }
    else
    {
      fout << "-s1 -dline\n";
      delml = true;
    }
    plot_ngpid(fout, id);
    plot_ngplimit(fout, opt, var, limpred);
  }
  if (opt.ngp("ce"))            // 信頼楕円
  {
    if (delml) fout << NGPdelm << "\n";
    else fout << "-s1 -dpolygon\n";
    plot_ngpid(fout, id);
    plot_ngpelp(fout, opt, var, limconf);
  }
  if (opt.ngp("te"))            // 許容楕円
  {
    if (delml) fout << NGPdelm << "\n";
    else fout << "-s1 -dpolygon\n";
    plot_ngpid(fout, id);
    plot_ngpelp(fout, opt, var, limtole);
  }
  if (opt.ngp("pe"))            // 予測楕円
  {
    if (delml) fout << NGPdelm << "\n";
    else fout << "-s1 -dpolygon\n";
    plot_ngpid(fout, id);
    plot_ngpelp(fout, opt, var, limpred);
  }
  if (opt.ngp("ma"))            // 主軸回帰直線
  {
    if (delml) fout << NGPdelm << "\n";
    else fout << "-s1 -dline\n";
    plot_ngpid(fout, id);
    plot_ngpma(fout, opt, var);
  }
}

void coregclass::plot_ngpid(ofstream& fout, idd& id) const
// Ngraph用データのタイトルとデータ識別情報を出力する
{
  fout << "%";  puttitle(fout, Title);
  fout << "%\n";
  id.put_title(fout, formngp);
  id.put_stra(fout, formngp);
  fout << "%説明変数x:" << id.item(0) << "\n";
  fout << "%目的変数y:" << id.item(1) << "\n";
}

void coregclass::plot_ngpdata(ofstream& fout, ldata& datalst) const
// Ngraph用に個々のデータを出力する
{
  fout << "%生データ\n";
  for (auto& cdatai : datalst)
  { fout << cdatai.fdata(0) << " " << cdatai.fdata(1) << "\n"; }
}

void coregclass::plot_ngpline(ofstream& fout, option& opt, const varclass& var) const
// Ngraph用に回帰直線を出力する
{
  fout << "%回帰直線\n";
  int n = var.n();
  double mx = var.mx(), my = var.my();
  double sxx = var.sxx(), syy = var.syy(), sxy = var.sxy();
  double b = (sxx > 0.0) ? sxy / sxx : 0.0;
  double a = my - b * mx;
  double xmin = 0.0, xmax = 0.0;
  double xrange = set_xrange(xmin, xmax, opt, var);
  fout << xmin << " " << a + b * xmin << "\n";
  fout << xmax << " " << a + b * xmax << endl;
}

double coregclass::set_xrange(double& xmin, double& xmax, option& opt, const varclass& var) const
// X軸の最小値と最大値を設定して範囲を返す
{
  if (opt.xmin() != BLANKf && opt.xmax() != BLANKf) // X軸の範囲指定あり
  { xmin = opt.xmin();  xmax = opt.xmax(); }
  else
  {
    xmin = var.xmin();  xmax = var.xmax();
    double xrange = xmax - xmin, xmg = 0.5;
    xmin -= xrange * xmg;  xmax += xrange * xmg;
  }
  return xmax - xmin;           // 範囲を返す
}

void coregclass::plot_ngplimit(ofstream& fout, option& opt, const varclass& var, limtype limt) const
// Ngraph用に信頼・許容・予測限界のデータを出力する
{
  int dfr = var.n() - 2;
  double mx = var.mx(), my = var.my();
  double sxx = var.sxx(), syy = var.syy(), sxy = var.sxy();
  double b1 = (sxx > 0.0) ? sxy / sxx : 0.0;
  double b0 = my - b1 * mx;
  double sr = syy - sxy * b1;
  if (sr < FDeps) sr = 0.0;
  double vr = (dfr > 0) ? sr / (double)dfr : 0.0;
  double vb = (sxx > 0.0) ? vr / sxx : 0.0, vc = vr;
  string title("許容限界(データの信頼区間)");
  switch (limt)
  {
  case limconf : title = "信頼限界(平均値の信頼区間)";
    vc = vr * 1.0 / (double)(var.n());
    break;
  case limpred : title = "予測限界(平均値+データの信頼区間)";
    vc = vr * (1.0 + 1.0 / (double)(var.n()));
  }
  double t = ptot(1.0 - opt.cc(), dfr);
  double ccpr = opt.cc() * 100.0;
  fout << "%回帰直線の" << ccpr << "%" << title << "上限\n";
  double xmin = 0.0, xmax = 0.0;
  double xrange = set_xrange(xmin, xmax, opt, var);
  double dx = xrange / (double)RDdef;
  double x = xmin, dm = 0.0, y = 0.0;
  if (limt == limtole)          // 許容限界
  {
    y = b0 + b1 * xmin + t * sqrt(vc);
    fout << xmin << " " << y << "\n";
    y = b0 + b1 * xmax + t * sqrt(vc);
    fout << xmax << " " << y << "\n";
  }
  else
  {
    for (uint i = 0; i <= RDdef; ++i, x += dx)
    {
      dm = x - mx;
      y = b0 + b1 * x + t * sqrt(vc + vb * dm * dm);
      fout << x << " " << y << "\n";
    }
  }
  fout << NGPdelm2 << "\n";
  fout << "%回帰直線の" << ccpr << "%" << title << "下限\n";
  if (limt == limtole)          // 許容限界
  {
    y = b0 + b1 * xmin - t * sqrt(vc);
    fout << xmin << " " << y << "\n";
    y = b0 + b1 * xmax - t * sqrt(vc);
    fout << xmax << " " << y << "\n";
  }
  else
  {
    x = xmin;
    for (uint i = 0; i <= RDdef; ++i, x += dx)
    {
      dm = x - mx;
      y = b0 + b1 * x - t * sqrt(vc + vb * dm * dm);
      fout << x << " " << y << "\n";
    }
  }
}

void coregclass::plot_ngpelp(ofstream& fout, option& opt, const varclass& var, limtype limt) const
// Ngraph用に信頼・許容・予測楕円のデータを出力する
{
  int dfr = var.n() - 2;
  double mx = var.mx(), my = var.my();
  double vxx = var.vxx(), vyy = var.vyy(), vxy = var.vxy();
  double cf = 1.0;
  string title("許容楕円");
  switch (limt)
  {
  case limconf : title = "信頼楕円";
    cf = 1.0 / (double)(var.n());
    break;
  case limpred : title = "予測楕円";
    cf = 1.0 + 1.0 / (double)(var.n());
  }
  fout << "%" << opt.cc() * 100.0 << "%" << title << "\n";
  double a = 0.0, b = 0.0, h = 0.0;
  calelp(a, b, h, opt.cc(), 2, dfr, vxx * cf, vyy * cf, vxy * cf);
  constexpr uint RD = 60;       // 楕円の分解数
  double dt = 2.0 * M_PI / (double)RD;
  double x0 = 0.0, y0 = 0.0, x = 0.0, y = 0.0, t = 0.0;
  for (uint i = 0; i<= RD; ++i, t += dt)
  {
    x0 = a * cos(t);  y0 = b * sin(t);
    x = mx + x0 * cos(h) - y0 * sin(h);  y = my + x0 * sin(h) + y0 * cos(h);
    fout << x << " " << y << "\n";
  }
}

void coregclass::plot_ngpma(ofstream& fout, option& opt, const varclass& var) const
// Ngraph用に主軸回帰直線を出力する
{
  int dfr = var.n() - 2;
  double mx = var.mx(), my = var.my();
  double vxx = var.vxx(), vyy = var.vyy(), vxy = var.vxy();
  double cf = 1.0 / (double)(var.n());
  double a = 0.0, b = 0.0, h = 0.0;
  calelp(a, b, h, opt.cc(), 2, dfr, vxx * cf, vyy * cf, vxy * cf);
  fout << "%信頼楕円の長軸=主軸回帰直線\n";
  double ba = tan(h);
  double aa = var.my() - ba * var.mx();
  double xmin = 0.0, xmax = 0.0;
  double xrange = set_xrange(xmin, xmax, opt, var);
  double y = aa + ba * xmin;
  fout << xmin << " " << y << "\n";
  y = aa + ba * xmax;
  fout << xmax << " " << y << endl;
}

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

void coregclass::put(option& opt, idd& id, ldata& datalst) const
// 集計結果を出力する
{
  varclass var;                 // 計算用変数クラス
  var.cal(datalst);
  if (opt.res("res"))           // 計算結果出力
  {
    put_resid(cout, id);
    if (var.n() < 1)
    {
      cout << "※有効なデータなし！\n\n" << endl;
      return;
    }
    if (opt.res("stat"))
    {
      cout << "\n・変数の基礎統計量\n";
      constexpr int LENline = 69; // 区切り線の長さ
      putline(cout, LENline);
      cout << "x:";  putmsd(cout, var.n(), var.mx(), var.vxx());
      cout << "y:";  putmsd(cout, var.n(), var.my(), var.vyy());
      putline(cout, LENline);
    }
    if (var.n() < 2)
    {
      cout << "※相関係数と回帰直線の計算不可能！\n\n" << endl;
      return;
    }
    if (opt.res("r")) put_r(cout, opt, var);
    if (opt.res("icc")) put_icc(cout, var);
    if (opt.res("reg")) put_reg(cout, opt, var);
    if (opt.res("dw")) put_dw(cout, opt, datalst,var);
    if (opt.res("r") || opt.res("reg")) put_regatab(cout, var);
    cout << "\n" << endl;
  }
  if (opt.ngpname() != "") plot_ngp(opt, id, datalst, var); // Ngraph用データ出力
}

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