CODEWITHSUNDEEP
netlogo
Ignacio
Ignacio

Reputation: 69

NetLogo computing beta distribution/function

I am using the stats extension in NetLogo to compute a beta function:

set z (stats:beta (H + 1) (T + 1))

Where H and T are numbers of heads and tails in a coin flip. The use of the stats extension is not essential, I have the same problem when using the factorial expression of the beta function.

The problem is that when H + T > 168, NetLogo reports z = 0 and there are some procedures I cannot perform [in particular, the beta distribution]

Is there any way to approximate the beta function (or distribution) in Netlogo such that it does not run into this problem?

Upvotes: 4

Views: 530

Answers (3)

Ignacio
Ignacio

Reputation: 69

I finally managed to solve my own question. I hope this is helpful for other people too. The main problem is to compute large factorials - as it came up in the comments. In particular, 171! is simply too largo to compute. So the optimal solution would be to minimize the need of computing such numbers.

For the discrete case, the beta function is the following:

B(alpha,beta) = (alpha - 1)!(beta - 1)! / (alpha + beta - 1)!

The mistake was to try to compute first nominator, then denominator, and then divide; because denominator gets too large. I found two solutions, one of them is better than the other but much more complicated - but actually I think it is close to optimal.

The less optimal but shorter is based on the idea of simplifying the denominator using the largest nominator:

to-report betafunc [alpha beta]
  ifelse (alpha >= beta)
  [let v1 (factorial (beta - 1) (0))
   let v2 (factorial (alpha + beta - 1) (alpha - 1))
   report (v1 / v2)]
  [let v1 (factorial (alpha - 1) (0))
   let v2 (factorial (alpha + beta - 1) (beta - 1))
  report (v1 / v2)]
end

where 

to-report factorial [n m]
 if (n = m) [report 1]
 report (n * factorial (n - 1) (m))
end

The problem with this one is that denominator still grows very fast. The idea for the second solution is to proceed by simplifying each step of the factorial.

For example, if we have (5!8!)/(12)! we can express it as the product of (5.8/12).(4.7/11).(3.6/10) ... (1.4/8).(1.3/7)...(1.1/5)(1.1/4)...(1.1/1)

By stepwise simplifying each term of the sequence we make sure that things do not grow up that fast. Here is my solution:

to-report stepwisefactorial3 [n1 n2 d]
  if (n1 = 0) [report stepwisefactorial2 (n2) (d)]
  if (n2 = 0) [report stepwisefactorial2 (n1) (d)]
  if (d = 0)  [report (stepwisefactorial1 (n1) * stepwisefactorial1 (n2))]
  report ((((n1 * n2) / d)) * stepwisefactorial3 (n1 - 1) (n2 - 1) (d - 1))
end

to-report stepwisefactorial2 [n d]
  if (n = 0) [report (1 / (stepwisefactorial1 (d)))]
  if (d = 0) [report stepwisefactorial1 (n)]
  report ((n / d)* stepwisefactorial2 (n - 1) (d - 1))
end

to-report stepwisefactorial1 [d]
  if d = 0 [ report 1 ]
  report d * stepwisefactorial1 (d - 1)
end

Still things blow up if alpha + beta > 345; but it is progress.

Upvotes: 2

Charles
Charles

Reputation: 4168

Building on Chris' answer, the stats extension does have a log of Gamma function, stats:loggamma. It handles an argument of considerably over 171 and thus a and b can add to well over 172 as well.

observer> show stats:beta 85 86
observer: 1.2864854397253604E-52
observer> show exp (stats:loggamma 85 + stats:loggamma 86 - stats:loggamma (85 + 86))
observer: 1.2864854397251408E-52
observer> show stats:beta 86 86
observer: 0
observer> show exp (stats:loggamma 86 + stats:loggamma 86 - stats:loggamma (86 + 86))
observer: 6.394810665301235E-53
observer> show exp (stats:loggamma 200 + stats:loggamma 200 - stats:loggamma (200 + 200))
observer: 9.713217247613997E-122

A new version of the stats extension (v1.4.0) has been posted with a "bigBeta" function that uses the logs.

Upvotes: 2

Chris
Chris

Reputation: 61

(This is a side comment, but the question was already answered, anyway.)

I wondered about using logs.

Beta(x,y) = Gamma(x) * Gamma(y) / Gamma(x + y)
= exp ( ln(Gamma(x)) + ln(Gamma(y)) - ln(Gamma(x + y)) )

Playing with that a bit shows that the Stats extension is smart enough to calculate Beta this way, and what it really comes down to is the overflow in Gamma.

observer> show stats:gamma 171
observer: 7.257415615308E306
observer> show stats:gamma 172
observer: Infinity

(I cribbed here http://www.stata.com/support/faqs/statistics/calculating-beta-function/ .)

Upvotes: 1

Related Questions