Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.

Random Numbers and Buildings

JohnNichols
Valued Contributor III
782 Views

There is a nice set of Daubechies wavelets coded in old Fortran by John Burkardt from Pittsburg.  I have a building that is shaking itself to pieces slowly, but surely. I was thinking of looking at the various acceleration signals from the building using these wavelets.  A blasted long shot I am sure. 

I enclose a zip of the code changed to F90.  I started with a 1600 long random set with D2.  The issue is firstly the input, as shown in the figure generated using the code random method. 

SI1.png

The mean is 0.496, the stdev is 0.289,   0.5/0.289 is 1.73 so it is a truncated Gaussian, the skew is 0.00176 and the kurt is -1.17 

For instance, the uniform distribution (i.e. one that is uniformly finite over some bound and zero elsewhere) is platykurtic. As this one is. 

And if I look at the first 50 numbers as a graph it is clearly a simple Fourier series. 

 

The code for the RNG, it is clearly not really random. 

 

!*********************************************************************72
!
! R8VEC_UNIFORM_01 returns a unit pseudorandom R8VEC.
!
!  Discussion:
!
!    An R8VEC is a vector of R8's.
!
!
!  Reference:
!
!    Paul Bratley, Bennett Fox, Linus Schrage,
!    A Guide to Simulation,
!    Springer Verlag, pages 201-202, 1983.
!
!    Bennett Fox,
!    Algorithm 647:
!    Implementation and Relative Efficiency of Quasirandom
!    Sequence Generators,
!    ACM Transactions on Mathematical Software,
!    Volume 12, Number 4, pages 362-376, 1986.
!
!    Peter Lewis, Allen Goodman, James Miller,
!    A Pseudo-Random Number Generator for the System/360,
!    IBM Systems Journal,
!    Volume 8, pages 136-143, 1969.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input/output, integer SEED, the "seed" value, which should NOT be 0.
!    On output, SEED has been updated.
!
!    Output, double precision R(N), the vector of pseudorandom values.
!
      implicit none

      integer n

      integer i
      integer k
      integer seed
      double precision r(n)

      do i = 1, n

        k = seed / 127773

        seed = 16807 * ( seed - k * 127773 ) - k * 2836

        if ( seed .lt. 0 ) then
          seed = seed + 2147483647
        end if

        r(i) = dble ( seed ) * 4.656612875D-10

      end do

      return
      end

JohnNichols_0-1747422600907.png

A Fourier series is not a random set, it cannot be, even if you pass the random number test. 

 

 

@mecej4  or @jimdempseythecove or one of the other gurus, please tell me why I am wrong?

 

 

Below is the D2 results. 

JohnNichols_1-1747422808202.png

I can change to the MKL randoms, but I cannot see that the result over a bounded width willnot be the same,  I will always have a platykurtic. distribution as we do from the acceleration data and FFT's for some components. 

 

 

 

0 Kudos
7 Replies
JohnNichols
Valued Contributor III
719 Views

You see, one goes along for years on this forum, and the same people provide the guiding light for a poor programmer.  And then they just go quiet.  Even though you have never met them, it leaves a hole in your soul.  

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

The problem with using Fortran as a RNG is that the language has strict structure, it is designed to solve engineering problems, not generate RNGs.  It repeats in RNG's, there is always a repeat, this has plagued Fortran for the last 80 years.  Basica always had a better built in random number generator, one only has to look at the most basic Genetic Algorithm from the SciAmer in the 1980's to see that.

A real RNG from nature is thermal or similar,  in the picture everything up to 2.5 Hz is purely random and Gaussian, you cannot get any information on the underlying structure, which is what you want in an RNG, but above that we can determine structure and hence it is not a RNG generator.  The thermal line in the FFT data is always the same shape up to a point of structural interaction. 

Histogram-Z-FFT-LOG.png

I realize this is not about the Intel Fortran, but an underlying issue. 

 

 

0 Kudos
jimdempseyatthecove
Honored Contributor III
624 Views

