メイン svdcmp.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="io.js">document.write("io.js ファイルが見つかりません??<br>")</script>
<script language="JavaScript">
<!--
function square(a)
{
return a == 0 ? 0 : a*a
}
function pythag(a, b)
{
var absa, absb
absa = Math.abs(a)
absb = Math.abs(b)
if (absa > absb) {
return absa * Math.sqrt(1 + square(absb / absa))
}
else if (absb == 0) {
return 0
}
else {
return absb * Math.sqrt(1 + square(absa / absb))
}
}
function sign(a, b)
{
return b >= 0 ? Math.abs(a) : -Math.abs(a)
}
// 以下のプログラム svdcmp は，
// William H. Press, Saul A. Teukolski, William T. Vetterling and Bryan P. Flannery
// Numerical Recipes in C
// Cambridge University Press 1992
// の 67-70 ページにある svdcmp を JavaScript に書き直したものである S. Aoki 2001/09/15
function svdcmp(a, n, m, w, v)
{
var flag, i, its, j, jj, k, l, nm
var anorm, c, f, g, h, s, scale, x, y, z, rv1
rv1 = makeVector(m)
g = scale = anorm = 0
for (i = 0; i < m; i++) {
l = i + 1
rv1[i] = scale * g
g = s = scale = 0
if (i < n) {
for (k = i; k < n; k++) {
scale += Math.abs(a[k][i])
}
if (scale) {
for (k = i; k < n; k++) {
a[k][i] /= scale
s += a[k][i] * a[k][i]
}
f = a[i][i]
g = -sign(Math.sqrt(s), f)
h = f * g - s
a[i][i] = f - g
for (j = l; j < m; j++) {
s = 0
for (k = i; k < n; k++) {
s += a[k][i] * a[k][j]
}
f = s / h
for (k = i; k < n; k++) {
a[k][j] += f * a[k][i]
}
}
for (k = i; k < n; k++) {
a[k][i] *= scale
}
}
}
w[i] = scale * g
g = s = scale = 0
if (i < n && i != m-1) {
for (k = l; k < m; k++) {
scale += Math.abs(a[i][k])
}
if (scale) {
for (k = l; k < m; k++) {
a[i][k] /= scale
s += a[i][k] * a[i][k]
}
f = a[i][l]
g = -sign(Math.sqrt(s), f)
h = f * g - s
a[i][l] = f - g
for (k = l; k < m; k++) {
rv1[k] = a[i][k] / h
}
for (j = l; j < n; j++) {
s = 0
for (k = l; k < m; k++) {
s += a[j][k] * a[i][k]
}
for (k = l; k < m; k++) {
a[j][k] += s * rv1[k]
}
}
for (k = l; k < m; k++) {
a[i][k] *= scale
}
}
}
anorm = Math.max(anorm, (Math.abs(w[i]) + Math.abs(rv1[i])))
}
for (i = m-1; i >= 0; i--) {
if (i < m-1) {
if (g) {
for (j = l; j < m; j++) {
v[j][i] = (a[i][j] / a[i][l]) / g
}
for (j = l; j < m; j++) {
s = 0
for ( k = l; k < m; k++) {
s += a[i][k] * v[k][j]
}
for (k = l; k < m; k++) {
v[k][j] += s * v[k][i]
}
}
}
for (j = l; j < m; j++) {
v[i][j] = v[j][i] = 0
}
}
v[i][i] = 1
g = rv1[i]
l = i
}
for (i = Math.min(n, m)-1; i >= 0; i--) {
l = i + 1
g = w[i]
for (j = l; j < m; j++) {
a[i][j] = 0
}
if (g) {
g = 1 / g
for (j = l; j < m; j++) {
s = 0
for (k = l; k < n; k++) {
s += a[k][i] * a[k][j]
}
f = (s / a[i][i]) * g
for (k = i; k < n; k++) {
a[k][j] += f * a[k][i]
}
}
for (j = i; j < n; j++) {
a[j][i] *= g
}
}
else {
for (j = i; j < n; j++) {
a[j][i] = 0
}
}
++a[i][i]
}
for (k = m-1; k >= 0; k--) {
for (its = 1; its <= 30; its++) {
flag = 1
for (l = k; l >= 0; l--) {
nm = l - 1
if (Math.abs(rv1[l]) + anorm == anorm) {
flag = 0
break
}
if (Math.abs(w[nm]) + anorm == anorm)
break
}
if (flag) {
c = 0
s = 1
for (i = l; i < k; i++) {
f = s * rv1[i]
rv1[i] = c * rv1[i]
if (Math.abs(f) + anorm == anorm) {
break
}
g = w[i]
h = pythag(f, g)
w[i] = h
h = 1 / h
c = g * h
s = -f * h
for (j = 0; j < n; j++) {
y = a[j][nm]
z = a[j][i]
a[j][nm] = y * c + z * s
a[j][i] = z * c - y * s
}
}
}
z = w[k]
if (l == k) {
if (z < 0) {
w[k] = -z
for (j = 0; j < m; j++) {
v[j][k] = -v[j][k]
}
}
break
}
if (its == 30) {
printf("no convergence in 30 svdcmp iterations\n")
}
x = w[l]
nm = k - 1
y = w[nm]
g = rv1[nm]
h = rv1[k]
f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2 * h * y)
g = pythag(f, 1)
f = ((x - z) * (x + z) + h * ((y / (f + sign(g, f))) - h)) / x
c = s = 1
for (j = l; j <= nm; j++) {
i = j + 1
g = rv1[i]
y = w[i]
h = s * g
g = c * g
z = pythag(f, h)
rv1[j] = z
c = f / z
s = h / z
f = x * c + g * s
g = g * c - x * s
h = y * s
y *= c
for (jj = 0; jj < m; jj++) {
x = v[jj][j]
z = v[jj][i]
v[jj][j] = x * c + z * s
v[jj][i] = z * c - x * s
}
z = pythag(f, h)
w[j] = z
if (z) {
z = 1 / z
c = f * z
s = h * z
}
f = c * g + s * y
x = c * y - s * g
for (jj = 0; jj < n; jj++) {
y = a[jj][j]
z = a[jj][i]
a[jj][j] = y * c + z * s
a[jj][i] = z * c - y * s
}
}
rv1[l] = 0
rv1[k] = f
w[k] = x
}
}
}
function outmat(r, n, m, str)
{
var j, k
printf("\n%s\n\n", str)
for (j = 0; j < n; j++) {
for (k = 0; k < m; k++) {
printf(" %12.6g", r[j][k])
}
printf("\n")
}
}
function sort(a, diag, v, n, m)
{
var i, j, maxi, maxv, temp
for (i = 0; i < m-1; i++) {
maxi = i
maxv = diag[i]
for (j = i+1; j < m; j++) {
if (diag[j] > maxv) {
maxi = j
maxv = diag[j]
}
}
if (maxi != i) {
diag[maxi] = diag[i]
diag[i] = maxv
for (j = 0; j < n; j++) {
temp = a[j][i]
a[j][i] = a[j][maxi]
a[j][maxi] = temp
}
for (j = 0; j < m; j++) {
temp = v[j][i]
v[j][i] = v[j][maxi]
v[j][maxi] = temp
}
}
}
}
function calc(data_string)
{
var a, n, m, diag, v
var debug = false
if ((a = getdata(data_string, 0)) != false) {
if ((n = a.length) == 0) {
printf("データが入力されていません\n")
return
}
m = a[0].length
diag = makeVector(m)
v = makeMatrix(m, m)
outmat(a, n, m, "★ 入力された行列 ★")
if (svdcmp(a, n, m, diag, v) == 1) {
printf("特異値分解ができませんでした\n")
}
else {
sort(a, diag, v, n, m)
sum = 0
for (i = 0; i < m; i++) {
sum += diag[i]
}
cum = 0
printf("\n%4s %10s %12s %12s\n\n", "番号", "固有値", "寄与率", "累積寄与率")
for (i = 0; i < m; i++) {
cum += diag[i]/sum*100
printf("%4i %10.6g %11.1f%% %11.1f%%\n", i+1, diag[i], diag[i]/sum*100, cum)
}
for (i = 0; i < n; i++) {
for (j = 0; j < m; j++) {
a[i][j] *= diag[j]
}
}
outmat(a, n, m, 'ＵΛ')
outmat(v, m, m, 'Ｖ')
if (debug) {
printf("\n検算\n\n")
for (i = 0; i < n; i++) {
for (j = 0; j < m; j++) {
sum = 0
for (k = 0; k < m; k++) {
sum += a[i][k]*v[j][k]
}
printf("%10.5g", sum)
}
printf("\n")
}
}
}
}
sep(40)
}
//-->
</script>
</head>
<body bgcolor="#ffffff">
<font size="+2"><b>特異値分解</b></font> <a href="src/svdcmp.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>
<input type="button" name="calcurate" value="計算開始" onClick="calc(this.form.data.value)">
<input type="button" name="clear" value="入力欄クリア" onClick="this.form.data.value=''">
<input type="button" name="clear" value="出力欄クリア" onClick="this.form.result.value=''">
<img src="../gra/button3.png" width=9 height=9 alt="・"> <A HREF="exa/svdcmp.html">使用法</a><p>
入力欄（データ行列を<a href="exa/kugirimoji.html">区切り文字</a>で区切って，行単位に入力）<br><textarea name="data" ROWS=10 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="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>
サブ io.js Last modified: Mar 25, 2004
直前のページへ戻る E-mail to Shigenobu AOKI
