3.6. Examples: Lecture 3

3.6.1. Drawing a histogram

• drawing_03.py

#!/usr/bin/python
'''
example on the use of matplotlib.pyplot
to draw a histogram with automatic binning setup
'''

import matplotlib.pyplot as plt
import numpy as np


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

    # the sample to be represented as an histogram
    sample = [1.3977702773941403, 0.74021968956204, 0.559129026598526, 0.034722468927208275, 1.472659179569878, 1.2276631640103624, 1.0483359119412394, 0.8904343375674812, 0.14771386083901775, 1.8655990499580688, 1.3776373531042858, 1.3348079223655975, -0.6326424273210922, -0.612527851941395, 1.594263621951928, 0.510681111246625, 0.7427516727970715, 1.6250978742928908, 0.6773261829611841, 1.4922122786883083]

    fig, ax = plt.subplots (nrows = 1, ncols = 1)
    ax.hist (sample,
             color = 'orange',
            )
    ax.set_title ('Histogram example', size=14)
    ax.set_xlabel ('variable')
    ax.set_ylabel ('event counts per bin')

    plt.savefig ('drawing_03.png')
    plt.show ()
    

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

    
if __name__ == "__main__":
    main ()

3.6.2. Drawing a histogram with predefined binning

• drawing_04.py

#!/usr/bin/python
'''
example on the use of matplotlib.pyplot
with predefined binning
'''

import matplotlib.pyplot as plt
import numpy as np


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

    # the sample to be represented as an histogram
    sample = [1.3977702773941403, 0.74021968956204, 0.559129026598526, 0.034722468927208275, 1.472659179569878, 1.2276631640103624, 1.0483359119412394, 0.8904343375674812, 0.14771386083901775, 1.8655990499580688, 1.3776373531042858, 1.3348079223655975, -0.6326424273210922, -0.612527851941395, 1.594263621951928, 0.510681111246625, 0.7427516727970715, 1.6250978742928908, 0.6773261829611841, 1.4922122786883083]

    N_bins = 10
    xMin = 0.
    xMax = 1. 
    bin_edges = np.linspace (xMin, xMax, N_bins)
    print ('length of the bin_edges container:', len (bin_edges))

    fig, ax = plt.subplots (nrows = 1, ncols = 1)
    ax.hist (sample,
             bins = bin_edges,
             color = 'orange',
            )
    ax.set_title ('Histogram example', size=14)
    ax.set_xlabel ('variable')
    ax.set_ylabel ('event counts per bin')

    plt.savefig ('drawing_04.png')
    plt.show ()
    

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

    
if __name__ == "__main__":
    main ()

3.6.3. Drawing a histogram and choice of bin number

• drawing_05.py

#!/usr/bin/python
'''
example on the use of matplotlib.pyplot
with predefined binning and a reasonable choice
of the bin number
'''

import matplotlib.pyplot as plt
import numpy as np
from math import ceil


def sturges (N_events) :
    return ceil (1 + 3.322 * np.log (N_events))
    

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

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

    # the sample to be represented as an histogram
    with open('dati.txt') as f:
        sample = [float (x) for x in f.readlines()]

    xMin = 0.
    xMax = 1. 
    N_bins = sturges (len (sample))
    bin_edges = np.linspace (xMin, xMax, N_bins)
    fig, ax = plt.subplots (nrows = 1, ncols = 1)
    ax.hist (sample,
             bins = bin_edges,
             color = 'orange',
            )
    ax.set_title ('Histogram example', size=14)
    ax.set_xlabel ('variable')
    ax.set_ylabel ('event counts per bin')

    plt.savefig ('drawing_05.png')
    plt.show ()
    

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

    
if __name__ == "__main__":
    main ()

3.6.4. Comparing different binning choices

• drawing_06.py

#!/usr/bin/python
'''
example on the use of matplotlib.pyplot
with the comparison among different binning choices
'''

import matplotlib.pyplot as plt
import numpy as np
from math import ceil, floor


