Skip to content

Files

Latest commit

 

History

History
239 lines (180 loc) · 8.53 KB

random_number.md

File metadata and controls

239 lines (180 loc) · 8.53 KB

random number

There's actually no way for software to produce a truly random number.

Pseudo random numbers: a series of numbers in a particular range that seems unpredicatable but is actually produced by a straightforward mathematical process.

Linear Congruential Generator

  • Linear Congruential Generator has 4 inputs:

    1. modulus m , 0<m
    2. multiplier a , 0<=a<m
    3. increment c , 0<=c<m
    4. seed X , starting value , 0<=X<m
  • linear congruential sequence

    • Xn+1 = ( a * Xn + c ) mod m , n>=0
    • for example, the sequence obtained when m = 10 and X₀ = a = c = 7 is
      • 7,6,9,0,7,6,9,0...
    • As this example shows, the sequence is not always "random" for all choices of m,a,c,and X.
  • Above example illustrates the fact that the congruential sequences always get into a loop.

    • This property is common to all sequences having the general form Xn+1 = f( Xn ) .
    • The repeating cycle is called the period. A useful sequence will of course have a relatively long period.
  • if we select good values for m , a ,c , then the period will be m.

Choice of modulus

  • We want m to be rather large, since the period can not have more than m elements.
    • Even if we intend to generate only randOm zeros and ones, we should not take m = 2, for then the sequence would at best have the form..., 0,1,0,1,0,1,...
  • Another factor that inffuences our choice of m is speed of generation:
    • We want to pick a value so that the computation of (aXn + c ) mod m is fast.
    • like 2³².
    • Another alternative is to let m be the largest prime number less than w .

Choice of multiplier

  • Theorem A. The linear congruential sequence defined by m, a, c, and X₀ has period length m if and only if // 序列有周期长度m当且仅当

    1. c is relatively prime to m; // c 与 m 互质
    2. b = a-1, is a multiple of p , for every prime p dividing m; // 对于整除m的每个素数p,b=a-1是p的倍数;
    3. b is a multiple of 4, if m is a multiple of 4. // 如果m是4的倍数,则b也是4的倍数
  • if m is the word size zᵉ , the multiplier

    • a = zᵏ + 1, 2 ≤ k ≤ e
    • satisfies these conditions.
    • Note: z = 2 for a binary computer, and z = 10 for a demical compute
  • We may take c=1,

    • Xn+1 = ( ( zᵏ + 1 ) Xn + 1 ) mod zᵉ
    • then we can avoid multiplication ; merely shifting and adding will suffice.
      • e.g, for z=2, X = ((X<<k) + X + 1) & ( zᵉ - 1 )

An implementation : No guaranteed of correctness

# python3
class Random(object):
    def __init__(self, e, x=None):
        if x is None:
            import time

            self.seed = int(time.time())
        else:
            self.seed = x
        self._seed = self.seed

        self.e = e
        self.k = (self.e + 1) // 2  # 2 <= k <= e
        self.m = 2**self.e
        # c could be any prime number, because is always coprime to m (power of 2)
        self.c = 11
        # self.a = (2** self.k) + 1  # use bit shift instead

    def random_LCG(self):
        """
        generate ranom in 0 ~ self.m-1

        X(k) = [a * X(k-1) + c] mod m
        """
        # m,a,c is referenced in GCC compiler
        # self._seed = (self.a * self._seed + 1) & ( self.m-1 )
        self._seed = (self._seed + (self._seed << self.k) + self.c) & (self.m - 1)
        return self._seed


for e in range(8, 20, 1):
    r = Random(e)
    d = {r.random_LCG() for _ in range(2**e)}
    print("e:{} , loop: {}".format(e, len(d)))
    assert len(d) == 2**e
# a simple implementation
class LCG():
    def __init__( self, seed, mul, inc , mod ):
        self.seed = seed
        self.mul = mul
        self.inc = inc
        self.mod = mod

    def next(self):
        self.seed = (( self.mul * self.seed + self.inc ) % self.mod)
        return self.seed


