メイン kw-exact.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="gxp.js">document.write("gxp.js ファイルが見つかりません??<br>")</script>
<script src="xxp.js">document.write("xxp.js ファイルが見つかりません??<br>")</script>
<script src="io.js">document.write("io.js ファイルが見つかりません??<br>")</script>

<script language="JavaScript">
<!--
//----------------------------------
var nr, nc, um, perm_table, t, nntrue, nntrue2, log_expmax, ntables, stat_val, nc1, nr1, gamma2, half, q, rank_sum
var EXPMAX = 1e100

function zfi(w, i)
{
  var ww = (""+i).length
  ww = ww > w ? 0 : w-ww
  return "000000000000000".substring(0, ww)+i
}

function vname(i, str)
{
  return "第"+zfi(2, i+1)+str
}
//----------------------------------
function permutation(n)
{
  var i
  perm_table = new Array(n+1)
  perm_table[0] = perm_table[1] = 0.0
  for (i = 2; i <= n; i++) {
    perm_table[i] = perm_table[i-1]+Math.log(i)
  }
}

function found()
{
  var i, j
  var hh = kw_h(um)
  if (hh >= stat_val) {
    var nprod = 0
    for (i = 0; i < nr; i++) {
      nprod += perm_table[t[i][nc]]
      for (j = 0; j < nc; j++) {
        nprod -= perm_table[um[i][j]]
      }
    }
      nntrue += Math.exp(nprod-nntrue2*log_expmax)
      while (nntrue >= EXPMAX) {
      nntrue /= EXPMAX
      nntrue2++
    }
  }
  ntables++
}

function search(x, y)
{
  var t
  if (y == 0) {
    found()
  }
  else if (x == 0) {
    if ((t = um[0][0]-um[y][0]) < 0) {
      return
    }
    um[0][0] = t
    search(nc1, y-1)
    um[0][0] += um[y][0]
  }
  else {
    search(x-1, y)
    while (um[y][0] && um[0][x]) {
      ++um[y][x]
      --um[y][0]
      --um[0][x]
      search(x-1, y)
    }
    um[y][0] += um[y][x]
    um[0][x] +=um[y][x]
    um[y][x] = 0
  }
}

function exact()
{
  var i, j, ntot, denom, denom2

  ntables = nntrue = nntrue2 = 0
  ntot = t[nr][nc]
  permutation(ntot)
  denom2 = 0
  log_expmax = Math.log(EXPMAX)
  denom = perm_table[ntot]
  for (j = 0; j < nc; j++) {
    denom -= perm_table[t[nr][j]]
  }
  while (denom > log_expmax) {
    denom -= log_expmax
    denom2++
  }
  denom = Math.exp(denom)

  um = new Array(nr)
  for (i = 0; i < nr; i++) {
    um[i] = new Array(nc)
    for (j = 0; j < nc; j++) {
      um[i][j] = 0
    }
  }

// 列和
  for (i = 0; i < nc; i++) {
    um[0][i] = t[nr][i]
  }

// 行和  2...k
  for (i = 1; i < nr; i++) {
    um[i][0] =t[i][nc]
  }

  nc1 = nc-1
  nr1 = nr-1

  search(nc1, nr1)

  printf("正確な P 値 = %.10g\n", nntrue / denom * Math.pow(EXPMAX, nntrue2 - denom2))
  printf("査察した分割表の個数は %s 個\n", ntables)
}
//----------------------------------
function marginals()
{
  var i, j
  for (i = 0; i < nr; i++) {
    for (j = 0; j < nc; j++) {
      t[i][nc] += t[i][j]
    }
    for (j = 0; j <= nc; j++) {
      t[nr][j] += t[i][j]
    }
  }
}

function kw_h(u)
{
  var i, j
  var r, retv = 0.0
  for (i = 0; i < nr; i++) {
    r = -t[i][nc]*half
    for (j = 0; j < nc; j++) {
      r += rank_sum[j][u[i][j]]
    }
    retv += r*r/t[i][nc]
  }
  return retv
}

function asymptotic()
{
  var i, j, df, p

  printf("★ Kruskal-Wallis の検定(Exact test) ★\n      ")
  for (j = 0; j < nc; j++) {
    printf(" %6s", vname(j, "列"))
  }
  printf(" %6s\n", "合計")
  for (i = 0; i < nr; i++) {
    printf("%6s", vname(i, "群"))
    for (j = 0; j < nc; j++) {
      printf(" %6i", t[i][j])
    }
    printf(" %6i\n", t[i][j])
  }
  printf("%6s", "合計 ")
  for (j = 0; j < nc; j++) {
    printf(" %6i", t[i][j])
  }
  printf(" %6i\n\n", t[i][j])

  stat_val = kw_h(t)
  p = xxp(stat_val*gamma2, nr-1)
  printf(" カイ二乗値 = %.5f\n", stat_val*gamma2)
  printf("     自由度 = %i\n", nr-1)
  printf("       P 値 = %.5f\n\n", p)

}
//----------------------------------

function pow3(x)
{
  return x*x*x
}

function calc(data_string)
{
  var i, j, k, t0, ntot, tie

  if ((t0 = getdata(data_string, 0)) != false) {
  
    if ((nr = t0.length) == 0) {
      printf("分割表が入力されていません\n")
      return
    }
    nc = t0[0].length
    if (nr < 2 || nc < 2) {
      printf("2行2列以上の大きさの分割表を対象にします\n")
      return
    }
    t = new Array(nr+1)
    for (i = 0; i <= nr; i++) {
      t[i] = new Array(nc+1)
      for (j = 0; j <= nc; j++) {
        t[i][j] = 0
      }
    }
    for (i = 0; i < nr; i++) {
      for (j = 0; j < nc; j++) {
        t[i][j] = t0[i][j]
      }
    }
    marginals()
    q = new Array(nc)
    rank_sum = new Array(nc)
    gamma2 = tie = 0
    ntot = t[nr][nc]
    for (j = 0; j < nc; j++) {
      q[j] = tie+(t[nr][j]+1)*0.5
      gamma2 +=pow3(t[nr][j])-t[nr][j]
      tie += t[nr][j]
      rank_sum[j] = new Array(t[nr][j]+1)
      for (k = 0; k <= t[nr][j]; k++) {
        rank_sum[j][k] = k*q[j]
      }
    }
    gamma2 = 12/(ntot*(ntot+1)*(1-gamma2/(pow3(ntot)-ntot)))
    half = (ntot+1)*0.5
    asymptotic()
    exact()
  }
  sep(40)
}
//-->
</script>
</head>

<body bgcolor="#ffffff">
<font size="+2"><b>分割表形式でのクラスカル・ウォリスの検定(Exact test)</b></font> <a  href="src/kw-exact.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>
注意:計算時間が長くなる場合があります(データの個数などからは事前に予測できません)。<p>
<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/kw-exact.html">使用法</a><p>

<table>
<tr>
<td>入力欄<br><textarea name="data" ROWS=20 COLS=15></textarea></td>
<td>出力欄<br><textarea name="result" ROWS=20 COLS=80></textarea></td>
</tr>
</table>
</form>

<p><hr noshade>
<img src="../gra/button3.png" width=9 height=9 alt="・"> <a href="../exact/exact.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>

サブ gxp.js   Last modified: Mar 25, 2004
サブ xxp.js   Last modified: Mar 25, 2004
サブ io.js   Last modified: Mar 25, 2004

・ 直前のページへ戻る  ・ E-mail to Shigenobu AOKI

Made with Macintosh