サブ fit.js Last modified: Mar 25, 2004
// S-JIS コードでなくてはならない
var n, ip, loop
var resid, cterm, d, x, y, ff, p, p0, p1, delta, x_value, fj
var MAXIT = 500
var EPSILON = 1e-7
function residual()
{
var i, temp
resid = 0.0
for (i = 0; i < n; i++) {
x_value = x[i]
temp = y[i]-Function(0, i)
resid += temp*temp
}
}
function solinv_fit()
{
var i, j, k, l, nold
var f, g, h, s, sl, sum, temp
for (i = 0; i < ip; i++) {
p[i] = p0[i]
}
residual()
for (i = 0; i < n; i++) {
x_value = x[i]
f = Function(1, i)
fj[i][ip] = y[i]-f
}
for (j = 0; j < ip; j++) {
temp = 0
for (i = 0; i < n; i++) {
temp += fj[i][j]*fj[i][ip]
}
ff[j] = temp
}
nold = n
n += ip
for (j = 0; j < ip; j++) {
for (i = 0; i < ip; i++) {
fj[nold+j][i] = (j == i) ? Math.sqrt(cterm) : 0
}
fj[nold+j][ip] = 0
}
sl = 1e-30
for (j = 0; j <= ip; j++) {
h = 0
for (k = j; k < n; k++) {
d[k] = fj[k][j]
h += d[k]*d[k]
}
if (h < sl) {
g = 0
}
else {
g = Math.sqrt(h)
f = fj[j][j]
if (f >= 0) {
g = -g
}
d[j] = f-g
h -= f*g
for (k = j+1; k <= ip; k++) {
s = 0
for (l = j; l < n; l++) {
s += d[l]*fj[l][k]
}
s /= h
for (l = j; l < n; l++) {
fj[l][k] -= d[l]*s
}
}
}
fj[j][j] = g
}
n = nold
for (i = ip-1; i >= 0; i--) {
sum = fj[i][ip]
for (j = i+1; j < ip; j++) {
sum -= delta[j]*fj[i][j]
}
delta[i] = (fj[i][i] == 0) ? 0 : sum/fj[i][i]
}
}
function output(title, method)
{
var yval, i
var p_name
p_name = new Array("a", "b", "c", "d", "e", "f", "g", "h", "i", "j")
printf("%s %s 法\n\n", title, method == 1 ? "Marquardt" : "Simplex")
for (i = 0; i < ip; i++) {
printf("%s ・・・・・・・・・・・・・%22.14g\n", p_name[i], p[i])
}
printf("残差平方和 ・・・・%22.14g\n\n", resid)
printf("★ 予測値および残差\n\nケース 独立変数 従属変数 予測値 残差\n")
for (i = 0; i < n; i++) {
x_value = x[i]
yval = Function(0, i)
printf("%6i%15.7g%15.7g%15.7g%15.7g\n",i+1, x_value, y[i], yval, y[i]-yval)
}
}
function marquardt(ip0, inival, title)
{
var errset = 0, i, diverg = 0, flag = 0
var check, check0, rchange, res0, res1
ip = ip0
fj = makeMatrix(n+ip, ip)
ff = makeVector(ip)
p = makeVector(ip)
p0 = makeVector(ip)
p1 = makeVector(ip)
delta = makeVector(ip)
for (i = 0; i < ip; i++) {
p0[i] = inival[i]
}
d = makeVector(n)
document.Result.result.value = "非線形最小二乗法によるあてはめ収束状況"+newline
cterm = 100
for (loop = 1; loop <= MAXIT; loop++) {
solinv_fit()
res0 = resid
flag = 0
for (i = 0; i < ip; i++) {
p1[i] = p0[i]+delta[i]
if (Math.abs(delta[i]/p1[i]) > EPSILON) {
flag = 1
}
}
check = 0
for (i = 0; i < ip; i++) {
check -= (ff[i]+cterm*delta[i])*delta[i]
}
check0 = check
for (i = 0; i < ip; i++) {
p[i] = p1[i]
}
residual()
printf("residual = %g\n", resid)
res1 = resid
if (res1 == 0) {
break
}
rchange = res1-res0
if (rchange <= 0) {
diverg = 0
if (errset == 1) {
errset = 0
}
check = rchange/check0
if (-rchange/res1 < EPSILON && flag == 0) {
break
}
if (check > 0.75) {
cterm *= 0.5
}
else if (check < 0.25) {
cterm *= 5
}
for (i = 0; i < ip; i++) {
p0[i] = p1[i]
}
}
else {
diverg++
errset = 1
if (diverg > 10) {
if (flag == 1) {
printf("\n発散してしまいました。初期値が不適当です。\n")
break
}
for (i = 0; i < ip; i++) {
p0[i] = p1[i]
}
}
else {
cterm *= 5
loop--
}
}
}
document.Result.result.value = "非線形最小二乗法によるあてはめ結果"+newline
for (i = 0; i < ip; i++) {
p[i] = p1[i]
}
if (loop > MAXIT || flag) {
printf("\n収束しませんでした！初期値が不適当です。\n")
}
output(title, 1)
}
function simplex(ip0, inival, title)
{
var i, j, l1, l2, l5, loop0, m, mi, mx, fi9, fj9, pa, res, resid_old, s, s1, s2, z1, z2
document.Result.result.value = "非線形最小二乗法によるあてはめ収束状況"+newline
ip = ip0
p0 = makeVector(ip)
for (i = 0; i < ip; i++) {
p0[i] = inival[i]
}
pa = makeMatrix(ip+3, ip+3)
res = makeVector(ip+3)
s = makeVector(ip)
p = makeVector(ip)
resid_old = 1e33
delta = new Array(ip+1)
for (i = 0; i < ip; i++) {
delta[i] = p0[i]*(Math.random()*0.1-0.05)
}
l5 = ip+1
l1 = l5+1
l2 = l5+2
fi9 = l5
fj9 = ip
s1 = 2/fj9
s2 = 3/fj9
z1 = (fi9+1)/fj9
z2 = (2*fi9+1)/fj9
for (i = 0; i < l5; i++) {
for (j = 0; j < ip; j++) {
pa[j][i] = p0[j]
}
if (i <= ip) {
pa[i][i] = p0[i]+delta[i]
}
for (m = 0; m < ip; m++) {
p[m] = pa[m][i]
}
residual()
printf("residual = %g\n", resid)
res[i] = resid
}
loop0 = 0
for (loop = 1; loop <= MAXIT; loop++) {
for (j = 0; j < ip; j++) {
s[j] = 0
}
mx = mi = 0
for (i = 0; i < l5; i++) {
if (res[i] > res[mx]) {
mx = i
}
if (res[i] < res[mi]) {
mi = i
}
for (j = 0; j < ip; j++) {
s[j] += pa[j][i]
}
}
if (++loop0 > MAXIT) {
break
}
if (res[mx] < EPSILON || res[mi] < EPSILON) {
break
}
if ((res[mx]-res[mi])/res[mi] <= EPSILON) {
break
}
loop0 = 0
i = l1
for (j = 0; j < ip; j++) {
pa[j][l1] = s1*s[j]-z1*pa[j][mx]
}
for (m = 0; m < ip; m++) {
p[m] = pa[m][l1]
}
residual()
printf("residual = %g\n", resid)
res[l1] = resid
if (res[l1] < res[mi]) {
for (j = 0; j < ip; j++) {
pa[j][l2] = s2*s[j]-z2*pa[j][mx]
}
for (m = 0; m < ip; m++) {
p[m] = pa[m][l2]
}
residual()
printf("residual = %g\n", resid)
res[l2] = resid
if (res[l2] <= res[l1]) {
i = l2
}
}
else if (res[l1] > res[mx]) {
for (j = 0; j < ip; j++) {
pa[j][l2] = s[j]/fi9
}
for (m = 0; m < ip; m++) {
p[m] = pa[m][l2]
}
residual()
printf("residual = %g\n", resid)
res[l2] = resid
if (res[l2] >= res[mx]) {
for (i = 0; i < l5; i++) {
if (i != mi) {
for (j = 0; j < ip; j++) {
pa[j][i] = (pa[j][i]+pa[j][mi])*0.5
}
for (m = 0; m < ip; m++) {
p[m] = pa[m][i]
}
residual()
printf("residual = %g\n", resid)
res[i] = resid
}
}
continue
}
i = l2
}
for (j = 0; j < ip; j++) {
pa[j][mx] = pa[j][i]
}
res[mx] = res[i]
}
document.Result.result.value = "非線形最小二乗法によるあてはめ結果"+newline
if (loop0 > MAXIT || loop > MAXIT) {
printf("収束しませんでした。\n")
}
for (j = 0; j < ip; j++) {
p[j] = pa[j][mi]
}
resid = res[mi]
output(title, 2)
}
直前のページへ戻る E-mail to Shigenobu AOKI
