メイン twowayANOVA.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="fxp.js">document.write("fxp.js ファイルが見つかりません??<br>")</script>
<script src="inv.js">document.write("inv.js ファイルが見つかりません??<br>")</script>
<script src="io.js">document.write("io.js ファイルが見つかりません??<br>")</script>
<script language="JavaScript">
<!--
var epsinv = 1e-8
function minf(x, n)
{
var i, retv
retv = x[0]
for (i = 1; i < n; i++) {
if (x[i] < retv) {
retv = x[i]
}
}
return retv
}
function maxf(x, n)
{
var i, retv
retv = x[0]
for (i = 1; i < n; i++) {
if (x[i] > retv) {
retv = x[i]
}
}
return retv
}
function mreg(data, nc, nv, result)
{
var mean, variance, r, sd, ss, i, j, k, nv1, se, syy, temp2
mean = new Array(nv)
r = makeMatrix(nv, nv)
sd = new Array(nv)
ss = new Array(nv)
for (i = 0; i < nv; i++) {
mean[i] = 0
for (j = 0; j < nv; j++) {
r[i][j] = 0
}
}
for (i = 0; i < nc; i++) {
for (j = 0; j < nv; j++) {
mean[j] += data[i][j]
}
}
for (j = 0; j < nv; j++) {
mean[j] /= nc
}
for (i = 0; i < nc; i++) {
for (j = 0; j < nv; j++) {
for (k = 0; k < nv; k++) {
r[j][k] += (mean[j]-data[i][j])*(mean[k]-data[i][k])
}
}
}
for (j = 0; j < nv; j++) {
ss[j] = r[nv-1][j]
variance = r[j][j]/(nc-1)
sd[j] = Math.sqrt(variance)
if (variance <= 0) {
printf("分散が０になる変数があります\n")
return
}
}
for (j = 0; j < nv; j++) {
for (k = 0; k < nv; k++) {
r[j][k] /= (nc-1)*sd[j]*sd[k]
}
}
nv1 = nv-1
if (inv(r, nv1, epsinv) == 0) {
se = syy = ss[nv1]
for (i = 0; i < nv1; i++) {
temp2 = 0
for (j = 0; j < nv1; j++) {
temp2 += r[i][j]*r[j][nv1]
}
se -= ss[i]*temp2*sd[nv1]/sd[i]
}
result[0] = syy-se
result[1] = se
result[2] = syy
}
}
function type4(x, is, js, nc, max_is, max_js)
{
var data, result, nvi, nvj, nvij, nv, i, j, ssabab, sse, sst, ssacb, ssab, ssadb, ssbda, dfadb, dfbda, dfab, dfe, dft, msadb, msbda, msab, mse, mst, fa, fb, fab, pa, pb, pab
result = new Array(3)
// full model
nvi = max_is-1
nvj = max_js-1
nvij = nvi* nvj
nv = nvi+nvj+nvij+1
data = makeMatrix(nc, nv)
for (i = 0; i < nc; i++) {
data[i][nv-1] = x[i]
for (j = 0; j < nv-1; j++) {
data[i][j] = 0
}
if (is[i] > 1) {
data[i][is[i]-2] = 1
}
if (js[i] > 1) {
data[i][js[i]-2+nvi] = 1
}
if (is[i] > 1 && js[i] > 1) {
data[i][(is[i]-1)*nvj+js[i]-2+nvi] = 1
}
}
mreg(data, nc, nv, result)
ssabab = result[0]
sse = result[1]
sst = result[2]
// factor A, B
nvi = max_is-1
nvj = max_js-1
nv = nvi+nvj+1
// data = makeMatrix(nc, nv)
for (i = 0; i < nc; i++) {
data[i][nv-1] = x[i]
for (j = 0; j < nv-1; j++) {
data[i][j] = 0
}
if (is[i] > 1) {
data[i][is[i]-2] = 1
}
if (js[i] > 1) {
data[i][js[i]-2+nvi] = 1
}
}
mreg(data, nc, nv, result)
ssacb = result[0]
ssab = ssabab-ssacb
// factor B
nv = max_is
for (i = 0; i < nc; i++) {
data[i][nv-1] = x[i]
for (j = 0; j < nv-1; j++) {
data[i][j] = 0
}
if (is[i] > 1) {
data[i][is[i]-2] = 1
}
}
mreg(data, nc, nv, result)
ssbda = ssacb-result[0]
// factor A
nv = max_js
for (i = 0; i < nc; i++) {
data[i][nv-1] = x[i]
for (j = 0; j < nv-1; j++) {
data[i][j] = 0
}
if (js[i] > 1) {
data[i][js[i]-2] = 1
}
}
mreg(data, nc, nv, result)
ssadb = ssacb-result[0]
dfadb = max_is-1
dfbda = max_js-1
dfab = dfadb*dfbda
dfe = nc-max_is*max_js
dft = nc-1
msadb = ssadb/dfadb
msbda = ssbda/dfbda
msab = ssab/dfab
mse = sse/dfe
mst = sst/dft
fa = msadb/mse
fb = msbda/mse
fab = msab/mse
pa = fxp(fa, dfadb, dfe)
pb = fxp(fb, dfbda, dfe)
pab = fxp(fab, dfab, dfe)
printf("%8s %15s %8s %15s %8s %8s\n", "変動要因", "平方和", "自由度", "平均平方", "Ｆ値", "Ｐ値")
printf("%8s %15.7g %8i %15.7g %8.3f %8.3f\n", "要因 A", ssadb, dfadb, msadb, fa, pa)
printf("%8s %15.7g %8i %15.7g %8.3f %8.3f\n", "要因 B", ssbda, dfbda, msbda, fb, pb)
printf("%8s %15.7g %8i %15.7g %8.3f %8.3f\n", "交互作用", ssab, dfab, msab, fab, pab)
printf("%8s %15.7g %8i %15.7g\n", "残差", sse, dfe, mse)
printf("%8s %15.7g %8i %15.7g\n", "全体", sst, dft, mst)
}
function get_gm(x, nc)
{
var gm, i
gm = 0
for (i = 0; i < nc; i++) {
gm += x[i]
}
return gm/nc
}
function get_gs(x, nc, gm)
{
var gs, i
gs = 0
for (i = 0; i < nc; i++) {
gs += Math.pow(x[i]-gm, 2)
}
return gs
}
function get_m(x, nc, level, level_a)
{
var i, n, m
n = 0
m = 0
for (i = 0; i < nc; i++) {
if (level_a[i] == level) {
n++
m += x[i]
}
}
return m/n
}
function get_s2(x, nc, level_i, level_a, level_j, level_b)
{
var i, n, m, s
n = 0
m = 0
s = 0
for (i = 0; i < nc; i++) {
if (level_a[i] == level_i && level_b[i] == level_j) {
n++
m += x[i]
}
}
m /= n
for (i = 0; i < nc; i++) {
if (level_a[i] == level_i && level_b[i] == level_j) {
s += Math.pow(x[i]-m, 2)
}
}
return s
}
function type1(x, is, js, nc, max_is, max_js)
{
var i, j, gm, sst, ssa, ssb, sse, dfa, dfb, dfe, msa, msb, mse, si, sj, fa, fb,
gm = get_gm(x, nc)
sst = get_gs(x, nc, gm)
si = 0
for (i = 0; i < max_is; i++) {
si += Math.pow(get_m(x, nc, i+1, is)-gm, 2)
}
sj = 0
for (j = 0; j < max_js; j++) {
sj += Math.pow(get_m(x, nc, j+1, js)-gm, 2)
}
ssa = si*max_js; dfa = max_is-1; msa = ssa/dfa
ssb = sj*max_is; dfb = max_js-1; msb = ssb/dfb
sse = sst-ssa-ssb; dfe = dfa*dfb; mse = sse/dfe
fa = msa/mse
fb = msb/mse
printf("%8s %15s %8s %15s %8s %8s\n", "変動要因", "平方和", "自由度", "平均平方", "Ｆ値", "Ｐ値")
printf("%8s %15.7g %8i %15.7g %8.3f %8.3f\n", "要因Ａ", ssa, dfa, msa, fa, fxp(fa, dfa, dfe))
printf("%8s %15.7g %8i %15.7g %8.3f %8.3f\n", "要因Ｂ", ssb, dfb, msb, fb, fxp(fb, dfb, dfe))
printf("%8s %15.7g %8i %15.7g\n", "残差", sse, dfe, mse)
printf("%8s %15.7g %8i %15.7g\n", "全体", sst, nc-1, sst/(nc-1))
}
function type3(x, is, js, nc, max_is, max_js, cases)
{
var i, j, gm, sst, ssa, ssb, sse, dfa, dfb, dfe, msa, msb, mse, si, sj, tn, ssab, dfab, msab, dfe, fa, fb, fab
gm = get_gm(x, nc)
sst = get_gs(x, nc, gm)
si = 0
for (i = 0; i < max_is; i++) {
tn = 0
for (j = 0; j < max_js; j++) {
tn += cases[i][j]
}
si += tn*Math.pow(get_m(x, nc, i+1, is)-gm, 2)
}
sj = 0
for (j = 0; j < max_js; j++) {
tn = 0
for (i = 0; i < max_is; i++) {
tn += cases[i][j]
}
sj += tn*Math.pow(get_m(x, nc, j+1, js)-gm, 2)
}
sse = 0
for (i = 0; i < max_is; i++) {
for (j = 0; j < max_js; j++) {
sse += get_s2(x, nc, i+1, is, j+1, js)
}
}
ssa = si; dfa = max_is-1; msa = ssa/dfa
ssb = sj; dfb = max_js-1; msb = ssb/dfb
ssab = sst-ssa-ssb-sse; dfab = dfa*dfb; msab = ssab/dfab
dfe = nc-max_is*max_js; mse = sse/dfe
fa = msa/mse
fb = msb/mse
fab = msab/mse
printf("\nモデル I（母数モデル）\n")
printf("%8s %15s %8s %15s %8s %8s\n", "変動要因", "平方和", "自由度", "平均平方", "Ｆ値", "Ｐ値")
printf("%8s %15.7g %8i %15.7g %8.3f %8.3f\n", "要因Ａ", ssa, dfa, msa, fa, fxp(fa, dfa, dfe))
printf("%8s %15.7g %8i %15.7g %8.3f %8.3f\n", "要因Ｂ", ssb, dfb, msb, fb, fxp(fb, dfb, dfe))
printf("%8s %15.7g %8i %15.7g %8.3f %8.3f\n", "交互作用", ssab, dfab, msab, fab, fxp(fab, dfab, dfe))
printf("%8s %15.7g %8i %15.7g\n", "残差", sse, dfe, mse)
printf("%8s %15.7g %8i %15.7g\n", "全体", sst, nc-1, sst/(nc-1))
fa = msa/msab
fb = msb/msab
fab = msab/mse
printf("\nモデル II（変量モデル）\n")
printf("%8s %15s %8s %15s %8s %8s\n", "変動要因", "平方和", "自由度", "平均平方", "Ｆ値", "Ｐ値")
printf("%8s %15.7g %8i %15.7g %8.3f %8.3f\n", "要因Ａ", ssa, dfa, msa, fa, fxp(fa, dfa, dfab))
printf("%8s %15.7g %8i %15.7g %8.3f %8.3f\n", "要因Ｂ", ssb, dfb, msb, fb, fxp(fb, dfb, dfab))
printf("%8s %15.7g %8i %15.7g %8.3f %8.3f\n", "交互作用", ssab, dfab, msab, fab, fxp(fab, dfab, dfe))
printf("%8s %15.7g %8i %15.7g\n", "残差", sse, dfe, mse)
printf("%8s %15.7g %8i %15.7g\n", "全体", sst, nc-1, sst/(nc-1))
}
function calc(data_string)
{
var data, is, js, x, nc, i, j, min_is, min_js, max_is, max_js, cases, same, hirei, ratio
if ((data = getdata(data_string, -1)) != false) {
printf("★ 二元配置分散分析\n\n")
nc = 0
is = new Array()
js = new Array()
x = new Array()
for (i = 0; i < data.length; i++) {
if (data[i].length < 3) {
printf("データの記述に誤りがあります\n")
return
}
for (j = 2; j < data[i].length; j++) {
is[nc] = data[i][0]
js[nc] = data[i][1]
x[nc++] = data[i][j]
}
}
min_is = minf(is, nc)
min_js = minf(js, nc)
if (min_is != 1 || min_js != 1) {
printf("水準を表す数値は1から始めてください %i %i\n", min_is, min_js)
return
}
max_is = maxf(is, nc)
max_js = maxf(js, nc)
cases = makeMatrix(max_is, max_js)
for (i = 0; i < max_is; i++) {
for (j = 0; j < max_js; j++) {
cases[i][j] = 0
}
}
for (i = 0; i < nc; i++) {
cases[is[i]-1][js[i]-1]++
}
same = true
for (i = 0; i < max_is; i++) {
for (j = 0; j < max_js; j++) {
if (cases[i][j] == 0) {
printf("観察データがない，水準の組合せがあります\n")
return
}
else if (cases[0][0] != cases[i][j]) {
same = false
}
}
}
if (same == true) {
if (cases[0][0] == 1) {
type1(x, is, js, nc, max_is, max_js)
}
else {
type3(x, is, js, nc, max_is, max_js, cases)
}
}
else {
hirei = true
for (i = 0; i < max_is-1; i++) {
ratio = cases[i+1][0]/cases[i][0]
for (j = 1; j < max_js; j++) {
if (cases[i+1][j]/cases[i][j] != ratio) {
hirei = false
break
}
}
if (hirei == false) {
break
}
}
if (hirei == true) {
type3(x, is, js, nc, max_is, max_js, cases)
}
else {
type4(x, is, js, nc, max_is, max_js)
}
}
sep(40)
}
}
//-->
</script>
</head>
<body bgcolor="#ffffff">
<font size="+2"><b>二元配置分散分析</b></font> <a href="src/twowayANOVA.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/twowayANOVA.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="../lecture/TwoWayANOVA/TwoWay.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>
サブ fxp.js Last modified: Mar 25, 2004
サブ inv.js Last modified: Mar 25, 2004
サブ io.js Last modified: Mar 25, 2004
直前のページへ戻る E-mail to Shigenobu AOKI
