Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Highlighted
New Contributor II
66 Views

Bayesian Statistics - standard example and Fortran

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

0 Kudos
29 Replies
Highlighted
Honored Contributor I
55 Views

OP seems to referring to the numerical approximation of PI using a Monte Carlo simulation: OP can refer to the links below for further details:

http://mathfaculty.fullerton.edu/mathews/n2003/montecarlopimod.html

http://mathfaculty.fullerton.edu/mathews/n2003/MonteCarloPiProg.html

 

0 Kudos
Highlighted
Black Belt Retired Employee
55 Views

Amusingly, this was the same approximation I turned into a coarray example for Intel Fortran before I retired. It only recently made its way into the samples bundle, though some of the explanatory text was omitted (I have asked that it be reinstated,) I've attached the completed example.

0 Kudos
Highlighted
55 Views

The issue may have more to do with the random number generator. Using a variation of: http://www.unm.edu/~hguo/MC.f

     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),allocatable :: x(:,:), r(:)
        real(8) :: pi
        allocate(x(n,2), r(n))
        call random_number(x)
        r = x(:,1)**2 + x(:,2)**2
        pi = 4.0_8 * dble( count( r <= 1.0_8 ) ) / n
        deallocate(r, x)
        end function

          57    63000000
          57    63000000   3.14145295238095      -9.498310291844447E-005
   3.14152921801776
          58    64000000
          58    64000000   3.14164743750000       1.944851190476271E-004
   3.14158832775888
          59    65000000
          59    65000000   3.14145310769231      -1.943298076922950E-004
   3.14152071772559
          60    66000000
          60    66000000   3.14185769696970       4.045892773891779E-004
   3.14168920734765

Jim Dempsey

0 Kudos
Highlighted
55 Views

This appears to oscillate about PI depending on the number of random numbers as well as values thereof.

If we change the calculation of averpi to become the true average of all tests:

     program mc
        integer :: n,i
        real(8) :: pi, pi1, pi2, averpi,sumpi
        n=10000
        pi1 = 0.d0
        averpi = 0.d0
        do i=1,60
            print*,i,n
            pi2 = pi(n)
            sumpi = sumpi + pi2
            averpi = sumpi / i
          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),allocatable :: x(:,:), r(:)
        real(8) :: pi
        allocate(x(n,2), r(n))
        call random_number(x)
        r = x(:,1)**2 + x(:,2)**2
        pi = 4.0_8 * dble( count( r <= 1.0_8 ) ) / n
        deallocate(r, x)
        end function

 

          57    63000000
          57    63000000   3.14145295238095      -9.498310291844447E-005
   3.14141584646072
          58    64000000
          58    64000000   3.14164743750000       1.944851190476271E-004
   3.14141983940967
          59    65000000
          59    65000000   3.14145310769231      -1.943298076922950E-004
   3.14142040327887
          60    66000000
          60    66000000   3.14185769696970       4.045892773891779E-004
   3.14142769150705

