Path: ccsf.homeunix.org!ccsf.homeunix.org!news1.wakwak.com!nf1.xephion.ne.jp!onion.ish.org!news.daionet.gr.jp!news.yamada.gr.jp!news.media.kyoto-u.ac.jp!not-for-mail
From: "Funatsu Kunihiro" <tmpfuna@netscape.net>
Newsgroups: fj.comp.applications.excel
Subject: Re: =?iso-2022-jp?B?RVhDRUwgGyRCJEclWSE8JT9KLElbJHI3Vzs7JDkka0p9GyhC?= =?iso-2022-jp?B?GyRCSyEbKEI=?=
Date: Mon, 23 Aug 2004 10:41:23 +0900
Organization: Public NNTP Service, Kyoto University, JAPAN
Lines: 98
Sender: tmpfuna@netscape.net
Message-ID: <opsc5wa9ig263y00@pecd218>
References: <040822211825.M0101342@myhost.SoftHome.net>
NNTP-Posting-Host: 221.251.57.237
Mime-Version: 1.0
Content-Type: text/plain; format=flowed; delsp=yes; charset=iso-2022-jp
Content-Transfer-Encoding: 8bit
X-Trace: caraway.media.kyoto-u.ac.jp 1093225286 29289 221.251.57.237 (23 Aug 2004 01:41:26 GMT)
X-Complaints-To: news@news.media.kyoto-u.ac.jp
NNTP-Posting-Date: Mon, 23 Aug 2004 01:41:26 +0000 (UTC)
User-Agent: Opera M2/7.53 (Win32, build 3864)
Xref: ccsf.homeunix.org fj.comp.applications.excel:184

On Sun, 22 Aug 2004 12:18:25 GMT, Taku <tuc@SoftHome.net> wrote:

> 正規分布は組み込み関数としてあるので，簡単に計算できるのですが，
> 今回，ベータ分布を計算する必要が出てきました。

ベータ分布も組み込み関数としてあるでしょ。
どうしても自分で作りたいって言うのなら
Function betai(x As Double, a As Double, b As Double) As Double
のような感じで...

参考: Numerical Recipes　(っていうかそのまま)
Function gammln(xx As Double) As Double
Dim x As Double, y As Double, tmp As Double, ser As Double
Dim cof(5) As Double
cof(0) = 76.1800917294715
cof(1) = -86.5053203294168
cof(2) = 24.0140982408309
cof(3) = -1.23173957245015
cof(4) = 1.20865097386618E-03
cof(5) = -5.395239384953E-06
Dim j As Long
x = xx
y = x
tmp = x + 5.5
tmp = tmp - (x + 0.5) * Log(tmp)
ser = 1.00000000019001
For j = 0 To 5
     y = y + 1#
     ser = ser + cof(j) / y
Next j
gammln = -tmp + Log(2.506628274631 * ser / x)
End Function

Function betacf(a As Double, b As Double, x As Double) As Double
Dim m As Long, m2 As Long
Dim aa As Double, c As Double, d As Double, del As Double, h As Double
Dim qab As Double, qam As Double, qap As Double
Const MAXIT As Long = 100
'Const EPS As Double = 0.0000003
Const EPS As Double = 0.00000000000001
Const FPMIN As Double = 1E-30
qab = a + b
qap = a + 1#
qam = a - 1#
c = 1#
d = 1# - qab * x / qap
If Abs(d) < FPMIN Then d = FPMIN
d = 1# / d
h = d
For m = 1 To MAXIT
     m2 = 2 * m
     aa = m * (b - m) * x / ((qam + m2) * (a + m2))
     d = 1# + aa * d
     If Abs(d) < FPMIN Then d = FPMIN
     c = 1# + aa / c
     If Abs(c) < FPMIN Then c = FPMIN
     d = 1# / d
     h = h * d * c
     aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
     d = 1# + aa * d
     If Abs(d) < FPMIN Then d = FPMIN
     c = 1# + aa / c
     If Abs(c) < FPMIN Then c = FPMIN
     d = 1# / d
     del = d * c
     h = h * del
     If Abs(del - 1#) < EPS Then Exit For
Next m
If m > MAXIT Then
     MsgBox "a or b too big, or MAXIT too small in betacf"
     Exit Function
End If
betacf = h
End Function

Function betai(x As Double, a As Double, b As Double) As Double
Dim bt As Double
If x < 0# Or x > 1# Then
     MsgBox "Bad x in routine betai"
     Exit Function
End If
If x = 0# Or x = 1# Then
     bt = 0#
Else
     bt = Exp(gammln(a + b) - gammln(a) - gammln(b) + a * Log(x) + b *  
Log(1# - x))
End If
If x < (a + 1#) / (a + b + 2#) Then
     betai = bt * betacf(a, b, x) / a
Else
     betai = 1# - bt * betacf(b, a, 1# - x) / b
End If
End Function

-- 
Funatsu Kunihiro
nospamtmpfuna@netscape.net
Using M2, Opera's revolutionary e-mail client: http://www.opera.com/m2/