def sturges (N_events) :
    return ceil (1 + np.log2 (N_events))
    

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

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

    # the sample to be represented as an histogram
    with open('dati.txt') as f:
        sample = [float (x) for x in f.readlines()]

    xMin = floor (min (sample))
    xMax = ceil (max (sample))
    N_bins = sturges (len (sample))
    fig, axes = plt.subplots (nrows = 3, ncols = 1)
    bin_edges = np.linspace (xMin, xMax, N_bins)
    axes[0].hist (sample,
             bins = bin_edges,
             color = 'orange',
             label = 'Sturges',
            )
    bin_edges = np.linspace (xMin, xMax, N_bins * 10)
    axes[1].hist (sample,
             bins = bin_edges,
             color = 'blue',
             label = 'Sturges x 10',
            )
    bin_edges = np.linspace (xMin, xMax, int (N_bins / 10))
    axes[2].hist (sample,
             bins = bin_edges,
             color = 'red',
             label = 'Sturges / 10',
            )
    
    for ax in axes :
        ax.set_xlabel ('variable')
        ax.set_ylabel ('event counts')
        ax.legend ()

    plt.savefig ('drawing_06.png')
    plt.show ()
    

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

    
if __name__ == "__main__":
    main ()

3.6.5. Testing logarithmic scales

• drawing_07.py

#!/usr/bin/python
'''
example on the use of matplotlib.pyplot
and the use of the logarithmic scale on the vertical axis
'''

import matplotlib.pyplot as plt
import numpy as np
from math import sqrt


def Gaussian (x, mean, sigma) :
    if sigma == 0 : return float (x == 0.)
    return np.exp (-0.5 * pow ((x - mean) / sigma, 2.)) / (sqrt (2 * np.pi) * sigma)
    

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

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

    # preparing the set of points to be drawn 
    x_coord = np.linspace (-3., 3., 10000)
    y_coord = np.arange (0., x_coord.size)
    for i in range (x_coord.size):
        y_coord[i] = Gaussian (x_coord[i], 0., 1.)

    fig, axes = plt.subplots (nrows = 1, ncols = 2)

    axes[0].plot (x_coord, y_coord)
    axes[0].set_title ('linear y scale', size=14)
    axes[0].set_xlabel ('x')
    axes[0].set_ylabel ('y')

    axes[1].plot (x_coord, y_coord)
    axes[1].set_title ('log y scale', size=14)
    axes[1].set_xlabel ('x')
    axes[1].set_yscale ('log')

    plt.savefig ('drawing_07.png')
    plt.show ()
    

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

    
if __name__ == "__main__":
    main ()

3.6.6. Writing data into text files

• saving.py

#!/usr/bin/python
'''
example on saving information on txt files
'''

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

    # the sample to be represented as an histogram
    sample = [1.3977702773941403, 0.74021968956204, 0.559129026598526, 0.034722468927208275, 1.472659179569878, 1.2276631640103624, 1.0483359119412394, 0.8904343375674812, 0.14771386083901775, 1.8655990499580688, 1.3776373531042858, 1.3348079223655975, -0.6326424273210922, -0.612527851941395, 1.594263621951928, 0.510681111246625, 0.7427516727970715, 1.6250978742928908, 0.6773261829611841, 1.4922122786883083]

    with open ('sample.txt', 'w') as fp :
        for item in sample:
            # write each item on a new line
            fp.write (str (item) + '\n')
        print ('Done')    


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

    
if __name__ == "__main__":
    main ()

3.6.7. Reading data from text files

• reading.py

#!/usr/bin/python
'''
example on saving information on txt files
'''

import sys


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

    if (len (sys.argv) < 2) :
        print ('usage:', sys.argv[0], 'inputfile.txt')
        sys.exit ()

    # read the file
    with open(sys.argv[1]) as f:
        sample = [float (x) for x in f.readlines()]

    print (len (sample))    


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

    
if __name__ == "__main__":
    main ()

3.6.8. Calculating sample moments

• stats_01.py

#!/usr/bin/python
'''
example on saving information on txt files
'''

import sys
from stats import mean, variance

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

    if (len (sys.argv) < 2) :
        print ('usage:', sys.argv[0], 'inputfile.txt')
        sys.exit ()

    # read the file
    with open(sys.argv[1]) as f:
        sample = [float (x) for x in f.readlines()]

    print (len (sample))    

    print ('sample mean:', mean (sample)) ;
    print ('sample variance:', variance (sample)) ; 


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

    
if __name__ == "__main__":
    main ()

• stats.py

from math import sqrt 

def mean (sample) :
    '''
    calculates the mean of the sample present in the object
    '''
    summ = sum (sample)
    N = len (sample)
    return summ / N

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

