Python-lezioni

Introduzione a markdown

Gli elenchi numerati sono fatti con un trattino: 1. prima voce 1. seconda voce 1. terza voce

I caratteri si possono gestire nle modo seguente: - corsivo oppure corsivo - grassetto oppure grassetto - grassetto in corsivo oppure grassetto - sottolineato: testo - monospazio (con l’apostrofo invero - backquote - fatto con ALT+GR+’)

Si può creare una linea orizzontale

---

Una citazione si può indicare con un rientro:

Sempre amai le rose che non colsi, sempre amai le cose che potevano essere non sono state

Si possono rappresentare tabelle nel modo seguente

Colonna 1	Colonna 2	Colonna 3
A	B	C
D	E	F

Le equazioni si possono scrivere mediante LaTeX. Se in linea si ua il “dollaro”: \(\int_a^b f(x) dx\) se, invece, sono al centro della riga, si usa il doppio dollar \[\mu=\frac{1}{n}\sum_{i=1}^{n}x_i.\]

Le prime operazioni

Con il software di base si può fare ben copo su python

2+5

## 7

4/3

## 1.3333333333333333

L’elevamento a potenza si calcola scrivendo due volte il simbolo del prodotto:

2**3

## 8

Vale l’elevamento a potenza con eponente reale (e con base negativa)

(-2)**(1/2)

## (8.659560562354934e-17+1.4142135623730951j)

Se vogliamo calcolare un logaritmo e diamo il comando

log(2)

riceviamo un errore perché Python non è in grado di fare questo calcolo serve un pacchetto apposito per fare matematica “seria”.

Il primo pacchetto

Quando si carica un pacchetto si utilizza la sintassi

import “nome pacchetto” as “nome a scelta”

Per effettuare i calcoli atematici serve il pacchetto numpy (numerical python)

import numpy as m
m.log(2)

## 0.6931471805599453

Notiamo che se usiamo la funzione “radice quadrata” del pacchetto numpy, perdiamo la possibilità di ottenere numeri immaginali/complessi:

m.sqrt(-2)

## nan
## 
## <string>:1: RuntimeWarning: invalid value encountered in sqrt

È possibile non dare alcun nome a un pacchetto. Si ottiene il risultato mettendo un asterisco al posto del nome e usando la seguente sintassi:

from numpy import *
log(2)

## 0.6931471805599453

Questa operazione non è consigliabile perché alcuni pacchetti hanno, al loro interno, funzion con lo stesso nome e non dando nomi diversi ai pacchetti, queste funzioni si sovrapporrebbero (cancellando la più vecchia). Usiamo, quindi, il nome standard utilizzato in Python che si trova nei suggerimenti di Colab:

import numpy as np
np.log(2)

## 0.6931471805599453

Lavorare con le matrici

Le matrici si possono inserire, all’interno del pacchetto numpy in modi diversi. Un metodo è mediante il coamndo matrix con la seguente sintassi

numpy.matrix([[riga1],[riga2],…,[rigan]])

dove ogni riga è un insieme di elementi divisi da una virgola:

A=np.matrix([[10,-4,5],[7,-2,1]])
A

## matrix([[10, -4,  5],
##         [ 7, -2,  1]])

Le matrici si possono inserire anche usando un’altra sintassi. Si inserisce una stringa in cui gli elementi sono divisi da una virgola per indicare le righe e da un punto e virgola per indicare le colonne.

A=np.matrix('10,-4,5;7,-2,1')
A

## matrix([[10, -4,  5],
##         [ 7, -2,  1]])

Per richiamare gli elementi di una matrice si mettono le “coordinate” all’interno di parentesi quadrate.

A[1,1]

## -2

Notiamo che l’elemento [1,1] non è quello che ci aspettiamo. La numerazione in Python, infatti, parte dallo zero. L’elemento della prima riga e della prima colonna, quindi, si ottiene nel modo seguente:

A[0,0]

## 10

Per avere tutti gli elementi di una riga o di una colonna si usano i due punti:

A[:,1]

## matrix([[-4],
##         [-2]])

Quando si inserisce un elemento di inizio e uno di fine, quello finale non è mai compreso

B=A[0:2,0:2]
B

## matrix([[10, -4],
##         [ 7, -2]])

Gli ultimi elementi si indicano inserendo come elmento della matrice un numero negativo

A[:,-1] # Ultima colonna della matrice A

## matrix([[5],
##         [1]])

A[:,-2] # Penultima colonna della matrice A

## matrix([[-4],
##         [-2]])

A[:,-3:-1] # Ultima e penultima colonna della matrice A (come sempre, non si comprende l'emento finale)

## matrix([[10, -4],
##         [ 7, -2]])

La matrice trasposta si indica così:

A.T

## matrix([[10,  7],
##         [-4, -2],
##         [ 5,  1]])

La matrice inversa si indica così:

B.I

## matrix([[-0.25 ,  0.5  ],
##         [-0.875,  1.25 ]])

Verifichiamo che la matrice inversa rispetti la proprietà per cui, moltipliata per se stessa non inversa, dia la matrice identità:

B.I*B

## matrix([[ 1.00000000e+00,  0.00000000e+00],
##         [-1.77635684e-15,  1.00000000e+00]])

Il determinante è:

np.linalg.det(B)

## 8.000000000000005

Il prodotto matriciale “classico” è

B*A

## matrix([[ 72, -32,  46],
##         [ 56, -24,  33]])

Avendo due vettori colonna \(x\) e \(y\) (di dimensione \(n\)), il loro prodotto può essere effettuato in forma matriciale moltipliando una riga per una colonna: \[x^\intercal y = \sum_{i=1}^n x_i y_i\] oppure in forma di prodotto “elemento per elemento” moltiplicando ogni elementi di \(x\) per ogni elemento di \(y\): \[\left[\begin{array}{c} x_{1}y_{1}\\ x_{2}y_{2}\\ ...\\ x_{n}y_{n} \end{array}\right]\]

Vediamo un esempio:

