Bayesian Statistics - standard example and Fortran

[JohnNichols]

Dear Gurus:

!  Console1.f90 
!
!  FUNCTIONS:
!  Console1 - Entry point of console application.
!

!****************************************************************************
!
!  PROGRAM: Console1
!
!  PURPOSE:  Entry point for the console application.
!
!****************************************************************************

     program mc
        integer :: n,i
        real(8) :: pi, pi1, pi2, averpi
        n=10000
        pi1 = 0.d0
        averpi = 0.d0
        do i=1,60
            print*,i,n
            pi2 = pi(n)
            averpi = (averpi + pi2)/2.d0
          print*,i,n,pi2,(pi2-pi1),averpi
          pi1=pi2
          if(i >= 4)then
              n=n + 1000000
          else
          n = n * 10
          endif
        end do
        end program
 
        function  pi(n)
        integer :: n
        real(8) :: x(2,n),pi
        call random_number(x)
        pi = 4.d0 * dble( count( hypot(x(1,:),x(2,:)) <= 1.d0 ) ) / n
        end function

There is a standard algorithm to demonstrate Bayesian Statistics using the calculation of PI.  The sample above was taken from a University Professor's website and I lost the link so I cannot attribute it.  But you see the same example at several locations on the WEB.  This follows on from the discussion on the HYPOT in this forum over the last little while.

There are two problems that I see with the algorithm and I was seeking the thoughts of experts,  firstly it consistently appears to favour an overestimation of PI,  The <= implies that a point on the circle exactly has zero radius, but the hypot function implicitly provides a radius to each point, we do not know what it is -- but it exists. 

I played with averaging the value of PI and that seems to help some or one could determine the nearness of the center of the point to the actual line and then adjust accordingly -- I would appreciate any thoughts on the matter.

John

[jimdempseyatthecove]

>>The FORTRAN's problem appears to be repeats

I think that the numbers do not exactly repeat, but rather they appear to cluster.

For a better random number generator

See: https://github.com/jannisteunissen/rng_fortran/blob/master/m_random.f90

https://github.com/jannisteunissen/rng_fortran is the up-level page

Jim Dempsey

[John_Campbell]

I have performed capacity modelling of materials handling systems for many years, using a Monte-Carlo approach. I have used many Fortran compilers and a variety of random number generators. I have never found a "random number" generator to be a significant cause of errors. What can be a problem with less accurate generators is the accuracy of the test "does the system change now" when applied to extreme probability events. However this test can be replaced by calculating "how many cycles till a change occurs". Other fields may require better accuracy of random number generators, but a Monte-Carlo approach has much scope for inaccuracy.

[mecej4]

Nichols, John wrote:
The FORTRAN random number generators - no matter which one you use are generally awful compared to the old MICROSOFT Basic random number generator.

I think that this is a not a well-thought-out conclusion. Monte Carlo methods are, in general, slow. When there is no other feasible algorithm, however, we may use Monte Carlo, and being able to obtain any solution at all makes us rejoice.

You will probably enjoy reading Weisstein's article about estimating π by simulating the Buffon Needle Drop, see http://mathworld.wolfram.com/BuffonsNeedleProblem.html . Note that 1000 needle drops barely gives you three good digits. If you took the same needle(s), a piece of chalk  and a long inelastic string, you could make the measurement by laying out a chiliagon on the floor of a large hall, and probably obtain better accuracy, see https://en.wikipedia.org/wiki/Chiliagon .

Random numbers do have some peculiar patterns, see Marsaglia's "Random numbers fall mainly in the planes", https://apps.dtic.mil/dtic/tr/fulltext/u2/685578.pdf. There are papers that expose related problems, e.g., Brent's "The Myth of Equidistribution for High-Dimensional Simulation", https://arxiv.org/pdf/1005.1320.pdf .

[JohnNichols]

Dear mecej4:

I would normally defer to your knowledge in numerical analysis, which is awesomely extensive, but one program in particular developed from some basic code provided in the beginning by Gardner, and termed soup. He provided the idea in Scientific American, Soup in Microsoft Basic was a great deal of fun to play with and get some interesting insights into the problem.  But translated into Microsoft Fortran and then played with over the years as it is like a nagging uncle or aunt, never performed as well as the Basic code, going into an endless loop rather than a random walk.  From a mathematical perspective you are correct, but from observational data, I would say there is an interesting issue with random numbers and the Fortran implementation. 

So, mathematical you are correct, but I just wish I could solve the Soup problem properly in Fortran it has bugged me since 1986.

John

[mecej4]

Any citations or available links to Gardner's Soup, please ? Even a verbal description would do. Much of Gardner's work was published before we had digital computers.

[JohnNichols]

I enclose the original BASIC program and my lousy Fortran knock off from 1987.  I have played with it a lot over the years and the BASIC program is great fun - I just used to be VERY annoyed with myself that I could not replicate the output in Fortran - the same random walk over and over is not a success.

It comes from Scientific American Recreational Mathematics in about July 1986, I used to keep the copy of the journal next to my bed but I lost it in a move. The BASIC program should run on GWBASIC -- I have a 32 bit version, but am on a 64 bit windows and I do not have the time to spin up a VM at the moment. 

Gardner was absolutely brilliant - along with the guy from Byte - Pournelle --

[JohnNichols]

Back in the early 19th century, Scotland Yard had a group of detectives known as the Big Four. Arthur Ransome in a brilliant children's book, if you have 8 to 12 year old's buying Christmas presents for, then the Big Six by AR will get you a book your child will not put down if they have any interest in science, played to the image in a excellent detective story with six dectiectives.

 I always think of Steve et al as the Big Six of Fortran.

PS I have no involvement in RANSOME other than as a reader.

[jimdempseyatthecove]

John,

Your *.BAS program is a compiled Basic program. Do you have the source?

Jim Dempsey

[mecej4]

jimdempseyatthecove wrote:
Your *.BAS program is a compiled Basic program. Do you have the source?

What he provided is actually a "tokenized" basic source. You can convert it to Ascii source by opening it in GWBasic (or compatible program), and resaving it with the /A option. For your convenience, here is the Ascii source (after saving the file, rename it to "soup1.bas").