メイン 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             偏回帰係数       標準誤差       t値     p値    標準化偏回帰係数\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, "t値の自由度 = ", 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("いずれかの群の有効ケース数が0です。解析は行いません。\n")
    return
  }
  else if (num[0]+num[1] == 2) {
    printf("各群の有効ケース数が1ずつです。解析は行いません。\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("ケース数が1以下です\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

Made with Macintosh