# Definizione di due vettori
x=np.matrix([[1],[3],[5]])
y=np.matrix([[-1],[0],[4]])
# Prodotto "interno" (inner product)
x.T*y

# Prodotto elemento per elemento

## matrix([[19]])

np.multiply(x,y)

## matrix([[-1],
##         [ 0],
##         [20]])

Creazione di numeri casuali

X=np.random.normal(0,1,10000)
print(np.mean(X))

## 0.0002145895411984796

print(np.std(X))

## 1.007373224363633

print(np.quantile(X,[0,0.5,1]))

## [-3.75088260e+00 -1.90756056e-03  4.27552805e+00]

print(np.quantile(X,[0,0.5,1])[0])

## -3.750882596420622

I numeri casuali differiscono ogni volta che li si crea

print(np.random.random(4))

## [0.0269277  0.53129832 0.56479739 0.0781205 ]

print(np.random.random(4))

## [0.1762087  0.88405007 0.50152587 0.86183923]

Se voglio fisare la “radice” dei numeri casuali uso il comando “seed” all’interno del sottopacchetto random

np.random.seed(1)
print(np.random.random(4))

## [4.17022005e-01 7.20324493e-01 1.14374817e-04 3.02332573e-01]

np.random.seed(1)
print(np.random.random(4))

## [4.17022005e-01 7.20324493e-01 1.14374817e-04 3.02332573e-01]

Creare funzioni

Una funzione inizia con il comando def e prosegue con il ‘nome della funzione’ e con gli argomenti tra parentesi; la linea che inizia la funzione temrina con i due punti. I comandi devono essere indicato con un rientro.

def nome(arg1,arg2):

comandi

I comandi devono essere tutti preceduti da un rientro. La prima riga senza rientro è la prima riga da intendersi al di fuori della funzione definita.

Creiamo una funzione per lanciare \(n\) dadi ciacuno dei quali abbia \(f\) facce e ne mostriamo i valori

def dadi(n,f):
  return np.random.random_integers(1,f,n)

Il comando return serve per indicare a Python quale valore deve “restituire” come output della funzione

dadi(2,6)

## array([1, 1])
## 
## <string>:2: DeprecationWarning: This function is deprecated. Please call randint(1, 6 + 1) instead

Si possono inserire dei valori di default per gli input in modo che la funzione restituisca comunque dei valori anche se non si insericono delle variabili in entrata:

def dadi(n=2,f=6):
  return np.random.random_integers(1,f,n)

dadi()

## array([2, 5])

Si possono inserire alcune variabili con un valore predefinito e altre senza. Tuttavia è necessario che le variabili con default seguano le altre:

def dadi(n,f=6):
  return np.random.random_integers(1,f,n)

dadi(2)

## array([6, 5])

Si possono inserire anche le variabili in ordine diverso da come sono state definite nella funzione e, tuttavia, in questo caso occorre utilizzare i nomi delle variabili:

dadi(f=6,n=3)

## array([2, 3, 5])

Si possono restituire più output della funzione inserendoli dopo il comando return e separati da una virgola. Se, per esempio, vogliamo vedere sia i valori dei dadi sia la loro somma, possiamo scrivere:

def dadi(n=1,f=6):
  D=np.random.random_integers(1,f,n)
  S=np.sum(D)
  return D, S

dadi(2,8)

## (array([7, 6]), 13)
## 
## <string>:2: DeprecationWarning: This function is deprecated. Please call randint(1, 8 + 1) instead

Si possono attribuire gli output della funzione a due variabili “esterne” inserendole prima della funzione, separate da una virgola:

facce,somma=dadi(3,8)
print(facce)

## [3 5 4]

print(somma)

## 12

Si può richiamare solo uno degli output utilizzando le parentesi quadre e ricordando che il primo elemento di un inieme è indicizzato con il numero zero:

print(dadi(3,8)[0])

## [5 3 5]

print(dadi(3,8)[1])

## 16

Il ciclo for

Quando si vuole ripetere lo stesso comando per più volte, camabiando ogni volta ua parte del comando, si deve utilizzare un cosiddetto “ciclo for” che funziona nel modo seguente:

for contatore in range: comandi

Tutti i comandi devono essere spostati di un rientro, “contatore” è una variabile che serve per tenere il conto del ciclo e “range” è un insieme a cui appartiene la variabile “contatore”. Vediao un esempio.

for i in [1,2,3]:
  print(i**2)

## 1
## 4
## 9

Vediamo come creare una funzione che calcoli un percorso randomico \(w_t\) (random walk) formato in questo modo: 1. il valore iniziale è nullo \(w_{t_0}=0\); 1. in ogni istante \(t\) il valore \(w_t\) è dato da \(w_{t-1}\) a cui si somma 1 con probabilità 0.5 o si sottrae 1 con la tea probabilità.

#Definisco il vettore w
n=10 # numero di simulazioni
w=np.zeros((n,1))
# Valore iniziale
print(w)

## [[0.]
##  [0.]
##  [0.]
##  [0.]
##  [0.]
##  [0.]
##  [0.]
##  [0.]
##  [0.]
##  [0.]]

for i in range(1,n):
  w[i]=w[i-1]+2*np.random.binomial(1,0.5)-1
w

## array([[ 0.],
##        [ 1.],
##        [ 0.],
##        [-1.],
##        [ 0.],
##        [ 1.],
##        [ 2.],
##        [ 3.],
##        [ 2.],
##        [ 3.]])

Creiamo una funzione per definire il random walk inserendo come input il numero di osservazioni \(n\)

def RW(n):
  w=np.zeros((n,1))
  w[0]=0
  for i in range(1,n):
    w[i]=w[i-1]+2*np.random.binomial(1,0.5)-1
  return w

RW(10)

## array([[0.],
##        [1.],
##        [2.],
##        [3.],
##        [2.],
##        [1.],
##        [2.],
##        [3.],
##        [4.],
##        [3.]])

Ora simuliamo, tutti insieme, \(N\) percorsi (ognuno di \(n\) passi)

def RW(n,N):
  w=np.zeros((n,N))
  w[0,:]=0
  for i in range(1,n):
    w[i,:]=w[i-1,:]+2*np.random.binomial(1,0.5,N)-1
  return w