John,

If you are going to roll your own RNG, try something along the lines of this:

program Prime
    implicit none
    integer(8) :: seed = 0, i8temp
    integer(8), parameter :: largeCircularPrime = 999331_8, mask=Z'FFFFFFFF'
    integer, parameter :: n=10000
    real(8) r(n), rTemp, rDivisor
    integer :: i
    rDivisor = mask + 1
    do i=1,n
        if(seed == 0) seed = 1
        seed = seed * largeCircularPrime
        i8temp = iand(ishft(seed, 16), mask) ! random number between 2^32-1 and 0
        rTemp = real(i8temp,8)
        r(i) = rTemp / rDivisor ! r(i) >= 0, r(i) < 1.0
    end do
end program Prime

Note, the above written for clarity of understanding not performance. You can choose a different largeCircular prime and/or change the mask width (with appropriate ishft bits).  IOW the above produces pseudo random numbers of 32 bits of precision.

To get more (say 48 bits), widen the mask by 16 bits and reduce the shift by 8 bits.

Jim Dempsey

0 Kudos
JohnNichols
Valued Contributor III
592 Views

JohnNichols_0-1747858638937.png

Your method has an underlying structure that repeats, the lines should all be along the bottom as a solid block, the popups are the repeats, this is in 4096 numbers, that may appear random but are not.  If you use this series looking for a Monte Carlo Analysis you will run into problems with the structure.  First 26 points show the structure.  

 

JohnNichols_1-1747858809354.png

If I measure a vibrating string and I get a non mode frequency < 1 or 2 depending on your luck, you have a random Gaussian, if the vibration is from the sun.  If you have an underlying math structure I can find the pattern.  It looks random, it may even pass a so called random test, but it is not truly random.

 

 

0 Kudos
JohnNichols
Valued Contributor III
591 Views

If you run Gardners genetic algorithm with Fortran and any RNG it always fails on the ones I have tested.

BASICA however does not fail, I just do not know the form of their RNG. 

0 Kudos
jimdempseyatthecove
Honored Contributor III
530 Views

>> the popups are the repeats, this is in 4096 numbers

Thanks. This would warrant investigation. 

What did you use to chart the pop-ups?

This may be due to bad choice of multiplier 999331, can you try a non-circular prime number.

Here is a list of different types of prime numbers:

https://en.wikipedia.org/wiki/List_of_prime_numbers

Try a balanced prime first. (e.g. 5387)

 

Jim

0 Kudos
JohnNichols
Valued Contributor III
500 Views

The plotting is on EXCEL, for simple stuff it is easier.  

My problem is I look for patterns in data.  Random is never random.

So your system gives a usable set of a Monte Carlo analysis, but it is a long way from truly random.  

 

If I take a bridge I have, it has a theoretical lowest frequency of 3 Hz, we find that nicely if we listen for an hour.   Non Gaussian in a Gaussian stream.

After three years of listening, I found it had a frequency of 1 Hz, that occurred 10 times in 3 million records.   If I pull out those ten records the 1 Hz data is Gaussian Random, but once a month it is non Gaussian.   10 in 3 million is a low probability but for my stuff it matters.  We now know we have low frequencies that are hard to find.  

A prime number generator is a set of circular steps,  anyone who says Prime numbers are random, really does not understand the relationship between the primes and the integer number line.  

Think of this number system

2                                                  1

3                                                10

5                                               100

7                                             1000

11                                         10000

So an composite number can come from the second list 1, 10, 100 etc. it is just a weird multiplication, one could say the next prime is 100000 and it is simple to find. 

so remove all factors of 2, then 3 pops up, remove all factors of 3 and 6 etc. this goes all the way.  There should be a representation system that is not the integer number line that makes it easy.   So 49 is just 2000, 9 is 20 etc. , 6 is just 11 and 25 is 200

0 Kudos
JohnNichols
Valued Contributor III
500 Views

This simple set proves the primes are infinite, but they get further apart in a non random way. 

0 Kudos
Reply