メイン pfa.html Last modified: Sep 01, 2009
<html>
<head>
<meta http-equiv="Content-Type" content="text/html;CHARSET=EUC-JP">
<link rel="shortcut icon" href="../favicon.ico">
<title>JavaScript</title>
<script src="inv.js">document.write("inv.js ファイルが見つかりません??<br>")</script>
<script src="eigen.js">document.write("eigen.js ファイルが見つかりません??<br>")</script>
<script src="io.js">document.write("io.js ファイルが見つかりません??<br>")</script>
<script language="JavaScript">
<!--
var maxrot = 500
var ne
function zfi(w, i)
{
var ww = (""+i).length
ww = ww > w ? 0 : w-ww
return "0000000000000".substring(0, ww)+i
}
function fname(i)
{
return "PCA"+zfi(2, i+1)
}
function vname(i)
{
return "Var"+zfi(2, i+1)
}
function facnumber(i)
{
return "第"+zfi(2, i+1)+"因子"
}
function or_r(dest, source, m)
{
var i, j
for (i = 0; i < m; i++) {
for (j = 0; j < m; j++) {
dest[i][j] = source[i][j]
}
}
}
function n_factor(m, eigen_value)
{
var i
if (isNaN(ne) || ne < 1 || ne > m) {
ne = 0
for (i = 0; i < m; i++) {
if (eigen_value[i] >= 1.0) {
ne++
}
}
}
}
function line_out(side, x, ne)
{
var i
printf("%12s", side)
for (i = 0; i < ne; i++) {
printf("%11.5f", x[i])
}
printf("\n")
}
function out_loadings(m, ne, fl, com, eig)
{
var i, j, x0, x1, temp2 = 0.0
x0 = makeVector(m)
x1 = makeVector(m)
for (i = 0; i < ne; i++) {
x1[i] = temp2 += x0[i] = eig[i]/m*100.0
}
printf(" 変数")
for (i = 0; i < ne; i++) {
printf("%11s", facnumber(i))
}
printf(" 共通性\n")
for (i = 0; i < m; i++) {
printf(" %5s", vname(i))
for (j = 0; j < ne; j++) {
printf("%11.5f", fl[j][i])
}
printf("%11.5f\n", com[i])
}
line_out("因子負荷量の二乗和", eig, ne)
line_out(" 寄与率％", x0, ne)
line_out(" 累積寄与率％", x1, ne)
}
function pfa(nv, r, fl, or, com, ocom, oocom, eigen_value)
{
var i, irotat, isw, j, jsw, ne1, difmax, eig, xx, temp
or_r(or, r, nv)
eigen(r, nv, nv, 1e-10, eigen_value, fl)
printf("\n初期固有値\n\n%4s %10s %10s\n", "番号", "固有値", "累積パーセント")
temp = 0
for (i = 0; i < nv; i++) {
temp += eigen_value[i]
printf("%4i %10.5f %10.2f\n", i+1, eigen_value[i], temp/nv*100)
}
n_factor(nv, eigen_value)
or_r(r, or, nv)
or_r(fl, or, nv)
if (inv(fl, nv, 1e-8) == 0) {
printf("\n対角要素にＳＭＣを使いました\n")
for (i = 0; i < nv; i++) {
r[i][i] = 1.0-1.0/fl[i][i]
}
}
else {
printf("\n対角要素に相関係数の最大値を使いました\n")
for (i = 0; i < nv; i++) {
xx = 0.0
for (j = 0; j < nv; j++) {
if (Math.abs(r[j][i]) > xx && i != j) {
xx = Math.abs(r[j][i])
}
}
r[i][i] = xx
}
}
for (i = 0; i < nv; i++) {
ocom[i] = 1.0
}
eigen(r, nv, ne, 1e-10, eigen_value, fl)
ne1 = ne
for (i = 0; i < nv; i++) {
if (eigen_value[i] <= 0.0) {
ne1 = i
break
}
}
if (ne1 < ne) {
printf("最初の共通性の推定で、因子数が縮小されました\n")
ne = ne1
}
irotat = 0
jsw = 0
for (irotat = 1; irotat <= maxrot+1; irotat++) {
var temp9
for (i = 0; i < nv; i++) {
com[i] = 0.0
}
for (i = 0; i < ne; i++) {
if (eigen_value[i] < 0) {
eigen_value[i] = 0.0
}
eig = Math.sqrt(eigen_value[i])
for (j = 0; j < nv; j++) {
fl[i][j] *= eig
com[j] += fl[i][j]*fl[i][j]
}
}
difmax = 0.0
isw = 0
for (i = 0; i < nv; i++) {
temp9 = Math.abs(com[i]-ocom[i])
if ((difmax = Math.max(temp9, difmax)) > 1e-10) {
isw = 1
}
}
if (irotat > maxrot) {
printf("共通性の推定が規定回数内にできませんでした\n")
jsw = 1
}
if (isw == 0 || jsw == 1) {
break
}
for (i = 0; i < nv; i++) {
if (com[i] > 1.0) {
jsw = 1
break
}
}
if (jsw == 1) {
printf("\n共通性が１を越えるものがあったので，直前までの解を表示します\n")
if (irotat == 1) {
printf("因子分析は行えません 抽出する因子数を少なくして再実行してみてください\n")
}
for (i = 0; i < nv; i++) {
com[i] = oocom[i]
}
}
for (i = 0; i < nv; i++) {
for (j = 0; j < i; j++) {
r[i][j] = r[j][i] = or[j][i]
}
r[i][i] = com[i]
oocom[i] = ocom[i]
ocom[i] = com[i]
}
eigen(r, nv, ne, 1e-10, eigen_value, fl)
}
if (jsw == 1) {
irotat--
}
printf("\n反復推定回数 ・・・・・・ %i\n\n", irotat)
printf("回転前の因子負荷量\n\n")
out_loadings(nv, ne, fl, com, eigen_value)
}
function rotat_fa(m, ne, vmxeps, eig, com, fl)
{
var ov = 0
var fnp, fljk, fljk1, theta, v, v0, x, xx, xy, y, yy, a, b, c, cs, d, dd, emax, si, temp, work, weight
var i, j, k, k1, irotat, mmax
weight = 1.0
fnp = m
for (i = 0; i < ne; i++) {
eig[i] = 0.0
}
for (i = 0; i < m; i++) {
for (j = 0; j < ne; j++) {
fl[j][i] /= Math.sqrt(com[i])
}
}
for (irotat = 1; irotat <= 9999; irotat++) {
for (k = 0; k < ne-1; k++) {
for (k1 = k+1; k1 < ne; k1++) {
a = b = c = d = 0.0
for (j = 0; j < m; j++) {
x = fl[k][j]
y = fl[k1][j]
xx = x*x
yy = y*y
xy = xx-yy
a += xy
b += 2.0*x*y
c += xy*xy-4.0*xx*yy
d += 4.0*x*y*xy
}
dd = d-2.0*weight*a*b/fnp
theta = Math.atan(dd/(c-weight*(a*a-b*b)/fnp))
if (Math.sin(theta)*dd < 0.0) {
if (theta > 0.0) {
theta -= Math.PI
}
else {
theta += Math.PI
}
}
theta *= 0.25
cs = Math.cos(theta)
si = Math.sin(theta)
for (j = 0; j < m; j++) {
fljk = fl[k][j]
fljk1 = fl[k1][j]
fl[k][j] = fljk*cs+fljk1*si
fl[k1][j] = fljk1*cs-fljk*si
}
}
}
v = 0.0
for (k = 0; k < ne; k++) {
v0 = 0.0
for (j = 0; j < m; j++) {
v0 += fl[k][j]*fl[k][j]
}
v0 *= weight/fnp
for (j = 0; j < m; j++) {
work = fl[k][j]*fl[k][j]-v0
v += work*work
}
}
if (Math.abs(v-ov) < vmxeps) {
break
}
ov = v
}
for (i = 0; i < m; i++) {
temp = Math.sqrt(com[i])
for (j = 0; j < ne; j++) {
fl[j][i] *= temp
eig[j] += fl[j][i]*fl[j][i]
}
}
for (i = 0; i < ne-1; i++) {
emax = eig[i]
mmax = i
for (j = i+1; j < ne; j++) {
if (eig[j] > emax) {
emax = eig[j]
mmax = j
}
}
if (mmax != i) {
var temp_p
eig[mmax] = eig[i]
eig[i] = emax
temp_p = fl[i]; /* ポインタ交換ですませられる */
fl[i] = fl[mmax]
fl[mmax] = temp_p
}
}
}
function out_coef(m, ne, coef)
{
var i, j
printf("\n因子得点係数\n\n")
printf(" 変数")
for (i = 0; i < ne; i++) {
printf("%11s", facnumber(i))
}
printf("\n")
for (i = 0; i < m; i++) {
printf("%6s", vname(i))
for (j = 0; j < ne; j++) {
printf("%11.5f", coef[j][i])
}
printf("\n")
}
}
function facscore(m, nc, data, mean, sd, ne, coef)
{
var i, j, k, x0, temp
x0 = makeVector(m)
for (k = 0; k < m; k++) {
sd[k] = sd[k]*Math.sqrt((nc-1)/nc)
}
printf("\n因子得点\n\n")
printf("ケース")
for (i = 0; i < ne; i++) {
printf("%11s", facnumber(i))
}
printf("\n")
for (i = 0; i < nc; i++) {
for (k = 0; k < m; k++) {
x0[k] = (data[i][k]-mean[k])/sd[k]
}
printf("%6i", i+1)
for (j = 0; j < ne; j++) {
temp = 0.0
for (k = 0; k < m; k++) {
temp +=coef[j][k]*x0[k]
}
printf("%11.5f", temp)
}
printf("\n")
}
}
function rotate(m, ne, eig, com, fl, or, coef)
{
var i, j, l
if (ne > 1) {
rotat_fa(m, ne, 1e-5, eig, com, fl)
printf("\nバリマックス回転後の因子負荷量\n\n")
out_loadings(m, ne, fl, com, eig)
}
if (inv(or, m, 1e-8) == 1) {
printf("逆行列が求まらないため、因子得点も求まりません。\n")
return
}
for (i = 0; i < m; i++) {
for (j = 0; j < ne; j++) {
var sum = 0.0
for (l = 0; l < m; l++) {
sum += or[l][i]*fl[j][l]
}
coef[j][i] = sum
}
}
out_coef(m, ne, coef)
}
function calc(data_string, ne0)
{
var data, nc, nv, i, j, k, mean, variance, sd, r, or, isw, coef, fl, com, ocom, oocom, eigen_value
if ((data = getdata(data_string, 0)) == false) return
nc = data.length
nv = data[0].length
if (nc < 2) {
printf("ケース数が１以下です\n")
return
}
ne = parseInt(ne0)
mean = new Array(nv)
variance = new Array(nv)
sd = new Array(nv)
r = makeMatrix(nv, nv)
or = makeMatrix(nv, nv)
coef = makeMatrix(nv, nv)
// Netscape communicator のバグ？配列ベクトルは初期化されていない
for (i = 0; i < nv; i++) {
mean[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] += data[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 < nv; k++) {
r[j][k] += (mean[j]-data[i][j])*(mean[k]-data[i][k])
}
}
}
isw = 0
for (j = 0; j < nv; j++) {
variance[j] = r[j][j]/(nc-1)
sd[j] = Math.sqrt(variance[j])
if (variance[j] <= 0) isw = 1
}
if (isw == 1) {
printf("分散が０になる変数があります\n")
return
}
for (j = 0; j < nv; j++) {
for (k = 0; k < nv; k++) {
r[j][k] /= (nc-1)*sd[j]*sd[k]
}
}
printf("★ 因子分析 ★\n\n")
printf("有効データ組数 = %i\n\n", nc)
printf("%5s%15s%15s%15s\n", "変数", "平均値", "不偏分散", "標準偏差")
for (j = 0; j < nv; j++) {
printf("%5s%15.7g%15.7g%15.7g\n", vname(j), mean[j], variance[j], sd[j])
}
printf("\n相関係数行列\n\n ")
for (j = 0; j < nv; j++) {
printf("%8s", vname(j))
}
printf("\n")
for (j = 0; j < nv; j++) {
printf("%5s", vname(j))
for (k = 0; k < nv; k++) {
printf("%8.3f", r[j][k])
}
printf("\n")
}
eigen_value = makeMatrix(nv, nv)
fl = makeMatrix(nv, nv)
com = new Array(nv)
ocom = new Array(nv)
oocom = new Array(nv)
pfa(nv, r, fl, or, com, ocom, oocom, eigen_value)
rotate(nv, ne, eigen_value, com, fl, or, coef)
facscore(nv, nc, data, mean, sd, ne, coef)
}
//-->
</script>
</head>
<body bgcolor="#ffffff">
<font size="+2"><b>因子分析</b></font> <a href="src/pfa.html"><img src="png/src.png" width=35 height=11 alt="src" align=top></a> Last modified: Jun 01, 2006<hr noshade><p>
<font color="#ff0000" size="+2">以下のプログラムのサポートは終了しました。自己責任でお使い下さい。</font>
<form name=Result>
<script language="JavaScript">
<!--
//-->JavaScript がサポートされていないブラウザですか？
</script>
<table><tr>
<td><input type="button" name="calcurate" value="計算開始" onClick="calc(this.form.data.value, this.form.ne.value)"> </td>
<td><input type="button" name="clear" value="入力欄クリア" onClick="this.form.data.value=''"> </td>
<td><input type="button" name="clear" value="出力欄クリア" onClick="this.form.result.value=''"> </td>
<td nowrap><img src="../gra/button3.png" width=9 height=9 alt="・"> <a href="exa/corr.html">使用法</a></td>
</tr></table>
<p>
抽出する因子の個数<input name="ne" value="" size=5><br>
入力欄<br><textarea name="data" ROWS=20 COLS=80></textarea><p>
出力欄<br><textarea name="result" ROWS=30 COLS=80></textarea>
</form>
<p><hr noshade>
<img src="../gra/button3.png" width=9 height=9 alt="・"> <a href="../lecture/PFA/index.html">手法の解説</a><br>
<img src="../gra/button3.png" width=9 height=9 alt="・"> <A HREF="javascript:history.go(-1)">直前のページへ戻る</A> <img src="../gra/button3.png" width=9 height=9 alt="・"> <a href="../mail.html">E-mail to Shigenobu AOKI</a>
<p><center><IMG SRC="../gra/ume5.png" width=121 height=37 ALT="Made with Macintosh"></center>
</body>
</html>
サブ eigen.js Last modified: Mar 25, 2004
サブ inv.js Last modified: Mar 25, 2004
サブ io.js Last modified: Mar 25, 2004
直前のページへ戻る E-mail to Shigenobu AOKI