RW(10,2)

## array([[ 0.,  0.],
##        [-1., -1.],
##        [-2.,  0.],
##        [-3., -1.],
##        [-2.,  0.],
##        [-3., -1.],
##        [-4., -2.],
##        [-5., -3.],
##        [-4., -4.],
##        [-3., -5.]])

Posso applicare alcune funzioni statistiche, come la media e la varianza, a una matrice indicando rispetto a quale dimension (axis) devo effettuare il calcolo

np.mean(RW(1000,3),axis=0) #si calcola la media dei valori sulle colonne

## array([ 2.436, -3.512, -5.908])

np.mean(RW(10,3),axis=1) #si calcola la media dei valori sulle righe

## array([ 0.        ,  1.        ,  0.66666667,  0.33333333,  0.66666667,
##         1.        ,  0.        ,  0.33333333,  0.        , -1.        ])

np.std(RW(1000,3),axis=0)

## array([9.14450195, 6.81435074, 8.50387559])

np.quantile(RW(1000,3),[0,0.25,0.5,0.74,1],axis=0)

## array([[-24. ,  -3. , -12. ],
##        [-14. ,  18. ,   4. ],
##        [-10. ,  27. ,   9.5],
##        [ -6. ,  34. ,  15. ],
##        [ 13. ,  58. ,  31. ]])

Si può applicare una funzione statistica (come la media, o la varianza) anche in un altro modo. Si fa seguire una matrice da un punto e dal nome della funzione. Il comando termina inserendo tra parentesi tonde la dimenzione (righe, \(0\), o colonne \(1\)) rispetto alle quali si vuole calcolare la funzione.

RW(1000,3).mean(0)

## array([10.562, 21.14 , -1.108])

RW(1000,3).std(0)

## array([18.48303427, 11.73523089,  6.66756897])

Creare figure

Per creare figure abbiamo bisogno, ovviamente, di un pacchetto! Si chiama matplotlib e, in particolare, ci serve, al suo interno, il pacchetto pyplot.

import matplotlib.pyplot as plt

Il comando plot prende in entrata, come elemento essenziale, solo il valore delle ordinate

plt.plot(RW(100,1))

Inseriamo valori a nostra scelta sull’asse delle ascisse

plt.plot(np.arange(0,1,1/250),RW(250,1))

Inseriamo le etichette sull’asse x e sull’asse y

plt.plot(np.arange(0,1,1/250),RW(250,1))
plt.xlabel('year')
plt.ylabel('Random Walk')

Quando si crea una figura si possono isolare tutti i comandi all’interno di due linee: - plt.figure() - plt.show()

Con questi comandi si comunica a Python che stiamo creando una figura e che, alla fine, deve mostrarla.

Questi comandi permettono a Python di creare la figura “tutta insieme” e gestendo nel modo migliore le dimensioni degli assi e la posizione degli elementi del grafico.

plt.figure()
plt.plot(np.arange(0,1,1/250),RW(250,3))
plt.xlabel('year')
plt.ylabel('RW')
plt.title('Random Walk')

# Inseriamo una legenda dando un nome ai percorsi 
plt.legend(['RW1','RW2','RW3'])
# La posizione della legenda è scelta autonomamente da Python nella zona del grafico più "libera"
plt.grid()
plt.show()

Voglio creare 100 simulazioni di percorsi giornalieri e indicare sul grafico: - tutti i percori in grigio chiaro - la linea del percoro medio in nero e tratteggiata - il quinto percentile in rosso - il novantacinquesimo percentile in verde - la legenda

X=RW(250,100)
asse=np.arange(0,1,1/250)
plt.figure()
plt.plot(asse,X,color='lightgray') # Scelgo il colore delle linee (le voglio tutte grige)
plt.xlabel('year')
plt.ylabel('RW')
plt.title('Random Walk')
plt.grid()
plt.plot(asse,np.mean(X,1),color='black',linestyle='dashed',label='Media') # La media è nera e tratteggiata
plt.plot(asse,np.quantile(X,0.05,1),color='red',label='Quantile 0.05')
plt.plot(asse,np.quantile(X,0.95,1),color='green',label='Quantile 0.95')
plt.legend()
plt.show()

La legenda si può posizionare al di fuori del grafico mediante un’opzione speciale del comando legend

plt.legend(bbox_to_anchor=(x,y),loc=‘center left’)

Le variabili \(x\) e \(y\) sono espresse in un’unità di miura pari alla “lunghezza del grafico”. Per esempio, se si sceglie

bbox_to_anchor=(1,0.5)

si intende posizionare la legenda a distanza di un grafico dal bordo di sinistra e al distanza di mezzo grafico dal bordo inferiore

Se si inserisce una delle variabili con numeri negativi, si ottiene il seguente risultato: - \(x<0\) si posiziona la legenda a sinistra del grafico - \(y<0\) si posiziona la legenda sotto il grafico

L’opzione loc indica dove posizionare il riferimento per la legenda. Il contenuto “centro sinistra” indica che il grafico si deve intendere agganciato alla sinistra della legenda.

plt.figure()
plt.plot(asse,X,color='lightgray') # Scelgo il colore delle linee (le voglio tutte grige)
plt.xlabel('year')
plt.ylabel('RW')
plt.title('Random Walk')
plt.grid()
plt.plot(asse,np.mean(X,1),color='black',linestyle='dashed',label='Media') # La media è nera e tratteggiata
plt.plot(asse,np.quantile(X,0.05,1),color='red',label='Quantile 0.05')
plt.plot(asse,np.quantile(X,0.95,1),color='green',label='Quantile 0.95')
plt.legend(bbox_to_anchor=(1,0.5),loc='center left')
plt.show()

Vogliamo creare un grafico con le seguenti caratteristiche - rappresentiamo solo i quantili 0.05 e 0.95 - rappresentiamo la media (nera tratteggiata) - coloriamo lo spazio tra i quantili

