Introduction to SciPy
The SciPy library (pronounced “Sigh Pie”) is built on NumPy and provides many user-friendly and efficient numerical routines, such as routines for numerical integration, interpolation, optimization, linear algebra, and statistics.
While NumPy provides the data structure (arrays) and basic operations, SciPy provides the algorithms.
Optimization and Root Finding
The scipy.optimize subpackage provides several commonly used optimization algorithms.
Minimizing a Scalar Function
Let’s find the minimum of the Rosenbrock function, a common test problem for optimization algorithms.
Interactive Lab
from scipy.optimize import minimize
import numpy as np
def rosen(x):
"""The Rosenbrock function"""
return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1.0-x[:-1])**2.0)
x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
res = minimize(rosen, x0, method='nelder-mead', options={'xatol': 1e-8, 'disp': True})
print(f"Minimum found at: {res.x}")
Expected output
Minimum found at: [1. 1. 1. 1. 1.]
1
2from scipy.optimize import minimize
3import numpy as np
4
5def rosen(x):
6 """The Rosenbrock function"""
7 return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1.0-x[:-1])**2.0)
8
9x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
10res = minimize(rosen, x0, method='nelder-mead', options={'xatol': 1e-8, 'disp': True})
11
12print(f"Minimum found at: {res.x}")
13
Numerical Integration
The scipy.integrate subpackage provides several integration techniques, including an ordinary differential equation (ODE) integrator.
General Integration (quad)
Interactive Lab
from scipy.integrate import quad
import numpy as np
# Integrate exp(-x^2) from 0 to infinity (The Gaussian integral)
func = lambda x: np.exp(-x**2)
result, error = quad(func, 0, np.inf)
print(f"Result: {result}")
print(f"Error Estimate: {error}")
Expected output
Result: 0.886226925452758
Error Estimate: 9.839070275988184e-11
1
2from scipy.integrate import quad
3import numpy as np
4
5
6func = lambda x: np.exp(-x**2)
7result, error = quad(func, 0, np.inf)
8
9print(f"Result: {result}")
10print(f"Error Estimate: {error}")
11
Interpolation
scipy.interpolate is useful for fitting a function to a set of data points and then evaluating that function at new points.
Interactive Lab
from scipy.interpolate import interp1d
import numpy as np
x = np.linspace(0, 10, num=11, endpoint=True)
y = np.cos(-x**2/9.0)
# Create linear and cubic interpolation functions
f_linear = interp1d(x, y)
f_cubic = interp1d(x, y, kind='cubic')
xnew = np.linspace(0, 10, num=41, endpoint=True)
# Evaluate at new points
print(f"Linear interp at 5.5: {f_linear(5.5)}")
print(f"Cubic interp at 5.5: {f_cubic(5.5)}")
Expected output
Linear interp at 5.5: -0.6695277561841323
Cubic interp at 5.5: -0.6212457813289297
1
2from scipy.interpolate import interp1d
3import numpy as np
4
5x = np.linspace(0, 10, num=11, endpoint=True)
6y = np.cos(-x**2/9.0)
7
8
9f_linear = interp1d(x, y)
10f_cubic = interp1d(x, y, kind='cubic')
11
12xnew = np.linspace(0, 10, num=41, endpoint=True)
13
14print(f"Linear interp at 5.5: {f_linear(5.5)}")
15print(f"Cubic interp at 5.5: {f_cubic(5.5)}")
16
Signal Processing
scipy.signal contains tools for filtering, spectral analysis, and LTI (linear time-invariant) system analysis.
Interactive Lab
from scipy import signal
import numpy as np
# Create a noisy signal
t = np.linspace(0, 1, 100)
sig = np.sin(2*np.pi*10*t) + np.random.randn(100)*0.1
# Apply a Butterworth filter
b, a = signal.butter(3, 0.05)
filtered_sig = signal.filtfilt(b, a, sig)
print("Signal filtered using Butterworth filter.")
Expected output
Signal filtered using Butterworth filter.
1
2from scipy import signal
3import numpy as np
4
5
6t = np.linspace(0, 1, 100)
7sig = np.sin(2*np.pi*10*t) + np.random.randn(100)*0.1
8
9
10b, a = signal.butter(3, 0.05)
11filtered_sig = signal.filtfilt(b, a, sig)
12
13print("Signal filtered using Butterworth filter.")
14
Statistics in SciPy
scipy.stats contains a huge number of probability distributions and statistical functions.
Interactive Lab
from scipy import stats
import numpy as np
# Normal distribution
data = np.random.normal(0, 1, 1000)
# Calculate parameters
mean, std = stats.norm.fit(data)
print(f"Fitted Mean: {mean:.2f}, Std: {std:.2f}")
# Perform a T-test
t_stat, p_val = stats.ttest_1samp(data, 0.0)
print(f"T-statistic: {t_stat:.2f}, P-value: {p_val:.2f}")
Expected output
Fitted Mean: 0.01, Std: 0.99
T-statistic: 0.32, P-value: 0.75
1
2from scipy import stats
3import numpy as np
4
5
6data = np.random.normal(0, 1, 1000)
7
8
9mean, std = stats.norm.fit(data)
10print(f"Fitted Mean: {mean:.2f}, Std: {std:.2f}")
11
12
13t_stat, p_val = stats.ttest_1samp(data, 0.0)
14print(f"T-statistic: {t_stat:.2f}, P-value: {p_val:.2f}")
15
In the final module, we will explore how to visualize our findings using Matplotlib and Seaborn.