if __name__ == "__main__":
    lcg = LCG( 7,7,7, 10  )
    print( [ lcg.next() for i in range(16) ] )

    lcg = LCG(7, 5, 1, 16)
    print([lcg.next() for i in range(16)])
    # [4, 5, 10, 3, 0, 1, 6, 15, 12, 13, 2, 11, 8, 9, 14, 7]

    lcg = LCG(7, 5, 3, 16)
    print([lcg.next() for i in range(16)])
    # [6, 1, 8, 11, 10, 5, 12, 15, 14, 9, 0, 3, 2, 13, 4, 7]

    lcg = LCG(9, 5, 91, 16)
    print([lcg.next() for i in range(16)])
    # [8, 3, 10, 13, 12, 7, 14, 1, 0, 11, 2, 5, 4, 15, 6, 9]

    # NOTE: the last number of sequence always be the (initial seed % m)

Python random implementation

def random(self):
    """Get the next random number in the range [0.0, 1.0)."""

    # Wichman-Hill random number generator.
    #
    # Wichmann, B. A. & Hill, I. D. (1982)
    # Algorithm AS 183:
    # An efficient and portable pseudo-random number generator
    # Applied Statistics 31 (1982) 188-190
    #
    # see also:
    #        Correction to Algorithm AS 183
    #        Applied Statistics 33 (1984) 123
    #
    #        McLeod, A. I. (1985)
    #        A remark on Algorithm AS 183
    #        Applied Statistics 34 (1985),198-200

    # This part is thread-unsafe:
    # BEGIN CRITICAL SECTION
    x, y, z = self._seed
    x = (171 * x) % 30269
    y = (172 * y) % 30307
    z = (170 * z) % 30323
    self._seed = x, y, z
    # END CRITICAL SECTION

    # Note:  on a platform using IEEE-754 double arithmetic, this can
    # never return 0.0 (asserted by Tim; proof too long for a comment).
    return (x/30269.0 + y/30307.0 + z/30323.0) % 1.0
  • Wichman-Hill算法通过线性合并不同短周期随机数发生器的输出来产生长周期的随机数序列。
    • 如果把周期N1和N2的波形加起来那么得到的波形周期为 N = lcm(N1,N2)
    • 这样,通过结合几个随机数发生器的输出,可以产生一个更长的序列。

Distribution

  • It's much more useful to have a sequence of numbers with the uniform distribution on the interval 0 to 1.
  • This is useful because with a sequence of numbers with a uniform distribution on [0,1], we can generate sequences of random numbers according to other distributions use something known as inverse transform sampling.
  • Let's say we want to generate a sequence of random numbers with the normal distribution with mean=162cm, and standard deviation=5cm , to represent the height of an American woman.
    • Probability Density Function --> Cumulative Distribution Function.
    • we uniformly randomly select a point on the y-axis ( Cumulative Distribution Function graph ) , and then determine the point on the x-axis, it gives that function value.

Python normalvariate implementation

    from math import log as _log, exp as _exp, e as _e
    from math import sqrt as _sqrt

    # NV_MAGICCONST = 4 * _exp(-0.5) / _sqrt(2.0)
    # self.assertAlmostEqual(random.NV_MAGICCONST, 1.71552776992141)
    # 4 * (math.e**-0.5) / math.sqrt(2.0)
    # 1.7155277699214135
    def normalvariate(self, mu, sigma):
        """Normal distribution.
        mu is the mean, and sigma is the standard deviation.
        """
        # mu = mean, sigma = standard deviation

        # Uses Kinderman and Monahan method. Reference: Kinderman,
        # A.J. and Monahan, J.F., "Computer generation of random
        # variables using the ratio of uniform deviates", ACM Trans
        # Math Software, 3, (1977), pp257-260.

        random = self.random
        while 1:
            u1 = random()
            u2 = 1.0 - random()
            z = NV_MAGICCONST*(u1-0.5)/u2
            zz = z*z/4.0
            if zz <= -_log(u2):  # math.e
                break
        return mu + z*sigma

Book "Probability Theory and Examples"