X=RW(250,100)
asse=np.arange(0,1,1/250)
plt.figure()
plt.xlabel('year')
plt.ylabel('RW')
plt.title('Random Walk')
plt.grid()
plt.plot(asse,np.mean(X,1),color='black',label='Media')
plt.plot(asse,np.quantile(X,0.05,1),color='red',label='Quantile 0.05',linestyle='dotted')
plt.plot(asse,np.quantile(X,0.95,1),color='green',label='Quantile 0.95',linestyle='dotted')
# Quando si colora lo spazio, si possno specificare sia il colore sia la "trasparenza" del colore (con il paramtero "alpha" tra 0 e 1)
plt.fill_between(asse,np.quantile(X,0.95,1),np.quantile(X,0.05,1),alpha=0.3,color='blue')
plt.legend()
plt.show()

In quest’altro esempio vogliamo colorare di verde (trasparente) l’area tra il quantile 0.95 e la media, mentre vogliamo colorare di rosso (trasparente) l’area tra la media e il quantile 0.05

X=RW(250,100)
asse=np.arange(0,1,1/250)
plt.figure()
plt.xlabel('year')
plt.ylabel('RW')
plt.title('Random Walk')
plt.grid()
plt.plot(asse,np.mean(X,1),color='black',label='Media')
plt.plot(asse,np.quantile(X,0.05,1),color='red',label='Quantile 0.05',linestyle='dotted')
plt.plot(asse,np.quantile(X,0.95,1),color='green',label='Quantile 0.95',linestyle='dotted')
plt.fill_between(asse,np.quantile(X,0.95,1),np.mean(X,1),alpha=0.3,color='green')
plt.fill_between(asse,np.mean(X,1),np.quantile(X,0.05,1),alpha=0.3,color='red')
plt.legend()
plt.show()

Altri formati di grafico sono, ovviamente, possibili. Per esempio vediamo di seguito come rappresentare 10’000 valori di una normale standard attraverso un istrogramma.

Y=np.random.normal(0,1,10000)
plt.figure()
plt.hist(Y,100) # Il secondo argomento è il numero di rettangoli dell'istogramma

## (array([  1.,   1.,   0.,   1.,   0.,   0.,   1.,   1.,   1.,   0.,   2.,
##          2.,   2.,   3.,   5.,   6.,   8.,   7.,   7.,  17.,  17.,  21.,
##         28.,  20.,  38.,  36.,  47.,  53.,  63.,  91.,  78., 105., 102.,
##        132., 133., 143., 170., 168., 199., 214., 236., 248., 277., 265.,
##        260., 295., 316., 311., 300., 306., 312., 336., 294., 319., 302.,
##        316., 301., 270., 273., 278., 223., 242., 191., 201., 192., 146.,
##        130., 129., 122., 102.,  85.,  87.,  77.,  55.,  39.,  41.,  41.,
##         31.,  27.,  17.,  20.,  18.,  12.,   4.,   5.,   6.,   5.,   3.,
##          2.,   1.,   3.,   1.,   1.,   0.,   1.,   0.,   0.,   1.,   0.,
##          1.]), array([-4.16704926, -4.08478728, -4.00252531, -3.92026333, -3.83800135,
##        -3.75573938, -3.6734774 , -3.59121542, -3.50895345, -3.42669147,
##        -3.34442949, -3.26216752, -3.17990554, -3.09764356, -3.01538158,
##        -2.93311961, -2.85085763, -2.76859565, -2.68633368, -2.6040717 ,
##        -2.52180972, -2.43954775, -2.35728577, -2.27502379, -2.19276182,
##        -2.11049984, -2.02823786, -1.94597588, -1.86371391, -1.78145193,
##        -1.69918995, -1.61692798, -1.534666  , -1.45240402, -1.37014205,
##        -1.28788007, -1.20561809, -1.12335612, -1.04109414, -0.95883216,
##        -0.87657019, -0.79430821, -0.71204623, -0.62978425, -0.54752228,
##        -0.4652603 , -0.38299832, -0.30073635, -0.21847437, -0.13621239,
##        -0.05395042,  0.02831156,  0.11057354,  0.19283551,  0.27509749,
##         0.35735947,  0.43962145,  0.52188342,  0.6041454 ,  0.68640738,
##         0.76866935,  0.85093133,  0.93319331,  1.01545528,  1.09771726,
##         1.17997924,  1.26224121,  1.34450319,  1.42676517,  1.50902714,
##         1.59128912,  1.6735511 ,  1.75581308,  1.83807505,  1.92033703,
##         2.00259901,  2.08486098,  2.16712296,  2.24938494,  2.33164691,
##         2.41390889,  2.49617087,  2.57843284,  2.66069482,  2.7429568 ,
##         2.82521878,  2.90748075,  2.98974273,  3.07200471,  3.15426668,
##         3.23652866,  3.31879064,  3.40105261,  3.48331459,  3.56557657,
##         3.64783854,  3.73010052,  3.8123625 ,  3.89462447,  3.97688645,
##         4.05914843]), <BarContainer object of 100 artists>)

plt.show()

Valgono anche per l’istogramma tutte le opzioni viste per i grafici precedenti. Nel seguente esempio vediamo come creare due istogrammi, uno sovrapposto all’altro, che rappresentano una normale e il suo quadrato. I grafici sono trasparenti in modo che si possano distinguere quando si sovrappongono.

plt.figure()
plt.hist(np.random.normal(0,1,10000),color='red',alpha=0.5,label='Gauss')

## (array([   3.,   47.,  314., 1173., 2444., 2945., 2044.,  816.,  189.,
##          25.]), array([-4.08348326, -3.31907183, -2.5546604 , -1.79024897, -1.02583753,
##        -0.2614261 ,  0.50298533,  1.26739676,  2.03180819,  2.79621962,
##         3.56063105]), <BarContainer object of 10 artists>)

plt.hist(np.random.normal(0,1,10000)**2,color='blue',alpha=0.5,label='Gauss^2')

