posted January 27, 2007 06:43 PM            

new a=0, i=1, s=1
{
a+(s/i); i+2; s=-s; if i LT 1000; then repeat
}
$ pi by MacLaurin series for atan(1) = `a*4`

and this one is far more efficient than the others: http://personal.bgsu.edu/~carother/pi/Pi3a.html
It starts with a square forming four cords in the circle.
dd is (0.5 * length of each chord)^2

new dd=0.5, i=4
{
dd=( power(1-power(1-dd,0.5),2)+dd )/4; i*2; if i LT 65536 then repeat
}
$ pi by Archimedes' Method of Exhaustion: = `power(dd,.5)*i`

Re:>The Computation of Certain Numbers Using a Ruler and Compass
Brad, that looks intriguing, will have to study it some more.

posted January 27, 2007 06:43 PM            

Plugging Stirling's log factorial into the program I put forth
earlier enables one to increase N into the trillions & still
get one's result in a flash. The bad news is that the Stirling
formula requires that you know Pi (or an appromixation to it,
such as PB's ATN(1)*4) -- which sounds like cheating, or at least
circularity.

The next question I have is one I asked a year or or more ago on
this forum: What is the closest approximation to Pi one could hope
to get by random simulations -- on a current PC, using PowerBASIC?

Here's one such widely-used simulation

N=1000 'set as you wish
seed=0 'set as you wish
RANDOMIZE seed
FOR I=1 to N
X=RND(1)
Y=RND(1)
r=SQR(X^2+Y^2)
IF r < 1 then incr Count
NEXT
' "expected" value of Count/N = ATN(1), so...
PRINT " PB Pi = "; ATN(1)*4
PRINT "estimated Pi = "; (Count/N)*4

note added Jan 29:
For a clear explanation the logic of the above program see: http://math.fullerton.edu/mathews/n2003/MonteCarloPiMod.html
or just Google "Monte Carlo pi"
------------------

posted January 28, 2007 01:27 PM            

This is more a function of the data structures you use than about
simulation. If you mean performing an experiment and finding the
ratio of successes (as your example does) then you are limited
by the size of the integer you use. The closest you can get to pi
using PowerBASIC's longs is 245850922/78256779 (actually that may
not be true but I got tired of searching any higher) so this
limits how close you can get to pi. If you use longer integers
you probably would have to write your own random number generator
since it appears PowerBASIC uses longs for this anyway.

posted January 28, 2007 03:02 PM            