And if we increase the DO I=1,600 and reduce advancement of n to n-n+100000 (~same number of random numbers but higher number of samples into averpi:

     program mc
        integer :: n,i
        real(8) :: pi, pi1, pi2, averpi,sumpi
        n=10000
        pi1 = 0.d0
        averpi = 0.d0
        do i=1,600
            print*,i,n
            pi2 = pi(n)
            sumpi = sumpi + pi2
            averpi = sumpi / i
          print*,i,n,pi2,(pi2-pi1),averpi
          pi1=pi2
          if(i >= 4)then
              n=n + 100000
          else
          n = n * 10
          endif
        end do
        end program
 
        function  pi(n)
        integer :: n
        real(8),allocatable :: x(:,:), r(:)
        real(8) :: pi
        allocate(x(n,2), r(n))
        call random_number(x)
        r = x(:,1)**2 + x(:,2)**2
        pi = 4.0_8 * dble( count( r <= 1.0_8 ) ) / n
        deallocate(r, x)
        end function

Now the results data shows by increasing the number calculations (as well as accumulating increasing random number harvests) that the value appears to converge PI:

         501    59700000   3.14156247906198       4.962283453036065E-006
   3.14155138619635
         502    59800000
         502    59800000   3.14115625418060      -4.062248813747615E-004
   3.14155059908078
         503    59900000
         503    59900000   3.14150824707846       3.519928978623632E-004
   3.14155051488197
         504    60000000
         504    60000000   3.14157133333333       6.308625486894925E-005
   3.14155055618842
         505    60100000
         505    60100000   3.14191134775374       3.400144204106503E-004
   3.14155127062716
         506    60200000
         506    60200000   3.14131880398671      -5.925437670328826E-004
   3.14155081120692
         507    60300000
         507    60300000   3.14175582089552       4.370169088114828E-004
   3.14155121556528
         508    60400000
         508    60400000   3.14175642384106       6.029455370537562E-007
   3.14155161951858
         509    60500000
         509    60500000   3.14144773553719      -3.086883038694310E-004
   3.14155141542432
         510    60600000
         510    60600000   3.14169828382838       2.505482911927537E-004
   3.14155170340158
         511    60700000
         511    60700000   3.14186919275124       1.709089228527638E-004
   3.14155232471146
         512    60800000
         512    60800000   3.14141394736842      -4.552453828146064E-004
   3.14155205444321
         513    60900000
         513    60900000   3.14161576354680       2.018161783770012E-004
   3.14155217863250
         514    61000000
         514    61000000   3.14151711475410      -9.864879269949611E-005
   3.14155211041484
         515    61100000
         515    61100000   3.14171554828151       1.984335274070048E-004
   3.14155242776992
         516    61200000
         516    61200000   3.14171568627451       1.379930041345290E-007
   3.14155274416237
         517    61300000
         517    61300000   3.14145324632953      -2.624399449828729E-004
   3.14155255171008
         518    61400000
         518    61400000   3.14140351791531      -4.972841421757579E-005
   3.14155226400005
         519    61500000
         519    61500000   3.14178217886179       3.786609464793145E-004
   3.14155270699593
         520    61600000
         520    61600000   3.14174155844156      -4.062042023011969E-005
   3.14155307017179
         521    61700000
         521    61700000   3.14156434359806      -1.772148435033039E-004
   3.14155309180984
         522    61800000
         522    61800000   3.14183365695793       2.693133598734576E-004
   3.14155362929097
         523    61900000
         523    61900000   3.14175185783522      -8.179912271044643E-005
   3.14155400831304
         524    62000000
         524    62000000   3.14201154838710       2.596905518785775E-004
   3.14155488148112
         525    62100000
         525    62100000   3.14179336553945      -2.181828476444103E-004
   3.14155533573647
         526    62200000
         526    62200000   3.14144372990354      -3.496356359153907E-004
   3.14155512355808
         527    62300000
         527    62300000   3.14177534510433       3.316152007970530E-004
   3.14155554143578
         528    62400000
         528    62400000   3.14146858974359      -3.067553607443152E-004
   3.14155537675454
         529    62500000
         529    62500000   3.14142860800000      -3.998174358965656E-005
   3.14155513711607
         530    62600000
         530    62600000   3.14188146964856       4.528616485623971E-004
   3.14155575283783
         531    62700000
         531    62700000   3.14124210526316      -6.393643854045727E-004
   3.14155516216443
         532    62800000
         532    62800000   3.14089980891720      -3.422963459605022E-004
   3.14155393029742
         533    62900000
         533    62900000   3.14176903020668       8.692212894798601E-004
   3.14155433386198
         534    63000000
         534    63000000   3.14179739682540       2.836661871974400E-005
   3.14155478903607
         535    63100000
         535    63100000   3.14159156893819      -2.058278872034691E-004
   3.14155485778355
         536    63200000
         536    63200000   3.14161588607595       2.431713775585820E-005
   3.14155497164230
         537    63300000
         537    63300000   3.14136859399684      -2.472920791087851E-004
   3.14155462457034
         538    63400000
         538    63400000   3.14162302839117       2.544343943267080E-004
   3.14155475171499
         539    63500000
         539    63500000   3.14145121259843      -1.718157927421693E-004
   3.14155455962015
         540    63600000
         540    63600000   3.14170465408805       2.534414896251391E-004
   3.14155483757287
         541    63700000
         541    63700000   3.14148169544741      -2.229586406405915E-004
   3.14155470237485
         542    63800000
         542    63800000   3.14176858934169       2.868938942830468E-004
   3.14155509700025
         543    63900000
         543    63900000   3.14227618153365       5.075921919535098E-004
   3.14155642496440
         544    64000000
         544    64000000   3.14156506250000      -7.111190336464013E-004
   3.14155644084223
         545    64100000
         545    64100000   3.14176162246490       1.965599648987570E-004
   3.14155681732227
         546    64200000
         546    64200000   3.14168186915888      -7.975330602016939E-005
   3.14155704635494
         547    64300000
         547    64300000   3.14112989113530      -5.519780235752059E-004
   3.14155626544960
         548    64400000
         548    64400000   3.14122776397516       9.787283985218664E-005
   3.14155566599435
         549    64500000
         549    64500000   3.14159503875969       3.672747845344659E-004
   3.14155573771159
         550    64600000
         550    64600000   3.14141578947368      -1.792492860057671E-004
   3.14155548326025
         551    64700000
         551    64700000   3.14194911901082       5.333295371352520E-004
   3.14155619766270
         552    64800000
         552    64800000   3.14199370370370       4.458469288426414E-005
   3.14155699024611
         553    64900000
         553    64900000   3.14173442218798      -2.592815157220763E-004
   3.14155731109953
         554    65000000
         554    65000000   3.14167132307692      -6.309911105839561E-005
   3.14155751689733
         555    65100000
         555    65100000   3.14143797235023      -2.333507266927271E-004
   3.14155730150175
         556    65200000
         556    65200000   3.14131441717791      -1.235551723164363E-004
   3.14155686465944
         557    65300000
         557    65300000   3.14161941807044       3.050008925300496E-004
   3.14155697696359
         558    65400000
         558    65400000   3.14161779816514      -1.619905306338154E-006
   3.14155708596215
         559    65500000
         559    65500000   3.14107572519084      -5.420729742979802E-004
   3.14155622485165
         560    65600000
         560    65600000   3.14187615853659       8.004333457458657E-004
   3.14155679616180
         561    65700000
         561    65700000   3.14139445966514      -4.816988714408943E-004
   3.14155650679193
         562    65800000
         562    65800000   3.14117282674772      -2.216329174240528E-004
   3.14155582408723
         563    65900000
         563    65900000   3.14190822458270       7.353978349806667E-004
   3.14155645002061
         564    66000000
         564    66000000   3.14169696969697      -2.112548857313712E-004
   3.14155669916898
         565    66100000
         565    66100000   3.14155558245083      -1.413872461379917E-004
   3.14155669719248
         566    66200000
         566    66200000   3.14143166163142      -1.239208194121488E-004
   3.14155647628160
         567    66300000
         567    66300000   3.14182063348416       3.889718527432251E-004
   3.14155694216732
         568    66400000
         568    66400000   3.14175524096386      -6.539252030757225E-005
   3.14155729128492
         569    66500000
         569    66500000   3.14178983458647       3.459362261093446E-005
   3.14155769997262
         570    66600000
         570    66600000   3.14157087087087      -2.189637155955992E-004
   3.14155772307946
         571    66700000
         571    66700000   3.14135580209895      -2.150687719200128E-004
   3.14155736945252
         572    66800000
         572    66800000   3.14169299401198       3.371919130255030E-004
   3.14155760655839
         573    66900000
         573    66900000   3.14172000000000       2.700598802363174E-005
   3.14155788996754
         574    67000000
         574    67000000   3.14163217910448      -8.782089552239469E-005
   3.14155801939113
         575    67100000
         575    67100000   3.14129794336811      -3.342357363700366E-004
   3.14155756708500
         576    67200000
         576    67200000   3.14172065476190       4.227113937971261E-004
   3.14155785022333
         577    67300000
         577    67300000   3.14164172362556      -7.893113634738214E-005
   3.14155799558451
         578    67400000
         578    67400000   3.14160243323442      -3.929039113570454E-005
   3.14155807246626
         579    67500000
         579    67500000   3.14176568888889       1.632556544675090E-004
   3.14155843104384
         580    67600000
         580    67600000   3.14186118343195       9.549454306378991E-005
   3.14155895303072
         581    67700000
         581    67700000   3.14156519940916      -2.959840227947375E-004
   3.14155896378180
         582    67800000
         582    67800000   3.14180424778761       2.390483784524555E-004
   3.14155938523198
         583    67900000
         583    67900000   3.14140618556701      -3.980622206003481E-004
   3.14155912245383
         584    68000000
         584    68000000   3.14144400000000       3.781443298978004E-005
   3.14155892532634
         585    68100000
         585    68100000   3.14157879588840       1.347958883997080E-004
   3.14155895929311
         586    68200000
         586    68200000   3.14177489736070       1.961014723041998E-004
   3.14155932778811
         587    68300000
         587    68300000   3.14166407027818      -1.108270825191937E-004
   3.14155950622506
         588    68400000
         588    68400000   3.14144701754386      -2.170527343250406E-004
   3.14155931491778
         589    68500000
         589    68500000   3.14169191240876       2.448948648994254E-004
   3.14155954004085
         590    68600000
         590    68600000   3.14144682215743      -2.450902513246866E-004
   3.14155934899359
         591    68700000
         591    68700000   3.14174358078603       2.967586285920198E-004
   3.14155966072251
         592    68800000
         592    68800000   3.14143209302326      -3.114877627705681E-004
   3.14155944523653
         593    68900000
         593    68900000   3.14139297532656      -3.911769669562304E-005
   3.14155916451156
         594    69000000
         594    69000000   3.14175124637681       3.582710502514352E-004
   3.14155948788170
         595    69100000
         595    69100000   3.14151068017366      -2.405662031503830E-004
   3.14155940585194
         596    69200000
         596    69200000   3.14147075144509      -3.992872857461194E-005
   3.14155925710294
         597    69300000
         597    69300000   3.14132756132756      -1.431901175252825E-004
   3.14155886900281
         598    69400000
         598    69400000   3.14180749279539       4.799314678276545E-004
   3.14155928476166
         599    69500000
         599    69500000   3.14166112230216      -1.463704932307408E-004
   3.14155945477425
         600    69600000
         600    69600000   3.14159408045977      -6.704184238826016E-005
   3.14155951248372
Press any key to continue . . .

Using r = hypot(x(:,1), x(:2))

         597    69300000
         597    69300000   3.14132756132756      -1.431901175252825E-004
   3.14155886900281
         598    69400000
         598    69400000   3.14180749279539       4.799314678276545E-004
   3.14155928476166
         599    69500000
         599    69500000   3.14166112230216      -1.463704932307408E-004
   3.14155945477425
         600    69600000
         600    69600000   3.14159408045977      -6.704184238826016E-005
   3.14155951248372

Jim Dempsey

0 Kudos
Highlighted
New Contributor II
55 Views

Dear Jim and Steve:

Thanks for the feedback.  I tried the coarrays, but I kept getting an error -- shown on Capture 1 -- I have included the zip file with the solution and the code -- I cannot see how to fix it - nor does the error pop up on a google search -- I could not see anything in the tutorial example. 

The MC code makes a difference and yes it appears to converge on the average nicely, Jim - thanks. 

 

 

0 Kudos
Highlighted
Black Belt Retired Employee
55 Views

For the Coarray project you need to set Fortran > Language > Use Coarrays to Yes.

0 Kudos
Highlighted
New Contributor II
55 Views

! Compiler options: /Qcoarray

! -coarray

!

I set it this way and with turning on the language feature which includes a shared -- word.  It still gave the error. 

The zip file in the previous post gives the complete code - it is a direct copy of the tutorial -- I am sorry but I cannot see any mistakes --

 

0 Kudos
Highlighted
Black Belt Retired Employee
55 Views

I took your Zip, turned on the Coarray feature (which was off in your project), built and ran it without errors. You don't need the command line option if you are using Visual Studio.

0 Kudos
Highlighted
New Contributor II
55 Views

This does not look like a good approximation to pi, as it does not appear to converge. I wrote a simple adaptation to show how the error varies as the number of random number pairs is increased, but there is no apparent convergence. The error typically remains at about 1.e-5. not much better than "22 on 7"

Or have I got this wrong ?

     integer*4 i,j,m
     integer*8 n, nc
     real*8    pii, pi, err, x(2), rr
     
     pii = 4 * atan (1.0d0)
     m  = 10000
     n  = 0
     nc = 0
     do i = 1,m
       do j = 1,m
         call random_number (x)
         rr = (x(1)*x(1)+x(2)*x(2))
         if ( rr <= 1 ) nc = nc+1
         n = n+1
       end do
       pi = 4.0d0 * dble(nc) / dble(n)
       err = pii-pi
       if ( m-i < 10 .or. abs (err) < 1.d-6 ) &
       write (*,*) i, pi, err
     end do
     end

 

0 Kudos
Highlighted
Black Belt Retired Employee
55 Views

That's the approach I took. It doesn't converge very fast, I'll grant you. I suspect part of the issue may be the random number generator - if it returns the space sparsely filled, you'll just be redoing the same calculations over time. I know that the ifort random_number algorithm is good regarding sequence length, but am not sure how "granular" it is. One might experiment with other generators.

0 Kudos
Highlighted
New Contributor II
55 Views

Dear Steve:

I took your downloaded code, opened a brand new empty project in VS 2015 and another one in VS2017 and loaded your code exactly into the program and I fixed the coarrays as shown in the attached figure and I still get the error shown on the second picture.  I have never seen this error before , unknown option - localonly -- any ideas where it comes from and man I followed your instructions exactly -- this is a stock standard installation of windows and Fortran.

The FORTRAN random number generator is not good, it has been not good since the 1980s -- it repeats to quickly.  

John

0 Kudos
Highlighted
New Contributor II
55 Views

Steve:

I created a blank project -- hello world - turned on coarrays and the error appears

John

0 Kudos
Highlighted
Black Belt Retired Employee
55 Views

I can't reproduce this error. 

0 Kudos
Highlighted
Black Belt Retired Employee
55 Views

Intel Fortran's RANDOM_NUMBER implements L'Ecuyer '91. It was supposed to be a good one at the time.

0 Kudos
Highlighted
55 Views

The "unknown option" is a clue.   Do you have another MPI installed on your system? 

By default, the Intel Fortran installation puts the Intel MPI directory at the end of PATH.   This means that if you have another MPI, it is being found first.

 

0 Kudos
Highlighted
55 Views

    integer*4 i,j,m
     integer*8 n, nc
     real*8    pii, pi, err, x(2), rr, errmin
     
     errmin = 1.0
     pii = 4 * atan (1.0d0)
     m  = 10000
     n  = 0
     nc = 0
     do i = 1,m
       do j = 1,m
         call random_number (x)
         rr = (x(1)*x(1)+x(2)*x(2))
         if ( rr <= 1 ) nc = nc+1
         n = n+1
       end do
       pi = 4.0d0 * dble(nc) / dble(n)
       err = pii-pi
       if(abs(err) < errmin) errmin = abs(err)
       if ( m-i < 10 .or. abs (err) < 1.d-6 ) &
       write (*,*) i, pii, pi, err, errmin
     end do
       10000   3.14159265358979        3.14171116000000
 -1.185064102067201E-004  1.184467972592529E-007

"Know when to hold-em, know when to fold-em" Kenny Rodgers - The Gambler

The best approximation had error of 1.845E-7

Jim Dempsey

0 Kudos
Highlighted
Honored Contributor I
55 Views

Steve Lionel (Ret.) wrote:

That's the approach I took. It doesn't converge very fast, I'll grant you. I suspect part of the issue may be the random number generator ..

It's well-documented in math literature the method in question to approximate the value of PI converges slowly and has an order of error of O( n**(-0.5) ).  The results shown using the PRNG in Intel Fortran appear consistent with this.

0 Kudos
Highlighted
New Contributor II
55 Views

Apologies for being so long to say thank you.  This example was just to get me going on Monte Carlo, but your assistance was great.  The FORTRAN random number generators - no matter which one you use are generally awful compared to the old MICROSOFT Basic random number generator. I base this opinion on trying for many years to code Martin Gardner's excellent genetic program from the late 1980s in Scientific American -- I enclose the basic routine and my attempt at Fortran.  I could never get the Fortran one to work properly and the basic one worked first time.  The FORTRAN's problem appears to be repeats, in this case in the Monte Carlo the numbers are clearly cycling with a set frequency, I do so much FFT work every day you can see it in the number output without putting them into an FFT routine.  I played with a lot of random routines over the years, but never found a satisfactory one for MG's program in Fortran and I really tried as hard as I knew.  The Soup Fortran is written for Microsoft 3 Fortran compiler.

The issue with the MC is not that it is slow, but knowing that the convergence to an answer is reasonable and the answer is reasonable and here it is. 

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

I am so sick of talking to engineers who can only use MATLAB and then wonder why I sit twiddling my thumbs while they spend weeks on problems of limited data size that can be handled in Fortran in a fraction of the time.  I HATE MATLAB it is taking the world of engineers into a cul-de-sac of no return. 

Of course I am biased I like Fortran and LISP.

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

 

 

0 Kudos
Highlighted
New Contributor II
55 Views

Lorri Menard (Intel) wrote:

The "unknown option" is a clue.   Do you have another MPI installed on your system? 

By default, the Intel Fortran installation puts the Intel MPI directory at the end of PATH.   This means that if you have another MPI, it is being found first.

 

Thank you for the reply.  I had 2018 update 4 installed, so I took it off completely, removed all of the directories and other crap and the reinstalled 2019.  I had to run the mpiexec routine then to set the user and the password into the registry and then it worked. 

I have no idea what hydra-services are but the MPI how to install outlines how to do the setup including hydra____, all steps had been completed except register a local user.

It worked a treat once I recreated a new solution as it would not open the old solution.  But it told me that it was doing the analysis on 6 images, but called each image (1) -- is this correct Steve, I have never used images before so I was a bit lost.

0 Kudos