Introduzione
Il pacchetto scipy utilizza la funzionalità di numpy per fornire un pacchetto di calcolo scientifico general purpose.
Scipy è in realtà una collezione enorme di sottopacchetti:
Le estensioni offerte
- scipy.constants: costanti matematiche e fisiche
- scipy.special: funzioni in uso in fisica matematica (ellittiche, Bessel, ipergeometriche)
- scipy.integrate: metodi di integrazione numerica(trapezoidale, Simpson), integrazione di equazioni differenziali
- scipy.optimize: metodi di ottimizzazione (minimi quadrati, gradiente, simulated annealing)
- scipy.linalg: estensione di numpy.linalg; soluzione di sistemi lineari, calcolo matriciale, decomposizione, fattorizzazione
- scipy.sparse: gestione di matrici sparse
- scipy.interpolate: metodi per l’interpolazione lineare e non di dati
- scipy.fftpack: Fast Fourier Transform
- scipy.signal: metodi di signal processing (filtraggio, correlazione, convoluzione, smoothing)
- scipy.stats: distribuzioni di probabilità continue e discrete, calcolo dei momenti, calcolo cumulative, statistica descrittiva, test
Vediamone alcune:
esempio: generazione di numeri casuali
import scipy.stats q = scipy.stats.beta(5, 5) # genera una beta(5,5) obs = q.rvs(2000) # produce 2000 osservazioni print(obs) ⇒ [0.5955691 0.55416895 0.34151925 ... 0.37947352 0.60492705 0.3533784 ] # Stampiamo statistiche sull'insieme print obs.min() 0.0749989919902 print obs.max() 0.919066721448 print obs.std() 0.152290115168 print obs.mean() 0.506227887253
Usiamo il pacchetto scipy.stats per effettuare una regressione lineare
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
np.random.seed(12345678)
x = np.random.random(10)
y = np.random.random(10)
slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)
print("Pendenza:", slope)
print("coefficiente angolare: " ,gradient) #: coefficiente angolare
print("intersezione con asse Y: ",intercept)#: intersezione con asse Y
print("radice quadrata del coefficiente di correlazione: ",r_value) #: radice quadrata del coefficiente di correlazione
print("Test statistico sull'ipotesi nulla “il coefficiente angolare della retta di regressione è zero: ",p_value) #: test statistico sull'ipotesi nulla “il coefficiente angolare della retta di regressione è zero”
print(" errore standard della stima: ",std_err) #: errore standard della stima
print("r-quadro:", r_value**2)
plt.plot(x, y, 'o', label='original data')
plt.plot(x, intercept + slope*x, 'r', label='fitted line')
plt.legend()
plt.show()
Output:
File I/O: modulo scipy.io
import numpy as np
from scipy import io
a = np.ones((3,3))
io.savemat("file.mat", {'a': a})
data = io.loadmat("file.mat", struct_as_record=True)
print(data["a"])
# leggere un'immagine
print("leggere un'immagine: ")
from matplotlib.pyplot import imread
im = imread("/content/provalettura.png")
print(im)
Calcolo del determinante di una matrice
from scipy import linalg
arr = np.array([[1, 2],[3, 4]])
print (" il Det. è: " , linalg.det(arr) )
arr = np.array([[3, 2], [6, 4]])
print(" il Det. è: ",linalg.det(arr))
Output:
il Det. è: -2.0
il Det. è: 0.0
Ovviamente se la matrice non è quadrata, restituisce errore:
ValueError: expected square matrix
Calcolo dell’inversa di una matrice
Matrice inversa
Ricordiamo che una matrice quadrata A si dice invertibile se esiste una matrice quadrata B, dello stesso ordine, tale che AB = BA = I.
from scipy import linalg arr = np.array([[1, 2],[3, 4]]) iarr = linalg.inv(arr) iarr Output: array([[-2. , 1. ], [ 1.5, -0.5]])
Soluzione di un sistema lineare
Data una matrice A(n,m) ed un vettore b(n), calcolare il vettore x(n) tale che Ax= b.
from scipy import linalg A = np.array([[1,2],[3,4]]) b = np.array([[5],[6]]) x = linalg.solve(A,b) x Output: array([[-4. ], [ 4.5]])
Ottimizzazione e fit: scipy.optimize
In generale, l’ottimizzazione consiste nella ricerca di una soluzione numerica a un problema di minimizzazione o nel calcolo degli zeri di una funzione
Ricerca del minimo di una funzione: definizione e grafico
def f(x): return x**2 + 10*np.sin(x) import matplotlib.pyplot as plt x = np.arange(-10, 10, 0.1) plt.plot(x, f(x)) plt.show() 
La funzione ha un minimo globale intorno a -1 ed un minimo locale intorno a 4!
Ricerca del minimo assoluto
Un metodo generale ed efficiente per la ricerca del minimo della funzione scelta consiste nel condurre una discesa lungo il gradiente, partendo da un punto iniziale. L’algoritmo BFGS è una valida soluzione del problema!
import numpy as np
from scipy import optimize
optimize.fmin_bfgs(f, 0)
Optimization terminated successfully.
Current function value: -7.945823
Iterations: 5
Function evaluations: 18
Gradient evaluations: 6
array([-1.30644012])
Un possibile problema di questo approccio è che, se la funzione presenta minimi locali, l’algoritmo può giungere ad uno di questi ultimi, invece che al minimo assoluto, per certe scelte dello starting point della discesa.
Ricerca del minimo con l’algoritmo BFGS (starting point x=3)
import numpy as np from scipy import optimize optimize.fmin_bfgs(f, 3, disp=0)
grid = (-10, 10, 0.1)
optimize.brute (f, (grid,),disp=True)
Optimization terminated successfully.
Current function value: -7.945823
Iterations: 11
Function evaluations: 22
array([-1.30641113]
Ricerca degli zeri di una funzione
La funzione scipy.optimize.fsolve consente di trovare lo zero di una funzione, ovvero una soluzione per l’equazione f(x) =0
Per semplicità utilizziamo la funzione vista per la ricerca del minimo assoluto.
Questo algoritmo richiede solamente un punto di partenza vicino a dove ci si aspetta che ci sia una radice (Non è però detto che il metodo converga).
from scipy import optimize x_1 = optimize.fsolve(f, 1) #starting point 1 x_2 = optimize.fsolve(f,-3) #starting point -3 print(x_1 ) print(x_2)
from scipy.optimize import fsolve
def f(x):
return x ** 3 - 2 * x ** 2
xstart = 3
x = fsolve(f, xstart) # one root is at x=2.0
# search starts at x=3
# fsolve returns a numpy array
print(f"The root x is approximately x={x[0]:21.19g}")
print(f"The exact error is {2 - x[0]:g}.")
The root x is approximately x= 2.000000000000006661 The exact error is -6.66134e-15.
[elementor-template id=”12586″]
(1228)