I think that precision in pi (from RND's) will also depend on the
precision of the numbers one uses for X, Y & r. I.e., to get
infinitely precise estimate of pi one would need an infinite
number of decimal places.

Maybe there's a way to figure these relationships precisely...
Maybe too we're now in the domain of "Frivolous Arithmetic", as
the Wolfram math website puts it. But thanks anyway!

posted January 28, 2007 10:03 PM            

new a=0, b=.707, c=1e5, n=c, i=c; randomize(timer())
{
a=random(1); a*a; if a+b GE 1 then c--
b=a; if i-- GT 0 then repeat
}
$ Estimate of PI by random numbers `4*c/n`

In terms of CPU time, I think this is comparable with some of the
other deterministic methods.

posted January 28, 2007 10:57 PM            

Probably so, if 4 or 5 decimal places is all you need. But your
deterministic ATN method can get 9 decimal places in a minute or
two on my 4-yr-old Dell, and to get that with the Randoms method
would be a fluke even with unlimited time &/or N.

By the way, with N=10, your ATN method gets the first decimal place
accurate (with very slight rounding error). With N=100 it gets
2 decimal places; with N=1000 it's 3 ... and so on, as far out as
I had patience to try. I think that that function would be the
theoretical "expected" for the mean of many analogous runs (with
different Ns) using the randomization method.

------------------

posted January 29, 2007 03:04 AM            

We should always divide at the latest to avoid multiplying any error

So, replace (Count/N)*4 with 4*Count/N

posted January 29, 2007 06:18 PM            

The code below tests Charles' "McClaurin" routine and a
randomization routine against each other and against PB's
ATN(1) function.

dim N as quad, seed as quad, p as ext
N=1000
p=ATN(1)   'theoretical "expected" probability, for Randy
seed=timer
randomize seed
CLS
PRINT "Estimating pi/4 or ATN(1)"
PRINT "With N =  ";N
PRINT "PB_ATN(1) ";ATN(1)
PRINT "McLaurin  ";McLaurin (N)
PRINT "Randy     ";Randy (N, p)

FUNCTION Randy (N as quad, p as ext) as ext
    dim tot as quad, i as quad, X as ext
    'RND is uniformly distributed beteen 0 and 1
    'so answer tot/N should be close to the "expected" p
    FOR i = 1 to N
        X = RND(1)
        if X <= p then incr tot
    next
    FUNCTION = tot/N
END FUNCTION
'McLaurin series for computing ATN(1)
FUNCTION McLaurin (N as quad) as ext
    dim tot as ext, s as quad, i as quad
    s=1
    i=1
    Do 
        incr tot, (s/i)
        incr i, 2
        s = -s
    loop while i <= N
    FUNCTION = tot
END FUNCTION

------------------

posted January 29, 2007 07:01 PM            

Of course Randy assumes that you already know PI, which is okay for Plato but not
for Darwin.

I wonder if there is a way of discovering e and phi through random processes,
without circularity.

Incidently, I was amazed to read that phi / 2 = sin( 666 ). with 666 in degrees.
Of course that is pentagon geometry. Sin(54).

posted January 29, 2007 07:21 PM            

That's not a documented claim for RND.

posted January 29, 2007 10:40 PM            

G. Lakoff & R. Nunez, in their book "Where does mathematics come from?"
have a whole chapter on e. I wonder what mathematicians have to
say about it -- or them.

Michael,
RND's statistical distribution might not be documented in PB manuals,
but it is easy to test for oneself.

The distribution is, of course, easily modified. E.g., to
generate (pseudo-) random normally distributed numbers with
theoretically expected mean=0 and standard deviation = 1 use:

z = cos(2 * pi * RND) * sqr ( -2 * log (RND))

A search of these PB forums for "random distribution" will give
more details; for example: http://www.powerbasic.com/support/forums/Forum6/HTML/004920.html

------------------

posted January 30, 2007 10:08 AM            

Generating Probability Distribution for p=.5 with 1000000 repeats of 26
Compared with binomial distribution for (.5+.5)^26

This is equivalent to doing 26 million coin flips.

The random distribution seems to concentrate slightly in the middle while
being lower at the extremes, when compared with the binomial distribution.
I have noticed this on several runs, so it may be systemic in the RND() function.

0
1.49011611938476E-8
#0
0
3.8743019104004E-7
#1
.000003
4.84287738800048E-6
#2
.000031
3.8743019104004E-5
#3
.000207
2.22772359848022E-4
#4
.000989
9.801983833313E-4
#5
_.003408
_3.43069434165955E-3
#6
____.009672
____9.80198383331297E-3
#7
_________.023039
_________2.32797116041185E-2
#8
__________________.046529
__________________4.65594232082365E-2
#9
_______________________________.079172
_______________________________7.91510194540031E-2
#10
_____________________________________________.115425
_____________________________________________.115128755569458
#11
_______________________________________________________.14355
________________________________________________________.143910944461822
#12
____________________________________________________________.15509
____________________________________________________________.154981017112731
#13
________________________________________________________.144439
________________________________________________________.143910944461823
#14
____________________________________________.114994
_____________________________________________.115128755569459
#15
_______________________________.079065
_______________________________7.91510194540031E-2
#16
__________________.046779
__________________4.65594232082366E-2
#17
_________.023282
_________2.32797116041185E-2
#18
____.009725
____9.80198383331297E-3
#19
_.003416
_3.43069434165954E-3
#20
.000938
9.801983833313E-4
#21
.000211
2.22772359848023E-4
#22
.000031
3.8743019104004E-5
#23
.000004
4.84287738800049E-6
#24
.000001
3.8743019104004E-7
#25
0
1.49011611938476E-8
#26

script used:

new x=26,g=1e4, p=.5, seed=1
new j=g, i, t, m=x+1, d[m]=0, pp=0, c, sc, mx=0
$ Generating Probability Distribution for p=`p` with `g` repeats of `x`
$ Compared with binomial distribution for (`p`+`1-p`)^`x`
randomize(seed)
{
i=x; t=0
{
if random(1) GE p then t++
if --i GT 0 then repeat
}
d[t+1]++; if --j GT 0 then repeat
}
i=m
{
mx=max(mx,d[i]); if --i GT 0 then repeat
}
i=m; c=1; sc=60/(mx/g)
new n=x,r=0, q=1-p // for binomial display
{
many ("_",d[c++]/g to pp *sc )+num(pp)?
many ("_",exp(lfactorial(n)-lfactorial(r)-lfactorial(n-r) ) _
*power(p,n-r)*power(q,r) to pp *sc)+num(pp)?; r++
$#`m-i`

if --i GT 0 then repeat
}

posted January 30, 2007 11:12 AM            

But your more specific question is how closely does a probability
distribution based on RND fit theoretical "expected" values of the
binomial expansion. The short answer is "not bad at all" by standard
statistical tests of "goodness of fit", but by no means exact
except by fluke regardless of the size of one's sample, so long
as it is finite in size and one's numbers are also finite in size.

The rest is something I had written before seeing your post.
===============================================================
This note might be more appropriate for the recent threads about
"Belief versus proof" or "How to argue" but I leave that problem
to you. If that sounds weird, check out Pascal. I have not as yet
studied him in detail, but maybe it is time for me to do so.
http://plato.stanford.edu/entries/pascal-wager/

--quote from The New York Times, Jan 30, 2007 re http://www.LongBets.org
Long Bets is a nonprofit foundation that calls itself an “arena for
competitive, accountable predictions.” It lets anyone make a
prediction and take wagers on it, with the proceeds going to a
charity named by the winner. The bets made so far are from $200 to
$10,000, on topics ranging from the driving habits of Americans
in 2010 to whether the universe will stop expanding. Mitchell Kapor,
the software guru, is betting that in 2029 no computer will have
passed the Turing test (by conversing so much like a human that you
couldn’t tell the difference). The physicist Freeman Dyson’s money
is on the first extraterrestrial life’s being found somewhere other
than a planet or its satellite.

Five years ago, Dr. Rees posted this prediction: “By 2020, bioterror
or bioerror will lead to one million casualties in a single event.”
He reasoned that “by 2020 there will be thousands — even millions —
of people with the capability to cause a catastrophic biological
disaster. My concern is not only organized terrorist groups, but
individual weirdos with the mindset of the people who now design
computer viruses.”

He didn’t get any takers on LongBets.org, which seems to me a missed
opportunity. So I’ve posted an offer there to bet him $200 — not a
huge sum, but enough to put both our reputations on the line. I
realize that betting on disaster may sound ghoulish, but neither
of us will personally profit (if I win, the money goes to the
International Red Cross). And I think bets like this serve a
purpose.

Besides stimulating public debate, they focus the issue and
discipline prophets. No matter how good their intentions, prophets
face strong temptations to hype...
--end of quote

------------------

posted January 31, 2007 06:29 PM            

Successive Approximation using random numbers

Here is an algorithm that uses randomised feedback to converge on any constant.
The random element makes it fairly robust. The feedback factor, fbr needs to be
adjusted to amplify the feedback. It is better set too high than too low.

I have tried it with sqroot, phi, pi and e.

Here is the algorithm for pi which usually converges to 12 decimal places in
60-100 loops. Not bad when you consider that a binary successive approximation
would take 40 loops to do the same.

new ml=1e3, i=ml, r, a=1, ee=1, ans=1, fbr=6, pr=".############", tg=num(atan(1)*4,pr)
{
r=random(1); r=(r*2-1)*ee+ans; if abs(sin(r/6)-.5) to e LT a then a=e; ans=r; ee=min(e*fbr,ee)
if num(ans,pr)=tg then exit
if --i GT 0 then repeat
$ Unable to complete in `ml` loops
}
$ Est of pi using converging random process to `len(pr)-1` decimal places in `ml-i+1` loops

Some Results

loop count: 92, 85, 80, 94, 60, 89, 90

------------------
www.pegge.net

posted February 01, 2007 12:51 PM            

I don't know your scripting language, & my attempt to translate
your program to PBCC 4.02 code compiled OK but is obviously in error.

The program looks neat. Of course, it ain't "evolution" except in
the original non-Darwinian sense, of convergence on the archetype.

Speaking of which:
How would you rate the chances that a truly Darwinian simulation can
compute the value of pi to as many decimal places as (say)
mathematicians can as of today? I am "morally certain" -- to use the
language of Poincare if not also Pascal -- that the probability is
in effect zero, and I would accordingly bet "my life, soul or even
a year's salary" on it, as physicist Leonard Susskind worded such
wagers more recently... But hey, what about a minimum bet, filed
with LongBets.org?

Kind readers:
Can you help with the following source code?

#COMPILE EXE
#DIM ALL
FUNCTION PBMAIN () AS LONG
    DIM m1 AS EXT, i AS QUAD, r AS EXT, a AS EXT, ee AS EXT, ans AS EXT
    DIM fbr AS EXT, pi AS EXT, e AS EXT, ml AS QUAD
    DIM pr AS STRING, tg AS STRING, seed AS EXT
    ml=1e3: i=ml: r=0: a=1: ee=1: ans=1: fbr=6: pi=ATN(1)*4: e=EXP(1)
    pr="#.############"
    tg=USING$(pr,pi)
    PRINT "Successive approximation using random numbers"
    PRINT "Target: ";tg
    PRINT
    seed=TIMER
    RANDOMIZE seed
    DO
        r=RND(1)
        r=(r*2-1)*ee+ans
        IF ABS(SIN(r/6)-.5)^ e < a THEN
            a=e
            ans=r
            IF ee > e*fbr THEN ee= e*fbr
        END IF
        IF USING$(pr, ans)=tg THEN EXIT DO
        i=i-2
    LOOP UNTIL i <= 0
    IF i<=0 THEN
        PRINT "Unable to complete in "; ml; "loops"
        PRINT ans
    ELSE
        PRINT "Hit target Est of pi to "; pr "decimal places in "; ml-i+1; "loops"
    END IF
    WAITKEY$
END FUNCTION