テストデータの生成     Last modified: May 15, 2002

 以下のスクリプトでは、"正確に" 指定された相関係数を持つ人工データを生成します。

 なお、位置および散布度の調整には対応していません。

# 与えられた相関係数行列(因子負荷行列)を再現するデータの生成
# 入力データファイルに相関係数行列(対角成分を含む下側三角行列)または因子負荷行列を用意する(これを仮に 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