def variance (sample, bessel = True) :
    '''
    calculates the variance of the sample present in the object
    '''
    summ = 0.
    sum_sq = 0.
    N = len (sample)
    for elem in sample :
       summ += elem
       sum_sq += elem * elem
    var = sum_sq / N - summ * summ / (N * N)
    if bessel : var = N * var / (N - 1)
    return var

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

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

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

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

3.6.9. Visualising sample moments

• stats_02.py

#!/usr/bin/python
'''
example on saving information on txt files
'''

import sys
from math import floor, ceil
from stats import mean, sigma
import numpy as np
import matplotlib.pyplot as plt


def sturges (N_events) :
    return ceil (1 + np.log2 (N_events))
    

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

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

    if (len (sys.argv) < 2) :
        print ('usage:', sys.argv[0], 'inputfile.txt')
        sys.exit ()

    # read the file
    with open(sys.argv[1]) as f:
        sample = [float (x) for x in f.readlines()]

    sample_mean = mean (sample) ;
    sample_sigma = sigma (sample) ;

    xMin = floor (min (sample))
    xMax = ceil (max (sample))
    N_bins = sturges (len (sample))
    bin_edges = np.linspace (xMin, xMax, N_bins)
    fig, ax = plt.subplots (nrows = 1, ncols = 1)
    ax.hist (sample,
             bins = bin_edges,
             color = 'orange',
            )
    ax.set_title ('Histogram example', size=14)
    ax.set_xlabel ('variable')
    ax.set_ylabel ('event counts per bin')

    vertical_limits = ax.get_ylim ()
    ax.plot ([sample_mean, sample_mean], vertical_limits, color = 'blue')
    ax.plot ([sample_mean - sample_sigma, sample_mean - sample_sigma], 
             vertical_limits, color = 'blue', linestyle = 'dashed')
    ax.plot ([sample_mean + sample_sigma, sample_mean + sample_sigma], 
             vertical_limits, color = 'blue', linestyle = 'dashed')

    plt.savefig ('stats_02.png')
    plt.show ()



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

    
if __name__ == "__main__":
    main ()

3.6.10. Calculating values of a probability density function

• pdf_01.py

#!/usr/bin/python
'''
example on the use of scipy
'''

from scipy.stats import norm


def main () :
    '''
    Funzione che implementa il programma principale
    '''
    
    mean = 1.
    sigma = 0.5
    print ('the maximum value of the Gaussian distribution with mean = ' +
           str (mean) + 
           ' and sigma = ' +
           str (sigma) + 
           ' is : ' +
           str (norm.pdf (mean, mean, sigma))
           )

    norm_fix = norm (mean, sigma)
    print (norm_fix.pdf (mean))

    x = mean + sigma / 2.
    print (norm.pdf (x, mean, sigma))

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

    
if __name__ == "__main__":
    main ()

3.6.11. Calculating the first four moments of a pdf

• pdf_02.py

#!/usr/bin/python
'''
example on the use of scipy
'''

from scipy.stats import norm


def main () :
    '''
    Funzione che implementa il programma principale
    '''
    
    mean = 1.
    sigma = 0.5
    norm_fix = norm (mean, sigma)

    ave, var, skew, kurt = norm_fix.stats (moments='mvsk')
    print (ave, var, skew, kurt)



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

    
if __name__ == "__main__":
    main ()

3.6.12. Calculating values of a cumulative distribution function

• cdf_01.py

#!/usr/bin/python
'''
example on the use of scipy
'''

from scipy.stats import norm


def main () :
    '''
    Funzione che implementa il programma principale
    '''
    
    mean = 1.
    sigma = 0.5
    print ('the value of the Gaussian distribution cumulative at its mean is: ' +
           str (norm.cdf (mean, mean, sigma))
           )

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

    
if __name__ == "__main__":
    main ()

3.6.13. Numerical integration

• integrate.py

#!/usr/bin/python
'''
example on the use of scipy.integrate
'''

from scipy.integrate import quad
import numpy as np
from math import exp

def expon (x) :
    return exp (-1 * x)
    

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

    
def main () :
    '''
    Funzione che implementa il programma principale
    '''
    # definition of a polinomial function
    polin = lambda x : x**2 + x + 1
    area = quad (polin, 0., 4.)
    print ('area = ', area[0])
    print ('absolute error estimate = ', area[1])

    area = quad (expon, 0, np.inf)
    print ('area = ', area[0])
    print ('absolute error estimate = ', area[1])
    

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

    
if __name__ == "__main__":
    main ()