メイン ld50.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 MAX_LOOP = 100
var M_SQ_PI2 = 0.398942280401432677940
var LN = 1, LOG = 2, NONE = 0
var useLog = NONE
var title
var lp1 = false
var alpha, beta, ollf, ww, xbar, ssnwx
function probit(p)
{
var a, b, c, d, e, f, g, h, t, u, v, y9, x9
a = 6.9362339826e-13
b = 3.657763036e-11
c = -3.231081277e-09
d = 0.00000008360937017
e = -0.00000104527497
f = 0.000005824238515
g = 0.000006841218299
h = -0.0002250947176
t = -0.0008364353589
u = 0.03706987906
v = 1.570796288
y9 = -Math.log(4 * p * (1 - p))
x9 = Math.sqrt((((((((((((y9 * a + b) * y9 + c) * y9 + d) * y9 + e) * y9 + f) * y9 + g) * y9 + h) * y9 + t) * y9 + u) * y9 + v) * y9))
if (p < 0.5) {
x9 = -x9
}
return x9 + 5
}
function estimation1(x, n, r, num)
{
var wx2, wy, wx, wyx, ww, w, p9, y, temp, i
wx2 = wy = wx = wyx = ww = ollf = 0.0
w = 1.0
for (i = 0; i < num; i++) {
p9 = r[i] / n[i]
if (p9 != 0 && p9 != 1) {
y = probit(p9)
wx2 += w * x[i] * x[i]
wy += w * y
wx += w * x[i]
wyx += w * y * x[i]
ww++
ollf += r[i] * Math.log(p9) + (n[i] - r[i]) * Math.log(1 - p9)
}
}
temp = ww * wx2 - wx * wx
alpha = (wx2 * wy - wx * wyx) / temp
beta = (ww * wyx - wx * wy) / temp
printf("***** 第一近似解 %s\n\n", title)
printf(" α = %g β = %g\n", alpha, beta)
}
function estimation2(x, n, r, num, eps)
{
var ybar, ssnwxy, beta1, alpha1, xi, x9, p9, y, w, wx2, wy, wx, wyx, cllf, z, loop, i
printf("\n***** 反復法によるパラメター推定\n\n")
if (lp1 == true) {
printf("%12s%17s%11s%15s\n", "反復回数", "対数尤度", "α", "β")
}
for (loop = 1; loop <= MAX_LOOP; loop++) {
wx2 = wx = wy = wyx = ww = cllf = 0
for (i = 0; i < num; i++) {
xi = x[i]
y = alpha + beta * xi
x9 = y - 5
p9 = 1 - gxp(x9)
z = Math.exp(-x9 * x9 / 2) * M_SQ_PI2
y += (r[i] / n[i] - p9) / z
w = z * z / (1 - p9) / p9 * n[i]
wx2 += w * xi * xi
wy += w * y
wx += w * xi
wyx += w * y * xi
ww += w
cllf += r[i] * Math.log(p9) + (n[i] - r[i]) * Math.log(1 - p9);
}
xbar = wx / ww
ybar = wy / ww
ssnwxy = wyx - wx * wy / ww
ssnwx = wx2 - wx * wx / ww
beta1 = ssnwxy / ssnwx
alpha1 = ybar - beta1 * xbar
if (lp1 == true) {
printf("%12i%17.7g%15.7g%15.7g\n", loop, cllf, alpha1, beta1)
}
if (Math.abs((alpha - alpha1) / alpha) < eps && Math.abs((beta - beta1) / beta) < eps) {
alpha = alpha1
beta = beta1
if (lp1 == true) {
printf("\n***** 最終解 %s\n\n", title)
}
printf(" α = %g β = %g\n", alpha, beta)
return
}
alpha = alpha1
beta = beta1
}
printf("%i 回反復しましたが，収束しませんでした\n", MAX_LOOP)
}
function goodnessOfFitness(x, n, r, num)
{
var chi, cllf, p9, pobs, q, i
printf("\n***** 直線性の検定（適合度検定）\n\n")
printf("%12s%11s%11s%11s%10s%10s\n", "用量", "個体数(n)", "反応数(r)", "期待値(e)", "p = r/n", "p = e/n")
chi = cllf = 0;
for (i = 0; i < num; i++) {
p9 = 1 - gxp(alpha + beta * x[i] - 5)
pobs = r[i] / n[i]
cllf = cllf + r[i] * Math.log(p9) + (n[i] - r[i]) * Math.log(1 - p9);
printf("%12.6g%8i%11i%14.4f%10.4f%10.4f\n", x[i], n[i], r[i], p9 * n[i], pobs, p9)
q = 1 - p9
if (q == 0) {
q = 0.000000000001
}
else if (q == 1) {
q = 1 - 0.000000000001
}
chi += n[i] * (pobs - p9) * (pobs - p9) / (q - q * q)
}
p9 = xxp(chi, num - 2)
printf("\n")
printf(" カイ二乗値 = %.5f 自由度 = %i\n", chi, num-2)
printf(" Ｐ値 = %.5f\n", p9)
printf(" 実測対数尤度 = %.6g\n", ollf)
printf(" 対数尤度 = %.6g\n", cllf)
}
function outResult(useLog, base)
{
var ed50, se, g, cl, pm, clu0, cll0
ed50 = (5 - alpha) / beta
se = Math.sqrt(1 / ww + Math.pow(ed50 - xbar, 2) / ssnwx) / Math.abs(beta)
g = Math.pow(1.96 / beta, 2) / ssnwx
cl = ed50 + g * (ed50 - xbar) / (1 - g);
pm = 1.96 / beta / (1 - g) * Math.sqrt((1 - g) / ww + Math.pow(ed50 - xbar, 2) / ssnwx)
if (useLog == NONE) {
clu0 = cl + pm
cll0 = cl - pm
}
else {
clu0 = Math.pow(base, cl + pm)
cll0 = Math.pow(base, cl - pm)
}
printf(" LD50 = %g\n", (useLog == NONE) ? ed50 : Math.pow(base, ed50))
printf(" S.E. = %g\n", (useLog == NONE) ? se : Math.pow(base, ed50 * (Math.exp(se) - 1)))
printf(" 95%% 信頼限界 %g < LD50 < %g\n", cll0, clu0)
}
function sub(num, x, n, r, eps)
{
var base, i
printf("★ プロビット法による ＥＤ50，ＬＤ50 などの計算 ★\n\n")
if (useLog == LN) {
printf("***** 自然対数変換を行った\n")
base = Math.exp(1.0)
title = "probit = α + ln(dose)・β"
}
else if (useLog == LOG) {
printf("***** 常用対数変換を行った\n")
base = 10.0
title = "probit = α + log(dose)・β"
}
else {
printf("***** 対数変換せず\n")
title = "probit = α + dose・β"
}
printf(" 収束判定値 = %g\n\n", eps)
if (useLog != NONE) {
for (i = 0; i < num; i++) {
if (x[i] <= 0) {
printf("用量が 0 のものがあるので，対数変換はできません\n")
sep(50)
return
}
x[i] = Math.log(x[i]) / Math.log(base)
}
}
estimation1(x, n, r, num)
estimation2(x, n, r, num, eps)
outResult(useLog, base)
goodnessOfFitness(x, n, r, num)
}
function calc(data_string, eps0)
{
var eps, err, data, i, num
var x = new Array()
var n = new Array()
var r = new Array()
eps = parseFloat(eps0)
if ((data = getdata(data_string, 3)) != false) {
num = data.length
err = false
for (i = 0; i < num; i++) {
x[i] = data[i][0]
n[i] = data[i][1]
r[i] = data[i][2]
if (n[i] != Math.floor(n[i]) || r[i] != Math.floor(r[i])) {
printf("動物数，反応数は整数値でなくてはなりません。\n")
err = true
break
}
else if (n[i] < 0 || r[i] < 0 || r[i] > n[i]) {
printf("動物数 = %i，反応数 = %i というのは変です。\n", n[i], r[i])
err = true
break
}
}
if (err == false) {
sub(num, x, n, r, eps)
}
}
sep(50)
}
//-->
</script>
</head>
<body bgcolor="#ffffff">
<font size="+2"><b>プロビット法による LD<sub>50</sub>，ED<sub>50</sub> などの計算</b></font> <a href="src/ld50.html"><img src="png/src.png" width=35 height=11 alt="src" align=top></a><br> 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, this.form.eps0.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><img src="../gra/button3.png" width=9 height=9 alt="・"> <a href="exa/ld50.html">使用法</a></td>
</tr></table>
<br>
<table>
<tr><td colspan=2><hr width="60%" align=left>用量を対数変換しますか
<input type="radio" name="uselog" checked onClick="useLog=NONE">しない
<input type="radio" name="uselog" onClick="useLog=LN">自然対数で変換
<input type="radio" name="uselog" onClick="useLog=LOG">常用対数で変換<hr width="60%" align=left></td></tr>
<tr><td colspan=2>収束判定条件 <input value="0.000001" name="eps0"><hr width="60%" align=left></td></tr>
<tr><td colspan=2>収束計算過程を出力しますか
<input type="radio" name="lp" onClick="lp1=true">はい
<input type="radio" name="lp" checked onClick="lp1=false">いいえ<hr width="60%" align=left></td></tr>
<tr>
<td colspan=2>入力欄には，用量，動物数，反応数を１行に一対ずつの値を<a href="exa/kugirimoji.html">区切り文字</a>で区切って入力</td>
</tr>
<tr>
<td>入力欄<br><textarea name="data" ROWS=25 COLS=15></textarea></td>
<td>出力欄<br><textarea name="result" ROWS=25 COLS=80></textarea></td>
</tr>
</table>
</form>
<p><hr noshade>
<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
