分散共分散行列からスタートすることもできる主成分分析です。以下のスクリプトの他に，固有値・固有ベクトルを求める関数が必要です。

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

BEGIN { # 環境設定
   FS = "[, \t]+" # 数値の区切りは，カンマ，空白，タブ
   covariance = 0 # 分散共分散行列から出発するときは 1 にする
   eps = 1e-10    # 固有値の精度（大きすぎないように）
}

{ # データファイル処理
   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++) {
       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 covariance ? "\n\n分散共分散行列\n\n" : "\n\n相関係数行列\n\n"
   for (i = 0; i < nv; i++) {
       printf "X%02i", i+1
       for (j = 0; j <= i; j++) {
           if (covariance) {
               a[i,j] /= nc-1 # 分散共分散行列
           }
           else 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"
   if (eigen2(a, nv, nv, eps, e, v) == 0) { # 固有値・固有ベクトル
       output("固有値・固有ベクトル", nv, e, 1)
       output("因子負荷量", nv, e, 0)
   }
}

function output(string, n, e, f_eval,     i, j) # 結果出力
{
   printf "\n%s\n\n   ", string
   for (i = 0; i < n; i++) {
       printf "%9i", i+1
   }
   printf "\n"
   for (j = 0; j < n; j++) {
       printf "X%02i", j+1
       for (i = 0; i < n; i++) {
           if (covariance) {
               printf "%9.5f", v[i,j]*(f_eval ? 1 : sqrt(e[i])/sd[j])
           }
           else {
               printf "%9.5f", v[i,j]*(f_eval ? 1 : sqrt(e[i]))
           }
       }
       printf "\n"
   }
   printf "EV "
   for (i = 0; i < n; i++) {
       printf "%9.5f", e[i]
   }
   printf "\n"
}