4.8. Solutions: Lecture 4

4.8.1. Exercise 4.1

• ex_4.1.py

#!/usr/bin/python

'''
Si scriva una funzione che implementi il generatore lineare congruenziale di numeri pseudo-casuali, 
utilizzando questi parametri:
M = 2147483647
A = 214013
C = 2531011
python3 es_4.1.py 10
'''

import sys


def main () :
    '''
    Funzione che implementa il programma principale
    '''
    if len(sys.argv) < 2 :
        print ('usage: ', sys.argv[0], 'numero')
        exit ()

    # fixed parameters
    A =  214013
    C =  2531011
    M =  2147483647
    # first number, i.e. the seed of the generation
    num = 10       

    randlist = []
    for i in range (int (sys.argv[1])):
        num = (A * num + C) % M
        randlist.append (num)
        print (i, randlist[-1])

    
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 


if __name__ == "__main__":
    main ()

4.8.2. Exercise 4.2

• ex_4.2.py

#!/usr/bin/python

'''
Si scriva una funzione che implementi il generatore lineare congruenziale di numeri pseudo-casuali, 
utilizzando questi parametri:
M = 2147483647
A = 214013
C = 2531011
python3 es_4.2.py 10
'''

import sys


class LCG :
    '''Linear congruential generator'''

    # fixed parameters
    A =  214013
    C =  2531011
    M =  2147483647
    # current sequence value

    def __init__ (self, value = 12345):
        self.num = value

    def seed (self, value) :
        self.num = value

    def random (self) :
        self.num = (self.A * self.num + self.C) % self.M
        return self.num
    

# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 


def main () :
    '''
    Funzione che implementa il programma principale
    '''
    if len(sys.argv) < 2 :
        print ('usage: ', sys.argv[0], 'numero')
        exit ()

    seed = 10
    my_LCG = LCG (seed)

    randlist = []
    for i in range (int (sys.argv[1])):
        randlist.append (my_LCG.random ())
        print (i, randlist[-1])

    
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 


if __name__ == "__main__":
    main ()

4.8.3. Exercise 4.3

• ex_4.3.py

#!/usr/bin/python

'''
Si mostri che inizializzare il seed di un generatore di numeri pseudo-casuali interi
equivale ad inserirsi in una sequenza di numeri pseudo-casuali ad un qualunque punto.
python3 es_4.3.py 10

nb: `random.seed` does not restart the sequence from the given value, 
rather it uses the given value to calculate a random *state* that is then used to generate the random numbers.
'''

import sys
import random


def main () :
    '''
    Funzione che implementa il programma principale
    '''
    if len(sys.argv) < 2 :
        print ('usage: ', sys.argv[0], 'numero')
        exit ()

    randlist = []
    randstates = []
    for i in range (int (sys.argv[1])):
        # Return the next random floating point number in the range 0.0 <= X < 1.0
        randlist.append (random.randint (0, 100))
        randstates.append (random.getstate ())
        print (i, randlist[-1])

    random.setstate (randstates[4])
    print (random.randint (0, 100), randlist[5])

    random.seed( randlist[4] )
    print (random.randint (0, 100), randlist[5])

    
# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 


if __name__ == "__main__":
    main ()

4.8.4. Exercise 4.4

• ex_4.4.py

#!/usr/bin/python
'''
Si implementi un generatore di numeri pseudo-casuali secondo una distribuzione uniforme
fra due estremi arbitrari.
  * Si utilizzi la libreria MatPlotLib per visualizzare la distribuzione dei numeri generati.
python3 es_4.4.py 2.2 4.3 2000 2.4
'''

import sys
import random
import matplotlib.pyplot as plt
import numpy as np
from math import floor

from myrand import generate_range


def main () :
    '''
    Funzione che implementa il programma principale
    '''

    if len(sys.argv) < 5 :
        print ('usage: ', sys.argv[0], 'xMin xMax numero seed')
        exit ()

    xMin = float (sys.argv[1])  # minimum of the histogram drawing range
    xMax = float (sys.argv[2])  # maximum of the histogram drawing range
    seed = float (sys.argv[4])
    N    = int (sys.argv[3])

    print (' -------- ')
    print (' minimum : ', xMin)
    print (' maximum : ', xMax)
    print (' seed    : ', seed)
    print (' N       : ', N)
    print (' -------- ')

    randlist = generate_range (xMin, xMax, N, seed)

    # plotting of the generated list of numbers in a histogram

    nBins = floor (len (randlist) / 20.) + 1     # number of bins of the hitogram
    bin_edges = np.linspace (xMin, xMax, nBins)  # edges o the histogram bins

    # disegno della funzione
    fig, ax = plt.subplots ()
    ax.set_title ('Histogram of random numbers', size=14)
    ax.set_xlabel ('random value')
    ax.set_ylabel ('events in bin')
    ax.hist (randlist,      # list of numbers
             bins = bin_edges,
             color = 'orange',
             # normed = True,
            )
    plt.savefig ('es_4.4.png')


# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 


if __name__ == "__main__":
    main ()

• myrand.py

#!/usr/bin/python

import random
from math import sqrt


def generate_uniform (N, seed = 0.) :
    '''
    generazione di N numeri pseudo-casuali distribuiti fra 0 ed 1
    a partire da un determinato seed
    '''
    if seed != 0. : random.seed (float (seed))
    randlist = []
    for i in range (N):
        # Return the next random floating point number in the range 0.0 <= X < 1.0
        randlist.append (random.random ())
    return randlist


# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 


def rand_range (xMin, xMax) :
    '''
    generazione di un numero pseudo-casuale distribuito fra xMin ed xMax
    '''
    return xMin + random.random () * (xMax - xMin)


# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 


def generate_range (xMin, xMax, N, seed = 0.) :
    '''
    generazione di N numeri pseudo-casuali distribuiti fra xMin ed xMax
    a partire da un determinato seed
    '''
    if seed != 0. : random.seed (float (seed))
    randlist = []
    for i in range (N):
        # Return the next random floating point number in the range 0.0 <= X < 1.0
        randlist.append (rand_range (xMin, xMax))
    return randlist


# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 


def rand_TAC (f, xMin, xMax, yMax) :
    '''
    generazione di un numero pseudo-casuale 
    con il metodo try and catch
    '''
    x = rand_range (xMin, xMax)
    y = rand_range (0, yMax)
    while (y > f (x)) :
        x = rand_range (xMin, xMax)
        y = rand_range (0, yMax)
    return x


# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 


def generate_TAC (f, xMin, xMax, yMax, N, seed = 0.) :
    '''
    generazione di N numeri pseudo-casuali
    con il metodo try and catch, in un certo intervallo,
    a partire da un determinato seed
    '''
    if seed != 0. : random.seed (float (seed))
    randlist = []
    for i in range (N):
        # Return the next random floating point number in the range 0.0 <= X < 1.0
        randlist.append (rand_TAC (f, xMin, xMax, yMax))
    return randlist


# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 


def rand_TCL (xMin, xMax, N_sum = 10) :
    '''
    generazione di un numero pseudo-casuale 
    con il metodo del teorema centrale del limite
    su un intervallo fissato
    '''
    y = 0.
    for i in range (N_sum) :
        y = y + rand_range (xMin, xMax)
    y /= N_sum ;
    return y ;


# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 


def generate_TCL (xMin, xMax, N, N_sum = 10, seed = 0.) :
    '''
    generazione di N numeri pseudo-casuali
    con il metodo del teorema centrale del limite, in un certo intervallo,
    a partire da un determinato seed
    '''
    if seed != 0. : random.seed (float (seed))
    randlist = []
    for i in range (N):
        # Return the next random floating point number in the range 0.0 <= X < 1.0
        randlist.append (rand_TCL (xMin, xMax, N_sum))
    return randlist


# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 


def rand_TCL_ms (mean, sigma, N_sum = 10) :
    '''
    generazione di un numero pseudo-casuale 
    con il metodo del teorema centrale del limite
    note media e sigma della gaussiana
    '''
    y = 0.
    delta = sqrt (3 * N_sum) * sigma
    xMin = mean - delta
    xMax = mean + delta
    for i in range (N_sum) :
        y = y + rand_range (xMin, xMax)
    y /= N_sum ;
    return y ;


# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 


def generate_TCL_ms (mean, sigma, N, N_sum = 10, seed = 0.) :
    '''
    generazione di N numeri pseudo-casuali
    con il metodo del teorema centrale del limite, note media e sigma della gaussiana,
    a partire da un determinato seed
    '''
    if seed != 0. : random.seed (float (seed))
    randlist = []
    delta = sqrt (3 * N_sum) * sigma
    xMin = mean - delta
    xMax = mean + delta
    for i in range (N):
        # Return the next random floating point number in the range 0.0 <= X < 1.0
        randlist.append (rand_TCL (xMin, xMax, N_sum))
    return randlist