## (array([8.387e+03, 1.141e+03, 3.150e+02, 9.900e+01, 3.700e+01, 1.400e+01,
##        4.000e+00, 2.000e+00, 0.000e+00, 1.000e+00]), array([1.03337384e-11, 1.94571980e+00, 3.89143960e+00, 5.83715940e+00,
##        7.78287920e+00, 9.72859900e+00, 1.16743188e+01, 1.36200386e+01,
##        1.55657584e+01, 1.75114782e+01, 1.94571980e+01]), <BarContainer object of 10 artists>)

plt.legend()
plt.show()

I grafici possono essere salvati, in un formato a piacere (png, jpg, pdf) mediante il comando:

plt.savefig(‘nome’)

Il grafico è salvato nella cartella di default su cui opera Colab. Vedremo tra poco come cambiare la cartella.

plt.figure()
plt.hist(np.random.normal(0,1,10000),color='red',alpha=0.5,label='Gauss')

## (array([   8.,   72.,  328., 1282., 2489., 2804., 1994.,  803.,  196.,
##          24.]), array([-3.90697028, -3.16627078, -2.42557128, -1.68487179, -0.94417229,
##        -0.20347279,  0.53722671,  1.2779262 ,  2.0186257 ,  2.7593252 ,
##         3.50002469]), <BarContainer object of 10 artists>)

plt.hist(np.random.normal(0,1,10000)**2,color='blue',alpha=0.5,label='Gauss^2')

## (array([7.727e+03, 1.421e+03, 4.800e+02, 2.050e+02, 9.200e+01, 3.300e+01,
##        2.800e+01, 9.000e+00, 3.000e+00, 2.000e+00]), array([3.41383351e-08, 1.47226570e+00, 2.94453137e+00, 4.41679704e+00,
##        5.88906271e+00, 7.36132838e+00, 8.83359405e+00, 1.03058597e+01,
##        1.17781254e+01, 1.32503911e+01, 1.47226567e+01]), <BarContainer object of 10 artists>)

plt.legend()
plt.savefig('istogrammi.pdf')

Salvare e richiamare i dati

Esistono diversi pacchetti per salvare e richiamare i dati. Il più “duttile” è pandas che utilizzeremo anche per richiamare dati da databases esterni. Il primo passo sarà quello di trasformare la variabile desiderata in un data frame del pacchetto pandas e poi salvarlo su un file.

import pandas as pd
df = pd.DataFrame(RW(100,3))
print(df)

##        0    1    2
## 0    0.0  0.0  0.0
## 1   -1.0  1.0  1.0
## 2    0.0  0.0  0.0
## 3   -1.0  1.0 -1.0
## 4    0.0  0.0  0.0
## ..   ...  ...  ...
## 95 -21.0  3.0 -7.0
## 96 -20.0  4.0 -8.0
## 97 -19.0  3.0 -7.0
## 98 -20.0  2.0 -6.0
## 99 -21.0  3.0 -5.0
## 
## [100 rows x 3 columns]

Per salvare il contenuto del dataframe usiamo il comando nella forma seguente:

variabile.to_csv('nome file.csv')

Al comando si può aggiungere l’opzione index=False che ipedisce all’operazione di salvataggio di usare la prima colonna come lista degli indici delle osservazioni (cioè i valori 0, 1, 2, e così via).

df.to_csv('RW.csv',index=False)

variabile = pd.read_csv('RW.csv')
print(variabile)

##        0    1    2
## 0    0.0  0.0  0.0
## 1   -1.0  1.0  1.0
## 2    0.0  0.0  0.0
## 3   -1.0  1.0 -1.0
## 4    0.0  0.0  0.0
## ..   ...  ...  ...
## 95 -21.0  3.0 -7.0
## 96 -20.0  4.0 -8.0
## 97 -19.0  3.0 -7.0
## 98 -20.0  2.0 -6.0
## 99 -21.0  3.0 -5.0
## 
## [100 rows x 3 columns]

plt.plot(variabile)

Sotto-grafici

Un caso classico è quello di più grafici che si vogliono rappresentare sullo tesso schema, dividendo il “plot” in “sub-plot”.

Il comando per gestire questo caso è subplot.

In particolare la sintassi è:

plt.subplot(r,c,a)

dove \(r\) indica il numero di righe, \(c\) il numero di colonne e \(a\) (compreso tra 1 e \(r\times c\)) è il numero che indica quale dei grafici è “attivo”. La numerazione dei grafici avviene dall’alto a sinistra verso il basso a destra. Nel caso di una matrice \(2\times 2\), per esempio, il numero \(a\) di ogni grafico è:

1 | 2 |
3 | 4 |

plt.figure()
plt.subplot(2,2,1)
plt.plot(X)
plt.plot(np.mean(X,1),color='black',linewidth=3,label='Media')
plt.subplot(2,2,2)
plt.plot(X**2)
plt.plot(np.mean(X**2,1),color='black',linewidth=3,label='Media')
plt.subplot(2,2,3)
plt.plot(X**3)
plt.plot(np.mean(X**3,1),color='black',linewidth=3,label='Media')
plt.show()

Possiamo inserire un grafico che occupi più di uno spazio mediante il comando:

plt.subplot(r,c,(a,b))

In questo caso si crea una griglia con \(r\) righe e \(c\) colonne e si attivano, contemporaneamente, i grafici che occupano le posizioni \(a\) e \(b\).

plt.figure()
plt.subplot(2,2,(1,2))
plt.plot(X)
plt.plot(np.mean(X,1),color='black',linewidth=3,label='Media')
plt.subplot(2,2,3)
plt.plot(X**2)
plt.plot(np.mean(X**2,1),color='black',linewidth=3,label='Media')
plt.subplot(2,2,4)
plt.plot(X**3)
plt.plot(np.mean(X**3,1),color='black',linewidth=3,label='Media')
plt.show()

Si possono aggiungere opzioni ai singoli grafici come già visto in precedenza (titoli, griglie, legende, etichette).

Esiste, tuttavia, la possibilità di dare un titolo “complessivo” al grafico mediante li comando suptitle.

plt.figure()

