データファイルは，「複数個の独立変数，１個の従属変数」の順とします。inv 関数を呼ぶまでは，pca.awk とほとんど同じです。このスクリプトをほんの少し手直しすると，多項式回帰分析を行うことができます。

　また，偏回帰係数の検定，分散分析表の検定の有意確率を求めるためには txp 関数，fxp 関数を使用できます。

　使用法は jgawk -f mreg.awk -f inv.awk test.dat のようになっています。

# mreg.awk

BEGIN { # 環境設定
   FS = "[, \t]+" # 数値の区切りは，カンマ，空白，タブ
   epsinv = 1e-8  # 特異度の判定値
}

{ # データファイル処理
   if (NR == 1) { # データファイルの 1 行目の変数の個数を記憶
       nv = NF
   }
   else if (NF != nv) { # 2 行目以降の変数の個数のチェック
       printf "１ケースあたりの変数の個数が違います！\n\n" > "/dev/stderr"
       error = 1
       exit(1)
   }
   nc++
   for (i = 0; i < nv; i++) { # 平均値，変動の積算
       temp[i] = vx = $(i+1)-mean[i]
       mean[i] += vx/nc
       for (k = 0; k <= i; k++) {
           a[i,k] += (nc-1)*temp[k]*vx/nc
       }
   }
}

END { # 前処理が終了したら，本処理
   if (error) {
       exit(1)
   }
   nc = NR
   printf "ケース数： %i　　変数の個数： %i\n\n", nc, nv
   printf "%4s%14s%15s%15s\n", "変数", "平均値", "不偏分散", "標準偏差"
   for (i = 0; i < nv; i++) {
       ss[i] = a[nv-1,i]
       sd[i] = sqrt(var = a[i,i]/(nc-1))
       printf "X%02i%15.7g%15.7g%15.7g\n", i+1, mean[i], var, sd[i]
   }
   printf "\n\n相関係数行列\n\n"
   for (i = 0; i < nv; i++) {
       printf "X%02i", i+1
       for (j = 0; j <= i; j++) {
           if (sd[i]*sd[j] == 0) {
               exit(1)
           }
           else {
               a[i,j] /= sd[i]*sd[j]*(nc-1)
           }
           printf "%9.5f", a[j,i] = a[i,j]
       }
       printf "\n"
   }
   printf "   "
   for (i = 1; i <= nv; i++) {
       printf "      X%02i", i
   }
   printf "\n\n"
   nv1 = nv-1
   nc1 = nc-1
   if (inv(a, nv1, epsinv) == 0) { # 逆行列
       b0 = mean[nv1]
       se = syy = ss[nv1]
       for (i = 0; i < nv1; i++) {
           temp2 = 0
           for (j = 0; j < nv1; j++) {
               temp2 += a[i,j]*a[j,nv1]
           }
           b1[i] = temp2
           b[i] = temp2*sd[nv1]/sd[i]
           b0 -= b[i]*mean[i]
           se -= b[i]*ss[i]
       }
       sr = syy-se
       rhoy = se/syy
       seb0 = 0
       df2 = nc-nv1-1
       for (i = 0; i < nv1; i++) {
           seb = a[i,i]*rhoy/df2*(sd[nv1]^2/sd[i]^2)
           seb0 += (mean[i]^2)*seb
           stde[i] = sqrt(seb)
       }
       seb = 0
       for (i = 0; i < nv1; i++) {
           temp2 = 2*mean[i]/sd[i]
           for (j = i+1; j < nv1; j++) {
               seb += a[j,i]*temp2*mean[j]/sd[j]
           }
       }
       seb0 = sqrt(seb0+(seb+nc1/nc)*sd[nv1]^2*rhoy/df2)
       printf "%6s%15s%15s%15s%18s\n",
              "変数", "偏回帰係数", "標準誤差", "t値", "標準化偏回帰係数"
       for (i = 0; i < nv1; i++) {
           printf "   X%02i%15.5f%15.5f%15.5f%15.5f\n",
                  i+1, b[i], stde[i], b[i]/stde[i], b1[i]
       }
       printf "定数項%15.5f%15.5f%15.5f\n", b0, seb0, b0/seb0
       printf "\n分散分析表\n\n%4s%15s%8s%15s%15s\n",
              "要因", "平方和", "自由度", "平均平方", "Ｆ値"
       printf "回帰%15g%8i%15g%15g\n", sr, nv1, sr/nv1, sr*df2/nv1/se
       printf "残差%15g%8i%15g\n", se, df2, se/df2
       printf "全体%15g%8i%15g\n\n", syy, nc1, syy/nc1
       printf "重相関係数 = %.5f   決定係数 = %.5f", sqrt(1-rhoy), 1-rhoy
   }
}