4.8.5. Exercise 4.5

• ex_4.4.py

#!/usr/bin/python
'''
Si implementi un generatore di numeri pseudo-casuali secondo una distribuzione uniforme
fra due estremi arbitrari.
  * Si utilizzi la libreria MatPlotLib per visualizzare la distribuzione dei numeri generati.
python3 es_4.4.py 2.2 4.3 2000 2.4
'''

import sys
import random
import matplotlib.pyplot as plt
import numpy as np
from math import floor

from myrand import generate_range


def main () :
    '''
    Funzione che implementa il programma principale
    '''

    if len(sys.argv) < 5 :
        print ('usage: ', sys.argv[0], 'xMin xMax numero seed')
        exit ()

    xMin = float (sys.argv[1])  # minimum of the histogram drawing range
    xMax = float (sys.argv[2])  # maximum of the histogram drawing range
    seed = float (sys.argv[4])
    N    = int (sys.argv[3])

    print (' -------- ')
    print (' minimum : ', xMin)
    print (' maximum : ', xMax)
    print (' seed    : ', seed)
    print (' N       : ', N)
    print (' -------- ')

    randlist = generate_range (xMin, xMax, N, seed)

    # plotting of the generated list of numbers in a histogram

    nBins = floor (len (randlist) / 20.) + 1     # number of bins of the hitogram
    bin_edges = np.linspace (xMin, xMax, nBins)  # edges o the histogram bins

    # disegno della funzione
    fig, ax = plt.subplots ()
    ax.set_title ('Histogram of random numbers', size=14)
    ax.set_xlabel ('random value')
    ax.set_ylabel ('events in bin')
    ax.hist (randlist,      # list of numbers
             bins = bin_edges,
             color = 'orange',
             # normed = True,
            )
    plt.savefig ('es_4.4.png')


# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 


if __name__ == "__main__":
    main ()

4.8.6. Exercise 4.6

• ex_4.6.py

#!/usr/bin/python
'''
Si implementi un generatore di numeri pseduo-casuali che utilizzi il metodo della funzione inversa
per generare eventi distribuiti secondo distribuzione di probabilita' esponenziale.
  * Si utilizzi la libreria MatPlotLib per visualizzare la distribuzione dei numeri generati.
python3 es_4.6.py 100 2.4
'''

import sys
import random
import matplotlib.pyplot as plt
import numpy as np
from math import floor

from myrand import generate_TAC


def inv_exp (y, lamb = 1) :
    '''
    Inverse of the primitive of the exponential PDF.
    pdf(x) = lambda * exp(-lambda x) x >= 0, 0 otherwise.
    F(x) = int_{0}^{x} pdf(x)dx = 1 - exp(-lambda * x) for x >= 0, 0 otherwise.
    F^{-1}(y) = - (ln(1-y)) / lambda
    '''
    return -1 * np.log (1-y) / lamb


# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 


def main () :
    '''
    Funzione che implementa il programma principale
    '''

    if len(sys.argv) < 3 :
        print ('usage: ', sys.argv[0], 'N tau [seed=10]')
        exit ()

    tau  = float (sys.argv[2])
    if tau <= 0 :
        print ('The tau parameter has to be positive')
        exit ()
    lamb = 1./tau ;    
    N    = int (sys.argv[1])
    seed = 10
    if len(sys.argv) >= 4 : 
        seed = float (sys.argv[3])

    random.seed (seed)
    randlist = []
    for i in range (N):
        randlist.append (inv_exp (random.random (), lamb))

    # plotting of the generated list of numbers in a histogram

    nBins = floor (len (randlist) / 100.)        # number of bins of the hitogram
    bin_edges = np.linspace (0., 3., nBins + 1)  # edges o the histogram bins

    # disegno della funzione
    fig, ax = plt.subplots ()
    ax.set_title ('Histogram of random numbers', size=14)
    ax.set_xlabel ('random value')
    ax.set_ylabel ('events in bin')
    ax.hist (randlist,      # list of numbers
             bins = bin_edges,
             color = 'orange',
             # normed = True,
            )

    plt.savefig ('es_4.6.png')


# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 


if __name__ == "__main__":
    main ()

4.8.7. Exercise 4.7

• stats.py

#!/usr/bin/python

from math import sqrt