plt.subplot(2,2,(1,2))
plt.plot(X,color='lightgray')
plt.plot(np.mean(X,1),color='black',linewidth=3,label='Media')
plt.plot(np.quantile(X,0.95,1),color='green',linewidth=3,label='Quantile 0.95')
plt.plot(np.quantile(X,0.05,1),color='red',linewidth=3,label='Quantile 0.05')
plt.title('RW')

plt.subplot(2,2,3)
plt.plot(X**2,color='lightgray')
plt.plot(np.mean(X**2,1),color='black',linewidth=3,label='Media')
plt.plot(np.quantile(X**2,0.95,1),color='green',linewidth=3,label='Quantile 0.95')
plt.plot(np.quantile(X**2,0.05,1),color='red',linewidth=3,label='Quantile 0.05')
plt.title('RW quadrato')

plt.subplot(2,2,4)
plt.plot(X**3,color='lightgray')
plt.plot(np.mean(X**3,1),color='black',linewidth=3,label='Media')
plt.plot(np.quantile(X**3,0.95,1),color='green',linewidth=3,label='Quantile 0.95')
plt.plot(np.quantile(X**3,0.05,1),color='red',linewidth=3,label='Quantile 0.05')
plt.title('RW cubo')

plt.tight_layout() # Comando per spaziare opportunament i grafici tra di loro
plt.legend(bbox_to_anchor=(1,1.2),loc='center left')

plt.savefig('subplot.pdf',bbox_inches='tight') # L'opzione "tight" serve per salvare la legenda dentro al PDF

Vediamo una solzione un po’ diversa per disporre i grafici

plt.figure()

plt.subplot(2,2,1)
plt.plot(X,color='lightgray')
plt.plot(np.mean(X,1),color='black',linewidth=3,label='Media')
plt.plot(np.quantile(X,0.95,1),color='green',linewidth=3,label='Quantile 0.95')
plt.plot(np.quantile(X,0.05,1),color='red',linewidth=3,label='Quantile 0.05')
plt.title('RW')

plt.subplot(2,2,2)
plt.plot(X**2,color='lightgray')
plt.plot(np.mean(X**2,1),color='black',linewidth=3,label='Media')
plt.plot(np.quantile(X**2,0.95,1),color='green',linewidth=3,label='Quantile 0.95')
plt.plot(np.quantile(X**2,0.05,1),color='red',linewidth=3,label='Quantile 0.05')
plt.title('RW quadrato')

plt.subplot(2,2,3)
plt.plot(X**3,color='lightgray')
plt.plot(np.mean(X**3,1),color='black',linewidth=3,label='Media')
plt.plot(np.quantile(X**3,0.95,1),color='green',linewidth=3,label='Quantile 0.95')
plt.plot(np.quantile(X**3,0.05,1),color='red',linewidth=3,label='Quantile 0.05')
plt.title('RW cubo')

plt.tight_layout() # Comando per spaziare opportunament i grafici tra di loro
plt.legend(bbox_to_anchor=(1.45,0.5),loc='center left')

plt.show()

Unità di misura per assi e griglie

plt.figure()
plt.plot(X)
plt.show()

plt.figure()
plt.plot(X)
plt.xticks(np.arange(0, 250, 10),rotation=90)

## ([<matplotlib.axis.XTick object at 0x7f02a1c98be0>, <matplotlib.axis.XTick object at 0x7f02a1c98bb0>, <matplotlib.axis.XTick object at 0x7f02a1c98a90>, <matplotlib.axis.XTick object at 0x7f02a1bc2680>, <matplotlib.axis.XTick object at 0x7f02a1bc29b0>, <matplotlib.axis.XTick object at 0x7f02a1bc3100>, <matplotlib.axis.XTick object at 0x7f02a1bc3160>, <matplotlib.axis.XTick object at 0x7f02a1bc3730>, <matplotlib.axis.XTick object at 0x7f02a1bc3e80>, <matplotlib.axis.XTick object at 0x7f02a1bfc610>, <matplotlib.axis.XTick object at 0x7f02a1bfcd60>, <matplotlib.axis.XTick object at 0x7f02a1bfd4b0>, <matplotlib.axis.XTick object at 0x7f02a1bc3ee0>, <matplotlib.axis.XTick object at 0x7f02a1bfdae0>, <matplotlib.axis.XTick object at 0x7f02a1bfc5e0>, <matplotlib.axis.XTick object at 0x7f02a1bfe1a0>, <matplotlib.axis.XTick object at 0x7f02a1bfe8f0>, <matplotlib.axis.XTick object at 0x7f02a1bff040>, <matplotlib.axis.XTick object at 0x7f02a1bff790>, <matplotlib.axis.XTick object at 0x7f02a1bffcd0>, <matplotlib.axis.XTick object at 0x7f02a1bfe650>, <matplotlib.axis.XTick object at 0x7f02a1bfdba0>, <matplotlib.axis.XTick object at 0x7f02a1c14430>, <matplotlib.axis.XTick object at 0x7f02a1c14c10>, <matplotlib.axis.XTick object at 0x7f02a1c15300>], [Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, '')])

plt.yticks(np.arange(-40, 50, 10))

## ([<matplotlib.axis.YTick object at 0x7f02a1c99ae0>, <matplotlib.axis.YTick object at 0x7f02a1c99360>, <matplotlib.axis.YTick object at 0x7f02a1c98ac0>, <matplotlib.axis.YTick object at 0x7f02a1bfe860>, <matplotlib.axis.YTick object at 0x7f02a1c150c0>, <matplotlib.axis.YTick object at 0x7f02a1c16170>, <matplotlib.axis.YTick object at 0x7f02a1c16650>, <matplotlib.axis.YTick object at 0x7f02a1c16da0>, <matplotlib.axis.YTick object at 0x7f02a1c174f0>], [Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, '')])

plt.grid(linestyle='dotted')
plt.show()

Per modificare la scala del grafico (per le ascisse o per le ordinate) possiamo utilizzare il comando xscale oppure yscale.

plt.figure()
plt.plot(X**2)
plt.yscale('log')
plt.show()

Soluzione numerica di equazioni

