|
Post by bitvast on May 27, 2014 14:21:36 GMT 1
I'm working on a binomial distribution calculator and have made a start writing some code to get the binomial coefficients, which you can get from the nCr function on a scientific calculator. The code works, but only for small values of n. FUNCTION binom(NUMBER n, NUMBER k) LOCAL coeff = 1 FOR i = n-k+1 TO n coeff = coeff * i NEXT i FOR i = 1 TO k coeff = coeff / i NEXT i RETURN coeff END FUNCTION
e.g. binom(17, 8) correctly returns 24,310. But when I use larger values, the returned value is incorrect. For example, binom(40,8) returns -4203 (it should be 76,904,685). I think this must be due to an overflow, but NUMBER is equivalent to long int which should be big enough? I've tried using double but still don't get the correct answers. Any advice?
|
|
|
Post by bigbass on May 27, 2014 15:21:50 GMT 1
Hey bitvast
used some C but it worked maybe there is still a better way to do it
FUNCTION binom(long long n, long long k) LOCAL coeff TYPE long long coeff = 1 FOR i = n-k+1 TO n coeff = coeff * i NEXT i FOR i = 1 TO k coeff = coeff / i NEXT i RETURN coeff END FUNCTION
PRINT binom(40,8)
|
|
|
Post by bitvast on May 27, 2014 16:25:01 GMT 1
Thanks Joe.
I still get an error for larger values e.g. binom(100,50), but my HP 35s calculator gives 1.0089 x 1029. There must be a smarter algorithm.
Here's the C version on Rosetta Code (evaluate binomial coefficients) -
#include <stdio.h> #include <limits.h> typedef unsigned long ULONG; ULONG binomial(ULONG n, ULONG k) { ULONG r = 1, d = n - k; /* choose the smaller of k and n - k */ if (d > k) { k = d; d = n - k; } while (n > k) { if (r >= UINT_MAX / n) return 0; /* overflown */ r *= n--; /* divide (n - k)! as soon as we can to delay overflows */ while (d > 1 && !(r % d)) r /= d--; } return r; } int main() { printf("%lu\n", binomial(5, 3)); return 0; }
But I still get 0 for binomial(100, 50).
|
|
|
Post by bigbass on May 27, 2014 16:33:14 GMT 1
Hey bitvast
Will look at the rosetta stuff
here is a way to get commas in the large numbers but a little ugly C to make it work
explained: To get the commas you need a header locale.h then you need these functions from C setlocale and printf so we PROTO them this allows us to use them now in bacon
PRAGMA INCLUDE locale.h
PROTO setlocale PROTO printf
FUNCTION binom(long long n, long long k) LOCAL coeff TYPE long long coeff = 1 FOR i = n-k+1 TO n coeff = coeff * i NEXT i FOR i = 1 TO k coeff = coeff / i NEXT i RETURN coeff END FUNCTION
PRINT binom(40,8)
'--- get some commas in there to make it readable setlocale(LC_NUMERIC, ""); printf("%'d\n", binom(40,8))
|
|
|
Post by bitvast on May 27, 2014 17:04:41 GMT 1
Thanks again, Joe.
The javacript version works fine with big numbers and it uses the same algorithm as the BaCon code:
<script> function binom(n, k) { var coeff = 1; for (var i = n-k+1; i <= n; i++) coeff *= i; for (var i = 1; i <= k; i++) coeff /= i; return coeff; } alert(binom(100,50)); </script>
1.0089134454556415e+29
|
|
|
Post by bitvast on May 27, 2014 17:12:40 GMT 1
Problem solved. I should have used double with FORMAT "%e".
FUNCTION binom(double n, double k) LOCAL coeff = 1 TYPE double FOR i = 1 TO k coeff = coeff * (n - i + 1) / i NEXT i RETURN coeff END FUNCTION
PRINT binom(100, 50) FORMAT "%e"
1.008913e+29
|
|
|
Post by bigbass on May 27, 2014 17:27:20 GMT 1
Hey bitvast
Excellent much better with e notation Glad to see you solved it
(The commas may come in handy for another project though)
Joe
|
|
|
Post by bitvast on May 27, 2014 17:39:16 GMT 1
(The commas may come in handy for another project though) Absolutely. So no C was necessary since FLOATING = double. Which is how it should be, I suppose.
|
|
|
Post by Pjot on May 27, 2014 17:42:49 GMT 1
Hi bitvast, On my 64-bit system your program gave the correct answer right away, because in 64 bit, a 'long' consists of 8 bytes. So I guess you use a 32 bit system, where a long consists of 4 bytes. But a 'double' is always 8 bytes, both on 32bit as on 64 bit. Some compilers accept the 'long long' type, this also will be 8 bytes on a 32bit system. So you can try this as well. And if you need to work with really long numbers, you can use the GMP interface also. BR Peter
|
|
|
Post by vovchik on May 27, 2014 18:03:33 GMT 1
Dear guys,
This will give you a pretty long one without the exponential notation:
PRINT binom(100, 50) FORMAT "%.0lf\n"
With kind regards, vovchik
|
|
|
Post by bitvast on May 27, 2014 18:25:33 GMT 1
|
|
|
Post by bitvast on May 28, 2014 12:12:04 GMT 1
type errors again. I added code to calculate the probability and this should output 0.24609375 FUNCTION binom(FLOATING n, FLOATING r, FLOATING p) LOCAL coeff = 1 TYPE FLOATING FOR i = 1 TO r coeff = coeff * (n - i + 1) / i NEXT i RETURN coeff * POW(p, r) * POW(1-p, n-r) END FUNCTION
PRINT binom(9, 5, 0.5) FORMAT "%.2f" But I get this error And the output is 5.00.
|
|
|
Post by bitvast on May 28, 2014 13:35:09 GMT 1
It definitely seems to be a FUNCTION related bug. Same code without the FUNCTION:
DECLARE p = 0.5, n = 9, r = 5, coeff = 1 TYPE FLOATING
FOR i = 1 TO r coeff = coeff * (n - i + 1) / i NEXT i
PRINT coeff * POW(p, r) * POW(1-p, n-r) FORMAT "%.8f"
Output is correct:
|
|
|
Post by bigbass on May 28, 2014 14:54:35 GMT 1
Hey bitvast only a very minor adjustment and it works
BTW having baconize math functions is a great thing all this type of work you are doing is very useful maybe later (some tutorial) you could show the standard math formula and what it looks like later in bacon
Joe
FUNCTION binom(FLOATING n, FLOATING r, FLOATING p) LOCAL coeff = 1 ,result TYPE FLOATING FOR i = 1 TO r coeff = coeff * (n - i + 1) / i NEXT i result = coeff * POW(p, r) * POW(1-p, n-r) RETURN result END FUNCTION
PRINT binom(9, 5, 0.5) FORMAT "%.8f"
0.24609375
|
|
|
Post by bitvast on May 28, 2014 16:20:49 GMT 1
Thanks, Joe.
I guess my lack of C knowledge is letting me down. I didn't know you had to use an intermediate variable.
Yep, working on a little statistics library which I'll post here with a tutorial when I'm done.
|
|