class stats :
    '''calculator for statistics of a list of numbers'''

    summ = 0.
    sumSq = 0.
    N = 0
    sample = []

    def __init__ (self, sample):
        '''
        reads as input the collection of events,
        which needs to be a list of numbers
        '''
        self.sample = sample
        self.summ = sum (self.sample)
        self.sumSq = sum ([x*x for x in self.sample])
        self.N = len (self.sample)

    # ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 

    def mean (self) :
        '''
        calculates the mean of the sample present in the object
        '''
        return self.summ / self.N

    # ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 

    def variance (self, bessel = True) :
        '''
        calculates the variance of the sample present in the object
        '''
        var = self.sumSq / self.N - self.mean () * self.mean ()
        if bessel : var = self.N * var / (self.N - 1)
        return var

    # ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 

    def sigma (self, bessel = True) :
        '''
        calculates the sigma of the sample present in the object
        '''
        return sqrt (self.variance (bessel))

    # ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 

    def sigma_mean (self, bessel = True) :
        '''
        calculates the sigma of the sample present in the object
        '''
        return sqrt (self.variance (bessel) / self.N)

    # ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 

    def append (self, x):
        '''
        add an element to the sample
        '''
        self.sample.append (x)
        self.summ = self.summ + x
        self.sumSq = self.sumSq + x * x
        self.N = self.N + 1

4.8.8. Exercise 4.8

• ex_4.8.py

#!/usr/bin/python
'''
Partendo dal lavoro svolto durante la Lezione 3,
si implementi un oggetto chiamato ```stats```, 
che calcola le statistiche associate ad un campione di numeri
salvato in una lista di Python.
  * quali diverse opzioni di progettazione sono possibili per questo oggetto?
  * che variabili è necessario aggiungere alla classe per garantirne la funzionalità?
  * che valore devono avere queste variabili in fase di inizilizzazione?
Si collaudi l'oggetto ```stats``` con ciascuono degli agoritmi di generazione implementati.
In particolare, poi:
  * si verifichi che il valore della varianza della distribuzione uniforme corrisponde alle attese
    (quale è l'incertezza associata al numero ottenuto?)
  * si verifichi che il valore della varianza ottenuta con la tecnica del teorema centrale del limite
    corrisponda a quello atteso

python3 es_4.8.py 0. 3. 20000 2.4 10
'''

import sys
import random
import matplotlib.pyplot as plt
import numpy as np
from math import floor

from myrand import generate_TCL_ms
from stats import stats


def main () :
    '''
    Funzione che implementa il programma principale
    '''

    if len(sys.argv) < 6:
        print ('usage: ', sys.argv[0], 'mean sigma numero seed numero_somme')
        exit ()

    mean  = float (sys.argv[1])  # minimum of the histogram drawing range
    sigma = float (sys.argv[2])  # maximum of the histogram drawing range
    seed  = float (sys.argv[4])
    N     = int (sys.argv[3])
    N_sum = int (sys.argv[5])

    print (' -------- ')
    print (' mean  : ', mean)
    print (' sigma : ', sigma)
    print (' seed  : ', seed)
    print (' N     : ', N)
    print (' N_sum : ', N_sum)
    print (' -------- ')

    randlist = generate_TCL_ms (mean, sigma, N, N_sum, seed)

    my_stats = stats (randlist)
    my_mean = my_stats.mean ()
    my_sigma = my_stats.sigma ()
    print ('media:    ', my_mean) ;
    print ('varianza: ', my_stats.variance ()) ;
    print ('sigma:    ', my_sigma) ;

    # plotting of the generated list of numbers in a histogram

    nBins = floor (len (randlist) / 400.)             # number of bins of the hitogram
    bin_edges = np.linspace (mean - 3 * sigma, mean + 3 * sigma, nBins + 1)  # edges o the histogram bins

    # disegno della funzione
    fig, ax = plt.subplots ()
    ax.set_title ('Histogram of random numbers', size=14)
    ax.set_xlabel ('random value')
    ax.set_ylabel ('events in bin')
    ax.hist (randlist,      # list of numbers
             bins = bin_edges,
             color = 'orange',
             # normed = True,
            )
    yrange = plt.ylim ()
    ax.plot ([my_mean, my_mean], yrange, color = 'red') 
    ax.plot ([my_mean - my_sigma, my_mean - my_sigma], yrange, linestyle='dashed', color = 'red') 
    ax.plot ([my_mean + my_sigma, my_mean + my_sigma], yrange, linestyle='dashed', color = 'red') 

    plt.savefig ('es_4.8.png')


# ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 


if __name__ == "__main__":
    main ()