Alcune equazioni, dette tracendenti, non sono risolvibili algebricamente. Un caso è la seguente equazione \[x=e^{-x}.\]

Nonostante non esista una soluzione algebrica, sappiamo che questa equazione ha una soluzione ed è unica.

Il primo passo è quello di definire una funzione in modo che il suo valore sia nullo e poi utilizzare un comando per trovare la soluzione.

# Il pacchetto per trovare le soluzioni è il seguente
import scipy.optimize as sci

def funzione(x):
  return x-np.exp(-x)

x=np.arange(-2,2,0.01)
plt.figure()
plt.plot(x,funzione(x))
plt.axhline(y = 0, color = 'red', linestyle = 'dotted')
plt.grid()
plt.show()

# Troviamo la soluzione
soluzione=sci.fsolve(funzione,0.5)
soluzione

## array([0.56714329])

x=np.arange(-2,2,0.01)
plt.figure()
plt.plot(x,funzione(x))
plt.axhline(y = 0, color = 'red', linestyle = 'dotted')
plt.axvline(x = soluzione, color = 'red', linestyle = 'dotted')
plt.grid()
plt.show()

Minimizzare una funzione

Supponiamo di avere la seguente funzione \[f(x)=e^{-x}+x,\] e di volerla minimizzare.

La derivata prima è \[\frac{\partial f(x)}{\partial x}=-e^{-x}+1,\] e la derivata seconda è \[\frac{\partial^2 f(x)}{\partial x^2}=e^{-x}>0.\] Poiché la derivata seconda è sempre positiva, allora il valore di \(x\) per cui la derivata prima è nulla è un minimo: \[x^*=0\]

Proviamo a verificare questo risultato numericamente.

# Definiamo la funzione
def funzione(x):
  return np.exp(-x)+x

# Rappresentiamo graficamente la funzione
asse=np.arange(-2,2,0.01)
plt.figure()
plt.plot(asse,funzione(asse))
plt.grid()
plt.show()

ottimo=sci.minimize(funzione,2.5)
ottimo

##       fun: 1.0000000000000004
##  hess_inv: array([[1.00059591]])
##       jac: array([-2.98023224e-08])
##   message: 'Optimization terminated successfully.'
##      nfev: 18
##       nit: 8
##      njev: 9
##    status: 0
##   success: True
##         x: array([-2.56952413e-08])

Nel risultato osserviamo le seguenti variabili: - fun: è il valore della funzione nel punto di minimo - hass_inv: è la matrice hessiana inversa (è la matrice delle derivate seconde della funzione, calcolate nel punto di minimo) - jac: è il vettore jacobiano che contiene le derivate prime calcolate nel punto di minimo (dovrebbe essere un vettore di zeri) - nit: è il numero di iterazioni effettuate - x: è il numero che minimizza la funzione obiettivo

ottimo.x

## array([-2.56952413e-08])

Per massimizzare una funzione è sufficiente minimizzarne l’opposto. Cioè il \[\max f(x)\] equivale al \[\min -f(x).\] Questo significa che la \(x\) che risolve il primo problema ha lo stesso valore della \(x\) che risolve il secondo problema.

Lavorare su database esterni

Utilizziamo nuovamente il pacchettp pandas per richiamare dall’esterno i dati da due database: FRED e Yahoo.

Per esempio su FRED prendiamo la serie dei tassi di interesse a 3 mesi dei Treasury Bill (il codice è DTB3).

import pandas_datareader as pdr
tasso=pdr.DataReader('DTB3','fred')
print(tasso)

##             DTB3
## DATE            
## 2017-09-25  1.04
## 2017-09-26  1.05
## 2017-09-27  1.05
## 2017-09-28  1.04
## 2017-09-29  1.04
## ...          ...
## 2022-09-14  3.14
## 2022-09-15  3.12
## 2022-09-16  3.10
## 2022-09-19  3.29
## 2022-09-20  3.26
## 
## [1302 rows x 1 columns]

plt.plot(tasso)

Con questi comandi si scaricano solo gli “ultimi” dati. Per scaricare più dati occorre inidcare da quale data partire. Per lavorare sulle date occorre caricare un nuovo pacchetto!

import datetime
tasso=pdr.DataReader('DTB3','fred',start=datetime.datetime(1900,1,1))
plt.plot(tasso)

Calcoliam il tasso medio

np.mean(tasso)

## DTB3    4.172805
## dtype: float64
## 
## /usr/lib/python3/dist-packages/numpy/core/fromnumeric.py:3438: FutureWarning: In a future version, DataFrame.mean(axis=None) will return a scalar mean over the entire DataFrame. To retain the old behavior, use 'frame.mean(axis=0)' or just 'frame.mean()'
##   return mean(axis=axis, dtype=dtype, out=out, **kwargs)

Posso fare la media SOLO della variabile DTB3

np.mean(tasso.DTB3)

## 4.172805147915208

Togliamo i valori NaN

tasso=tasso.dropna()/100

Inseriamo na linea orizzontale in corripondenza della media

plt.figure()
plt.plot(tasso)
plt.axhline(np.mean(tasso.DTB3),color='red',linestyle='dotted')
plt.axvline(datetime.datetime(2008,9,15),color='black',linestyle='dashed')
plt.text(datetime.datetime(2008,9,15),max(tasso.DTB3),'Lehman',horizontalalignment='center',fontstyle='italic',fontsize='small')
plt.show()

Si possono identificare gli elementi della serie storica che godono di una certa proprietà. Per esempio vediamo se (e dove) esistono dei tassi di interesse nulli, oppure negativi.

np.where(tasso==0)

## (array([13724, 13730, 13734, 14422, 14431, 14479, 14480, 14925, 15430,
##        15432, 15433, 15437, 15443]), array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]))

np.where(tasso<0)

## (array([15420, 15422, 15425, 15428, 15429, 15434, 16547, 16548]), array([0, 0, 0, 0, 0, 0, 0, 0]))

np.where(tasso<0)[0]

## array([15420, 15422, 15425, 15428, 15429, 15434, 16547, 16548])

