メイン mlogist.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="txp.js">document.write("txp.js ファイルが見つかりません??<br>")</script>
<script src="io.js">document.write("io.js ファイルが見つかりません??<br>")</script>
<script language="JavaScript">
<!--
var epsinv = 1e-8
var YES = 1
var NO = 0
function zfi(w, i)
{
var ww = (""+i).length
ww = ww > w ? 0 : w-ww
return "000000000000000".substring(0, ww)+i
}
function vname(i)
{
return "Var"+zfi(2, i+1)
}
function diff(x, m, n, diff1, diff2, coeff)
{
var i, j, k, mp1, p0, p1, pp, temp
mp1 = m+1
for (i = 0; i < mp1; i++) {
diff1[i] = 0
}
for (i = 0; i < mp1; i++) {
for (j = 0; j <= i; j++) {
diff2[j][i] = 0
}
}
for (i = 0; i < n; i++) {
temp = coeff[m]
for (j = 0; j < m; j++) {
temp += coeff[j]*x[i][j]
}
p0 = 1/(1+Math.exp(-temp))
p1 = 1-p0
pp = (x[i][m] == 1) ? p1 : -p0
diff1[m] += pp
diff2[m][m] -= p0*p1
for (j = 0; j < m; j++) {
temp = -x[i][j]*p0*p1
diff1[j] += x[i][j]*pp
for (k = 0; k <= j; k++) {
diff2[k][j] += x[i][k]*temp
}
diff2[j][m] += temp
}
}
for (j = 0; j < mp1; j++) {
for (k = 0; k <= j; k++) {
diff2[j][k] = diff2[k][j]
}
}
}
function llh(x, m, n, coeff)
{
var i, j
var llh_v, temp
llh_v = 0
for (i = 0; i < n; i++) {
temp = coeff[m]
for (j = 0; j < m; j++) {
temp += x[i][j]*coeff[j]
}
if (x[i][m] == 1) {
temp = -temp
}
llh_v -= Math.log(1+Math.exp(temp))
}
return llh_v
}
function newton_logist(x, m, n, sd, coeff, coeff0, diff1, diff2)
{
var mp1, converge, i, itr, k, p, t, temp
mp1 = m+1
for (itr = 0; itr < 500; itr++) {
diff(x, m, n, diff1, diff2, coeff)
if (inv(diff2, mp1, epsinv) == 1) {
converge = NO
return 1
}
for (i = 0; i < mp1; i++) {
temp = 0
for (k = 0; k < mp1; k++) {
temp += diff2[k][i]*diff1[k]
}
coeff0[i] = temp
}
converge = YES
for (i = 0; i < mp1; i++) {
if (Math.abs(coeff0[i]/coeff[i]) > 1e-10) {
converge = NO
}
coeff[i] -= coeff0[i]
}
// for (i = 0; i < m; i++) {
// printf("%9s%24.14g%24.14g\n", vname(i), coeff[i], -coeff0[i])
// }
// printf("%9s%24.14g%24.14g\n\n", "定数項", coeff[m], -coeff0[m])
if (converge == YES) {
break
}
}
if (converge == NO) {
printf("収束計算を途中で打ち切りました。\n")
}
for (i = 0; i < mp1; i++) {
diff1[i] = Math.sqrt(-diff2[i][i])
}
printf("\n対数尤度 = %.14g\n", llh(x, m, n, coeff))
printf("\n係数の推定値\n\n 偏回帰係数 標準誤差 ｔ値 ｐ値 標準化偏回帰係数\n")
for (i = 0; i < m; i++) {
t = Math.abs(coeff[i]/diff1[i])
p = txp(t, n-mp1)
printf("%8s %14.7g %14.7g %12.7g %8.5f %14.7g\n", vname(i), coeff[i], diff1[i], t, p, coeff[i]*sd[i])
}
t = Math.abs(coeff[m]/diff1[m])
p = txp(t, n-mp1)
printf("%8s %14.7g %14.7g %12.7g %8.5f\n\n%60s%i\n", "定数項", coeff[m], diff1[m], t, p, "ｔ値の自由度 = ", n-mp1)
return 0
}
function logist_sub(x, m, n)
{
var mp1, i, j, num, coeff0, coeff, diff1, diff2, mean, sd, xv
mp1 = m--
num = makeVector(2)
mean = makeVector(m)
sd = makeVector(m)
coeff0 = makeVector(mp1)
coeff = makeVector(mp1)
diff1 = makeVector(mp1)
diff2 = makeMatrix(mp1, mp1)
num[0] = num[1] = 0
for (i = 0; i < n; i++) {
j = x[i][m]
num[j]++
}
printf("★ 多重ロジスティックモデル(Walker-Duncan 法) ★\n\n")
printf("有効ケース数 %5i\n", n)
printf(" ネガティブケース %5i\n", num[0])
printf(" ポジティブケース %5i\n", num[1])
if (num[0] == 0 || num[1] == 0) {
printf("いずれかの群の有効ケース数が０です。解析は行いません。\n")
return
}
else if (num[0]+num[1] == 2) {
printf("各群の有効ケース数が１ずつです。解析は行いません。\n")
return
}
for (i = 0; i < mp1; i++) {
coeff[i] = 1e-15
}
for (j = 0; j < m; j++) {
mean[j] = sd[j] = 0
}
for (i = 0; i < n; i++) {
for (j = 0; j < m; j++) {
xv = x[i][j]-mean[j]
mean[j] += xv/(i+1)
sd[j] += i*xv*xv/(i+1)
}
}
for (j = 0; j < m; j++) {
sd[j] = Math.sqrt(sd[j]/(n-1))
}
printf("\n %15s %15s\n", "平均値", "標準偏差")
for (i = 0; i < m; i++) {
printf("%8s %15.7g %15.7g\n", vname(i), mean[i], sd[i])
}
// printf("\n係数の初期値\n\n")
// for (i = 0; i < m; i++) {
// printf("%8s%18.7g\n", vname(i), coeff[i])
// }
// printf(" 定数項%20.7g\n", coeff[m])
if (newton_logist(x, m, n, sd, coeff, coeff0, diff1, diff2) == 1) {
printf("逆行列が求まらず、係数の推定ができませんでした。初期値の設定が妥当でないためと考えられます。\n")
return
}
}
function calc(data_string)
{
var data, nc, nv
if ((data = getdata(data_string, 0)) == false) return
nc = data.length
nv = data[0].length
if (nc < 2) {
printf("ケース数が１以下です\n")
return
}
else if (nc <= nv) {
printf("データの組数は，独立変数の個数より2以上大きくないといけません\n")
return
}
logist_sub(data, nv, nc)
}
//-->
</script>
</head>
<body bgcolor="#ffffff">
<font size="+2"><b>多重ロジスティックモデル(Walker-Duncan 法)</b></font> <a href="src/mlogist.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)"> </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/mlogist.html">使用法</a></td>
</tr></table>
<p>
入力欄<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/Survival/mul-log.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>
サブ inv.js Last modified: Mar 25, 2004
サブ txp.js Last modified: Mar 25, 2004
サブ io.js Last modified: Mar 25, 2004
直前のページへ戻る E-mail to Shigenobu AOKI
