Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
Announcements
Welcome to the Intel Community. If you get an answer you like, please mark it as an Accepted Solution to help others. Thank you!
26749 Discussions

## Bayesian Statistics - standard example and Fortran

Valued Contributor II
213 Views

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

29 Replies
Honored Contributor II
187 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

Black Belt Retired Employee
187 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.

Black Belt
187 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

Black Belt
187 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

Valued Contributor II
187 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.

Black Belt Retired Employee
187 Views

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

Valued Contributor II
187 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 --

Black Belt Retired Employee
187 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.

New Contributor II
187 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```

Black Belt Retired Employee
187 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.

Valued Contributor II
187 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

Valued Contributor II
187 Views

Steve:

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

John

Black Belt Retired Employee
187 Views

I can't reproduce this error.

Black Belt Retired Employee
187 Views

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

Employee
187 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.

Black Belt
187 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

Honored Contributor II
187 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.

Valued Contributor II
187 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.

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

Valued Contributor II
187 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.

Black Belt Retired Employee
117 Views

I disagree significantly with your characterization of the Fortran RANDOM_NUMBER procedure. The old Microsoft generator you refer to was a simplistic "multiplicative congruential" generator with many flaws. Modern Fortran compilers use much better generators with significantly longer periods and uniformity.