以下のスクリプトでは、"正確に" 指定された相関係数を持つ人工データを生成します。
なお、位置および散布度の調整には対応していません。
# 与えられた相関係数行列（因子負荷行列）を再現するデータの生成
# 入力データファイルに相関係数行列（対角成分を含む下側三角行列）または因子負荷行列を用意する（これを仮に input というファイル名とする）。
# 例えば，以下のような5行からなるファイルを用意する。先頭の # は除く。
# 1.00000
# 0.40000 1.00000
# 0.30000 0.40000 1.00000
# 0.20000 0.30000 0.40000 1.00000
# 0.10000 0.20000 0.30000 0.40000 1.00000
# 下の nc = と書いた所の数字を希望する値に変更する。
# awk -f chol.awk input > output
# 結果として，output というファイルにデータが得られる。
# 参考：
# Newsgroups: comp.soft-sys.spss
# Message-ID: <304F10DD@msmailgw.spss.com>
# Subject: Re: Generating data to match a correlation matrix.
# ref. a technical note in Keywords( Aug. 1992, Number 49).
BEGIN {
nc = 100 # 標本の大きさ 希望する値に変更すること
nv = 0 # 変更不可 おまじない(^_^)
}
{
nfact = NF
for (i = 0; i < NF; i++) {
b[nv, i] = $(i+1)
}
nv++
}
END {
# 以下は，用意された入力データが，相関係数行列か因子負荷量行列かを判定する。
if (b[0, 0] == 1) { # 対角成分が1なら相関係数行列とみなし，対称行列を作る。
for (i = 0; i < nv; i++) {
for (j = 0; j <= i; j++) {
a[i, j] = a[j, i] = b[i, j]
}
}
}
else { # 因子負荷量行列のときは相関係数行列に変換する。
MakeCorrMatrix(b, a, nv, nfact)
}
# 相関係数行列のコレスキー分解を行う。
Cholesky(a, nv)
# 仮のデータ行列を作る。この時点では変数間の相関は近似的に0である。
GenerateData(z, nv, nc)
# データ行列の各変数の平均値，標準偏差を求める。
CalculateStatistics(z, nv, nc, mean, sd, r)
CopyCorrMatrix(r, or, nv)
inv(or, nv, 1e-8)
# 主成分分析を行い，主成分得点を求める。この時点で変数間の相関は完全に0である。
eigen2(r, nv, nv, 1e-10, eval, evec)
CalculateLoadings(eval, evec, fl, nv)
CalculateFactorScoreCoefficients(or, fl, coeff, nv)
Orthogonal(z, mean, sd, coeff, nv, nc)
# コレスキー分解の結果をもとに，指定された相関係数行列を持つように主成分得点を変換する。
TransformData(z, nv, nc, a)
# 生成されたデータ行列を書き出す。
for (i = 0; i < nc; i++) {
for (j = 0; j < nv; j++) {
printf "%g ", z[i, j]
}
printf "\n"
}
}
# ボックス・ミュラー法により正規乱数を発生させ，仮のデータ行列を作る。
function GenerateData(z, nv, nc, i, j, u1, u2)
{
for (i = 0; i < nc; i++) {
for (j = 0; j < nv; j++) {
u1 = sqrt(-2.0*log(rand()))
u2 = 6.283185307179586*rand() # 6.28318530717959=2*pi
z[i, j] = u1*cos(u2)
}
}
}
# データ行列の平均値，標準偏差，相関係数行列を求める。
function CalculateStatistics(z, nv, nc, mean, sd, r, i, j, k)
{
for (i = 0; i < nv; i++) {
mean[i] = sd[i] = 0
for (j = 0; j < nv; j++) {
r[i, j] = 0
}
}
for (i = 0; i < nc; i++) {
for (j = 0; j < nv; j++) {
mean[j] += z[i, j]
}
}
for (j = 0; j < nv; j++) {
mean[j] /= nc
}
for (i = 0; i < nc; i++) {
for (j = 0; j < nv; j++) {
for (k = 0; k <= j; k++) {
r[j, k] += (z[i, j]-mean[j])*(z[i, k]-mean[k])
}
}
}
for (j = 0; j < nv; j++) {
for (k = 0; k <= j; k++) {
r[j, k] /= nc
}
sd[j] = sqrt(r[j, j])
}
for (j = 0; j < nv; j++) {
for (k = 0; k <= j; k++) {
r[j, k] /= sd[j]*sd[k]
r[k, j] = r[j, k]
}
}
}
# 因子負荷行列を計算する。
function CalculateLoadings(eval, evec, fl, nv, i, j)
{
for (i = 0; i < nv; i++) {
for (j = 0; j < nv; j++) {
fl[i, j] = sqrt(eval[i])*evec[i, j]
}
}
}
# 主成分得点係数を計算する。
function CalculateFactorScoreCoefficients(or, fl, coeff, nv, i, j, k, x)
{
for (i = 0; i < nv; i++) {
for (j = 0; j < nv; j++) {
x = 0
for (k = 0; k < nv; k++) {
x += fl[i, k]*or[j, k]
}
coeff[i, j] = x
}
}
}
# 主成分得点を計算する。
function Orthogonal(z, mean, sd, coeff, nv, nc, i, j, k, x)
{
for (i = 0; i < nc; i++) {
for (j = 0; j < nv; j++) {
z[i, j] = (z[i, j]-mean[j])/sd[j]
}
for (j = 0; j < nv; j++) {
x[j] = 0
for (k = 0; k < nv; k++) {
x[j] += z[i, k]*coeff[j, k]
}
}
for (j = 0; j < nv; j++) {
z[i, j] = x[j]
}
}
}
# コレスキー分解に基づき，データ行列を変換する。
function TransformData(z, nv, nc, a, i, j, k, x)
{
for (i = 0; i < nc; i++) {
for (j = 0; j < nv; j++) {
x[j] = 0
for (k = 0; k <= j; k++) {
x[j] += z[i, k]*a[k, j]
}
}
for (j = 0; j < nv; j++) {
z[i, j] = x[j]
}
}
}
# 因子負荷量行列から相関係数行列を計算する。
function MakeCorrMatrix(b, a, nv, nfact, i, j, k, x)
{
for (i = 0; i < nv; i++) {
x = 0
for (j = 0; j < nfact; j++) {
x += b[i, j]^2
}
x = sqrt(x)
for (j = 0; j < nfact; j++) {
b[i, j] /= x
}
}
for (i = 0; i < nv; i++) {
for (j = 0; j <= i; j++) {
x = 0
for (k = 0; k < nfact; k++) {
x += b[i, k]*b[j, k]
}
a[j, i] = a[i, j] = x
}
}
}
# 行列の複写。
function CopyCorrMatrix(r, or, nv, i, j)
{
for (i = 0; i < nv; i++) {
for (j = 0; j < nv; j++) {
or[i, j] = r[i, j]
}
}
}
function Cholesky(aa, n, a, u, i, icol, ii, irow, j, k, kk, l, m, eta, eta2, w, x, nn)
{
nn = 0
for (i = 0; i < n; i++) {
for (j = 0; j <= i; j++) {
a[++nn] = aa[i,j]
}
}
eta = 1e-9
eta2 = eta*eta
j = 1
k = ii = 0
for (icol = 1; icol <= n; icol++) {
ii += icol
x = eta2*a[ii]
l = kk = 0
for (irow = 1; irow <= icol; irow++) {
kk += irow
w = a[++k]
m = j
for (i = 1; i <= irow; i++) {
l++
if (i == irow) break
w -= u[l]*u[m++]
}
if (irow == icol) break
if (u[l] != 0) {
u[k] = w/u[l]
}
else if (w*w > abs(x*a[kk])) {
exit
}
else {
u[k] = 0
}
}
if (abs(w) <= abs(eta*a[k])) {
u[k] = 0
nullty++
}
else if (w < 0) {
exit
}
else {
u[k] = sqrt(w)
}
j += icol
}
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
aa[i,j] = 0
}
}
ii = 0
for (i = 0; i < n; i++) {
for (j = 0; j <= i; j++) {
aa[j,i] = u[++ii]
}
}
}
# x の絶対値を求める
# 使用例： x = abs(-10.95)
function abs(x)
{
return (x < 0) ? -x : x
}
# x, y の小さい方を求める
# 使用例： x = min(10, 5)
function min(x, y)
{
return (x < y) ? x : y
}
# x, y の大きい方を求める
# 使用例： x = max(10, 5)
function max(x, y)
{
return (x > y) ? x : y
}
# 逆行列を求める関数 inv.awk
function inv(r, ip, epsinv, i, j, k, b, t, l, m)
{
for (k = 0; k < ip; k++) {
m[k] = l[k] = k; b = r[k, k]
for (j = k; j < ip; j++) {
for (i = k; i < ip; i++) {
if (abs(r[j, i]) > abs(b)) {
b = r[j, i]
l[k] = i
m[k] = j
}
}
}
if (abs(b) < epsinv) {
return 1
}
j = l[k]
if (j > k) {
for (i = 0; i < ip; i++) {
t = -r[i, k]
r[i, k] = r[i, j]
r[i, j] = t
}
}
i = m[k]
if (i > k) {
for (j = 0; j < ip; j++) {
t = -r[k, j]
r[k, j] = r[i, j]
r[i, j] = t
}
}
for (i = 0; i < ip; i++) {
if (i != k) {
r[k, i] /= -b
}
}
for (i = 0; i < ip; i++) {
for (j = 0; j < ip; j++) {
if (i != k && j != k) {
r[j, i] += r[k, i]*r[j, k]
}
}
}
for (j = 0; j < ip; j++) {
if (j != k) {
r[j, k] /= b
}
}
r[k, k] = 1.0/b
}
for (k = ip-2; k >= 0; k--) {
i = l[k]
if (i > k) {
for (j = 0; j < ip; j++) {
t = r[k, j]
r[k, j] = -r[i, j]
r[i, j] = t
}
}
j = m[k]
if (j > k) {
for (i = 0; i < ip; i++) {
t = r[i, k]
r[i, k] = -r[i, j]
r[i, j] = t
}
}
}
return 0
}
# 実対称行列の，固有値(e)，固有ベクトル(v)を求める関数 eigen.awk
function eigen2(a, n, ne, eps, e, v, sw, i, j, k, k1, m, nea, nm1,
nm2, nv, nva, a1, ajk, d, ei, eps1, eps2, f, q, r, s, sr, t, wj1, w, lw)
{
nv = nva = nea = abs(ne)
if (nea == 0 || nea > n) return 1
nm1 = n-1; nm2 = n-2
if (eps <= 0) {
eps = 1e-16
}
if (n == 1) {
e[0] = a[0,0]
if (nv != 0) {
v[0,0] = 1
}
return 0
}
if (n == 2) {
w[0,0] = a[0,1]; t = 0.5*(a[0,0]+a[1,1])
r = a[0,0]*a[1,1]-a[0,1]^2
d = t*t-r
q = abs(t)+sqrt(d)
if (t < 0) {
q = -q
}
t *= ne
if (t >= 0) {
e[0] = q
if (nea == 2) {
e[1] = r/q
}
}
else {
e[0] = r/q
if (nea == 2) {
e[1] = q
}
}
}
else {
for (k = 1; k <= nm2; k++) {
k1 = k-1
s = 0
for (i = k; i < n; i++) {
s += a[k1,i]^2
}
w[0,k1] = 0
if (s != 0) {
sr = sqrt(s)
a1 = a[k1,k]
if (a1 < 0) {
sr = -sr
}
w[0,k1] = -sr
r = 1/(s+a1*sr)
a[k1,k] = a1+sr
for (i = k; i < n; i++) {
s = 0
for (j = k; j <= i; j++) {
s += a[j,i]*a[k1,j]
}
if (i != nm1) {
for (j = i+1; j < n; j++) {
s += a[i,j]*a[k1,j]
}
}
w[0,i] = s*r
}
s = 0
for (i = k; i < n; i++) {
s += a[k1,i]*w[0,i]
}
t = 0.5*s*r
for (i = k; i < n; i++) {
w[0,i] -= t*a[k1,i]
}
for (j = k; j < n; j++) {
wj1 = w[0,j]
ajk = a[k1,j]
for (i = j; i < n; i++) {
a[j,i] -= a[k1,i]*wj1+w[0,i]*ajk
}
}
}
}
w[0,nm2] = a[nm2,nm1]
for (i = 0; i < n; i++) {
w[5,i] = a[i,i]
}
r = abs(w[5,0])+abs(w[0,0])
f = abs(w[0,nm2])+abs(w[5,nm1])
r = max(r, f)
for (i = 1; i < nm1; i++) {
f = abs(w[0,i-1])+abs(w[5,i])+abs(w[0,i])
r = max(r, f)
}
eps1 = r*1e-16
eps2 = r*eps
for (i = 0; i < nm1; i++) {
w[1,i] = w[0,i]^2
}
if (ne < 0) {
r = -r
}
f = r
for (i = 0; i < nea; i++) {
e[i] = -r
}
pass2 = 0
for (k = 1; k <= nea; k++) {
if (!pass2) {
k1 = k-1
d = e[k1]
}
pass2 = 0
t = 0.5*(d+f)
if (abs(d-f) > eps2 && t != d && t != f) {
j = 0
i = 1
for (;;) {
if (!pass1) {
q = w[5,i-1]-t
}
pass1 = 0
if (q >= 0) {
j++
}
if (q != 0) {
if (++i > n) {
break
}
q = w[5,i-1]-t-w[1,i-2]/q
pass1 = 1
continue
}
else {
i += 2
}
if (i > n) {
break
}
}
if (ne < 0) {
j = n-j
}
if (j < k) {
f = t
pass2 = 1
k--
continue
}
d = t
m = min(j, nea)
for (i = k-1; i < m; i++) {
e[i] = t
}
pass2 = 1
k--
continue
}
e[k1] = t
}
}
if (nv != 0) {
if (n == 2) {
w[5,0] = a[0,0]
w[5,1] = a[1,1]
}
w[0,nm1] = 0
for (i = 0; i < nva; i++) {
ei = e[i]
for (j = 0; j < n; j++) {
w[1,j] = w[5,j]-ei
w[2,j] = w[0,j]
v[i,j] = 1
}
sw = 0
for (j = 0; j < nm1; j++) {
if (abs(w[1,j]) >= abs(w[0,j])) {
if (w[1,j] == 0) {
w[1,j] = 1e-30
}
w[4,j] = w[0,j]/w[1,j]
lw[j] = 0
w[1,j+1] -= w[4,j]*w[2,j]
w[3,j] = 0
}
else {
w[4,j] = w[1,j]/w[0,j]
lw[j] = 1
w[1,j] = w[0,j]
t = w[2,j]
w[2,j] = w[1,j+1]
w[3,j] = w[2,j+1]
w[1,j+1] = t-w[4,j]*w[2,j]
w[2,j+1] = -w[4,j]*w[3,j]
}
}
if (w[1,nm1] == 0) {
w[1,nm1] = 1e-30
}
if (i > 0 && abs(e[i]-e[i-1]) < eps1) {
for (j = 0; j < n; j++) {
v[i,j] = rand()
}
}
for (;;) {
t = v[i,nm1]
r = v[i,nm2]
v[i,nm1] = t/w[1,nm1]
v[i,nm2] = (r-w[2,nm2]*v[i,nm1])/w[1,nm2]
for (k = nm2-1; k >= 0; k--) {
v[i,k] = (v[i,k]-w[2,k]*v[i,k+1]-w[3,k]*v[i,k+2])/w[1,k]
}
if (sw == 0) {
sw = 1
for (j = 0; j < nm1; j++) {
if (lw[j] == 1) {
t = v[i,j]
v[i,j] = v[i,j+1]
v[i,j+1] = t-w[4,j]*v[i,j+1]
}
else {
v[i,j+1] -= w[4,j]*v[i,j]
}
}
continue
}
else
break
}
}
if (n > 2) {
for (i = 0; i < nm2; i++) {
w[0,i] *= -a[i,i+1]
}
for (i = 0; i < nva; i++) {
for (k = nm2-1; k >= 0; k--) {
r = w[0,k]
if (r != 0) {
r = 1/r
s = 0
for (j = k+1; j < n; j++) {
s += a[k,j]*v[i,j]
}
r *= s
for (j = k+1; j < n; j++) {
v[i,j] -= r*a[k,j]
}
}
}
}
}
for (i = 0; i < nva; i++) {
if (i == 0) {
m = 1
}
else {
if (abs(e[i]-e[i-1]) >= eps1) {
m = i+1
}
else {
for (j = m-1; j < i; j++) {
s = 0
for (k = 0; k < n; k++) {
s += v[j,k]*v[i,k]
}
for (k = 0; k < n; k++) {
v[i,k] -= s*v[j,k]
}
}
}
}
t = s = 0
for (j = 0; j < n; j++) {
s += v[i,j]^2
}
if (s != 0) {
t = sqrt(1/s)
}
for (j = 0; j < n; j++) {
v[i,j] *= t
}
}
}
return 0
}
直前のページへ戻る E-mail to Shigenobu AOKI