tasso.DTB3[np.where(tasso<0)[0]]

## DATE
## 2015-09-18   -0.0001
## 2015-09-22   -0.0001
## 2015-09-25   -0.0001
## 2015-09-30   -0.0001
## 2015-10-01   -0.0002
## 2015-10-08   -0.0001
## 2020-03-25   -0.0004
## 2020-03-26   -0.0005
## Name: DTB3, dtype: float64

Per utilizzare i dati al di fuori della forma della serie storica, si può trasformare la serie in un oggetto “numpy”:

r=pd.Series.to_numpy(tasso.DTB3)
r

## array([0.0133, 0.0128, 0.0128, ..., 0.031 , 0.0329, 0.0326])

Si può modificare la frequenza dei dati.

# frequenza mensile "m" (metodo della media - si può usare anche la somma "sum" come per il PIL)
tasso.resample('m').mean()

##                 DTB3
## DATE                
## 1954-01-31  0.011800
## 1954-02-28  0.009689
## 1954-03-31  0.010322
## 1954-04-30  0.009664
## 1954-05-31  0.007630
## ...              ...
## 2022-05-31  0.009833
## 2022-06-30  0.014948
## 2022-07-31  0.022325
## 2022-08-31  0.026309
## 2022-09-30  0.030623
## 
## [825 rows x 1 columns]

Adesso scarichiamo due seire di dati con lo stesso comando.

tassi=pdr.DataReader(['IRLTLT01USM156N','TB3MS'],'fred',start=datetime.datetime(2000,1,1))
plt.figure()
plt.plot(tassi)
plt.legend(['Tasso 10a','Tasso 3m'])
plt.show()

Sapendo che la differenza tra i due tassi di interesse è una variabile che prevee le crisi, possiamo rapppresentare tale differenza su un grafico.

plt.figure()
plt.plot(tassi.IRLTLT01USM156N-tassi.TB3MS)
plt.axhline(y=0,color='red',linestyle='dotted')
plt.grid()
plt.show()

Regressione lineare

# Importiamo il pacchetto statistico
import statsmodels.api as stat

In questo caso prendo il PIL e i consumi italiani (dal database FRED) e faccio una regressione OLS del consumo sul PIL.

dati=pdr.DataReader(['CPMNACSCAB1GQIT','ITAPFCEQDSMEI'],'fred',start=datetime.datetime(2000,1,1))
dati

##             CPMNACSCAB1GQIT  ITAPFCEQDSMEI
## DATE                                      
## 2000-01-01         304165.2   1.835236e+11
## 2000-04-01         309195.6   1.868507e+11
## 2000-07-01         311720.7   1.891021e+11
## 2000-10-01         318195.8   1.920102e+11
## 2001-01-01         322663.4   1.926422e+11
## ...                     ...            ...
## 2021-04-01         439622.4   2.540863e+11
## 2021-07-01         452305.7   2.633786e+11
## 2021-10-01         453668.0   2.670711e+11
## 2022-01-01         460184.6   2.706622e+11
## 2022-04-01         471524.0   2.820270e+11
## 
## [90 rows x 2 columns]

I dati non hanno la stessa unità di misura. I dati della seconda colonna sono in Euro, mentre gli altri sono in milioni di Euro. Mi conviene, quindi, dividere per 1’000’000 i dati della seconda colonna:

dati.ITAPFCEQDSMEI=dati.ITAPFCEQDSMEI/1000000
dati

##             CPMNACSCAB1GQIT  ITAPFCEQDSMEI
## DATE                                      
## 2000-01-01         304165.2       183523.6
## 2000-04-01         309195.6       186850.7
## 2000-07-01         311720.7       189102.1
## 2000-10-01         318195.8       192010.2
## 2001-01-01         322663.4       192642.2
## ...                     ...            ...
## 2021-04-01         439622.4       254086.3
## 2021-07-01         452305.7       263378.6
## 2021-10-01         453668.0       267071.1
## 2022-01-01         460184.6       270662.2
## 2022-04-01         471524.0       282027.0
## 
## [90 rows x 2 columns]

# Definisco la variabile "consumo" C e la variabile "reddito" Y
C=dati.ITAPFCEQDSMEI
Y=dati.CPMNACSCAB1GQIT

# Chiedo a Python di aggiungere la costante di regressione
Y = stat.add_constant(Y)

# Disegno la nuvola di punti del consumo-reddito
plt.figure()
plt.plot(dati.CPMNACSCAB1GQIT,dati.ITAPFCEQDSMEI,marker='o',linestyle='')
plt.show()

# Ottengo i risultati della regressione lineare
model = stat.OLS(C,Y)
results = model.fit()
print(results.summary())

##                             OLS Regression Results                            
## ==============================================================================
## Dep. Variable:          ITAPFCEQDSMEI   R-squared:                       0.970
## Model:                            OLS   Adj. R-squared:                  0.969
## Method:                 Least Squares   F-statistic:                     2827.
## Date:                gio, 22 set 2022   Prob (F-statistic):           1.12e-68
## Time:                        18:40:57   Log-Likelihood:                -876.16
## No. Observations:                  90   AIC:                             1756.
## Df Residuals:                      88   BIC:                             1761.
## Df Model:                           1                                         
## Covariance Type:            nonrobust                                         
## ===================================================================================
##                       coef    std err          t      P>|t|      [0.025      0.975]
## -----------------------------------------------------------------------------------
## const             162.0178   4490.349      0.036      0.971   -8761.608    9085.644
## CPMNACSCAB1GQIT     0.5980      0.011     53.167      0.000       0.576       0.620
## ==============================================================================
## Omnibus:                       21.966   Durbin-Watson:                   0.171
## Prob(Omnibus):                  0.000   Jarque-Bera (JB):               32.244
## Skew:                          -1.074   Prob(JB):                     9.96e-08
## Kurtosis:                       4.996   Cond. No.                     4.11e+06
## ==============================================================================
## 
## Notes:
## [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
## [2] The condition number is large, 4.11e+06. This might indicate that there are
## strong multicollinearity or other numerical problems.