#Mathematical Operations and Functions
In this chapter, we'll dive deep into NumPy's mathematical operations and functions, including basic arithmetic, trigonometric functions, statistical functions, linear algebra, and more.
#Basic Arithmetic Operations
#Element-wise Operations
import numpy as np
# Create sample arrays
array1 = np.array([1, 2, 3, 4, 5])
array2 = np.array([10, 20, 30, 40, 50])
print("Array 1:", array1)
print("Array 2:", array2)
print()
# Basic arithmetic operations
print("Addition array1 + array2:", array1 + array2)
print("Subtraction array2 - array1:", array2 - array1)
print("Multiplication array1 * array2:", array1 * array2)
print("Division array2 / array1:", array2 / array1)
print("Integer division array2 // array1:", array2 // array1)
print("Modulo array2 % array1:", array2 % array1)
print("Power array1 ** 2:", array1 ** 2)
print()
# Operations with scalar
scalar = 10
print(f"Array plus scalar array1 + {scalar}:", array1 + scalar)
print(f"Array times scalar array1 * {scalar}:", array1 * scalar)
print(f"Array divided by scalar array1 / {scalar}:", array1 / scalar)
print(f"Scalar minus array {scalar} - array1:", scalar - array1)#Mathematical Functions
import numpy as np
# Create test array
array = np.array([1, 4, 9, 16, 25])
print("Original array:", array)
print()
# Basic math functions
print("Square root np.sqrt():", np.sqrt(array))
print("Cube root np.cbrt():", np.cbrt(array))
print("Square np.square():", np.square(array))
print("Absolute value np.abs():", np.abs(array - 10))
print()
# Exponential and logarithmic functions
small_array = np.array([1, 2, 3, 4, 5])
print("Small array:", small_array)
print("Natural exponential np.exp():", np.exp(small_array))
print("Base-2 exponential np.exp2():", np.exp2(small_array))
print("Natural logarithm np.log():", np.log(small_array))
print("Base-10 logarithm np.log10():", np.log10(small_array))
print("Base-2 logarithm np.log2():", np.log2(small_array))
print()
# Rounding functions
float_array = np.array([1.2, 2.7, -1.5, -2.8, 3.0])
print("Float array:", float_array)
print("Floor np.floor():", np.floor(float_array))
print("Ceiling np.ceil():", np.ceil(float_array))
print("Round np.round():", np.round(float_array))
print("Truncate np.trunc():", np.trunc(float_array))#Trigonometric Functions
#Basic Trigonometric Functions
import numpy as np
# Create angle array (radians)
angles_rad = np.array([0, np.pi/6, np.pi/4, np.pi/3, np.pi/2, np.pi])
print("Angles (radians):", angles_rad)
print("Angles (degrees):", np.degrees(angles_rad))
print()
# Basic trig functions
print("Sine np.sin():", np.sin(angles_rad))
print("Cosine np.cos():", np.cos(angles_rad))
print("Tangent np.tan():", np.tan(angles_rad))
print()
# Inverse trigonometric functions
values = np.array([0, 0.5, 0.707, 0.866, 1.0])
print("Values:", values)
print("Arcsine np.arcsin():", np.arcsin(values))
print("Arccosine np.arccos():", np.arccos(values))
print("Arctangent np.arctan():", np.arctan(values))
print()
# Hyperbolic functions
print("Hyperbolic sine np.sinh():", np.sinh([0, 1, 2]))
print("Hyperbolic cosine np.cosh():", np.cosh([0, 1, 2]))
print("Hyperbolic tangent np.tanh():", np.tanh([0, 1, 2]))#Angle Conversion and Special Functions
import numpy as np
# Angle and radian conversion
degrees = np.array([0, 30, 45, 60, 90, 180])
radians = np.radians(degrees)
print("Degrees:", degrees)
print("Converted to radians:", radians)
print("Converted back to degrees:", np.degrees(radians))
print()
# Special trigonometric functions
x = np.array([1, 2, 3])
y = np.array([1, 1, 4])
print("x:", x)
print("y:", y)
print("atan2(y, x):", np.arctan2(y, x)) # Quadrant-aware arctangent
print("hypot(x, y):", np.hypot(x, y)) # Calculate sqrt(x²+y²)
print()
# Angle normalization
angles = np.array([0, 90, 180, 270, 360, 450, -90])
print("Original angles:", angles)
print("Normalized to [0, 360):", angles % 360)
print("Normalized to [-180, 180):", ((angles + 180) % 360) - 180)#Statistical Functions
#Basic Statistics
import numpy as np
# Create test data
np.random.seed(42)
data = np.random.normal(50, 10, 20) # Normal distribution, mean 50, std 10
print("Test data:")
print(data)
print()
# Basic statistics
print(f"Mean np.mean(): {np.mean(data):.2f}")
print(f"Median np.median(): {np.median(data):.2f}")
print(f"Standard deviation np.std(): {np.std(data):.2f}")
print(f"Variance np.var(): {np.var(data):.2f}")
print(f"Minimum np.min(): {np.min(data):.2f}")
print(f"Maximum np.max(): {np.max(data):.2f}")
print(f"Range np.ptp(): {np.ptp(data):.2f}") # Peak to peak
print(f"Sum np.sum(): {np.sum(data):.2f}")
print(f"Product np.prod(): {np.prod(data):.2e}")
print()
# Percentiles
percentiles = [25, 50, 75, 90, 95]
print("Percentiles:")
for p in percentiles:
value = np.percentile(data, p)
print(f" {p}th percentile: {value:.2f}")#Multidimensional Array Statistics
import numpy as np
# Create 2D array
array_2d = np.random.randint(1, 10, (4, 5))
print("2D array:")
print(array_2d)
print()
# Global statistics
print(f"Global mean: {np.mean(array_2d):.2f}")
print(f"Global std: {np.std(array_2d):.2f}")
print()
# Statistics along axes
print("Row-wise statistics (axis=1):")
print(f" Mean per row: {np.mean(array_2d, axis=1)}")
print(f" Max per row: {np.max(array_2d, axis=1)}")
print(f" Min per row: {np.min(array_2d, axis=1)}")
print()
print("Column-wise statistics (axis=0):")
print(f" Mean per column: {np.mean(array_2d, axis=0)}")
print(f" Max per column: {np.max(array_2d, axis=0)}")
print(f" Min per column: {np.min(array_2d, axis=0)}")
print()
# Cumulative statistics
print("Cumulative sum np.cumsum():")
print(np.cumsum(array_2d, axis=1)) # Row-wise cumulative
print()
print("Cumulative product np.cumprod():")
print(np.cumprod(array_2d, axis=0)) # Column-wise cumulative#Correlation and Covariance
import numpy as np
# Create correlated data
np.random.seed(42)
x = np.random.normal(0, 1, 100)
y = 2 * x + np.random.normal(0, 0.5, 100) # y correlated with x
z = np.random.normal(0, 1, 100) # z uncorrelated with x, y
data = np.column_stack([x, y, z])
print(f"Data shape: {data.shape}")
print("First 5 rows:")
print(data[:5])
print()
# Correlation coefficient matrix
corr_matrix = np.corrcoef(data.T)
print("Correlation coefficient matrix:")
print(corr_matrix)
print()
# Covariance matrix
cov_matrix = np.cov(data.T)
print("Covariance matrix:")
print(cov_matrix)
print()
# Correlation between two variables
corr_xy = np.corrcoef(x, y)[0, 1]
print(f"Correlation between x and y: {corr_xy:.3f}")
print(f"Correlation between x and z: {np.corrcoef(x, z)[0, 1]:.3f}")
print(f"Correlation between y and z: {np.corrcoef(y, z)[0, 1]:.3f}")#Linear Algebra
#Basic Matrix Operations
import numpy as np
# Create matrices
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
B = np.array([[9, 8, 7], [6, 5, 4], [3, 2, 1]])
v = np.array([1, 2, 3])
print("Matrix A:")
print(A)
print("Matrix B:")
print(B)
print("Vector v:", v)
print()
# Matrix multiplication
print("Matrix multiplication A @ B:")
print(A @ B)
print()
print("Matrix multiplication np.dot(A, B):")
print(np.dot(A, B))
print()
# Matrix-vector multiplication
print("Matrix-vector multiplication A @ v:")
print(A @ v)
print()
# Transpose
print("Transpose of A:")
print(A.T)
print()
# Trace (sum of diagonal elements)
print(f"Trace of A: {np.trace(A)}")
print(f"Trace of B: {np.trace(B)}")#Advanced Linear Algebra Operations
import numpy as np
# Create invertible matrix
A = np.array([[2, 1, 1], [1, 3, 2], [1, 0, 0]])
print("Matrix A:")
print(A)
print()
# Determinant
det_A = np.linalg.det(A)
print(f"Determinant det(A): {det_A:.3f}")
print()
# Matrix inverse
if det_A != 0:
A_inv = np.linalg.inv(A)
print("Inverse of A:")
print(A_inv)
print()
# Verify A * A^(-1) = I
identity = A @ A_inv
print("A @ A^(-1) (should be identity matrix):")
print(identity)
print()
# Eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(A)
print("Eigenvalues:")
print(eigenvalues)
print("Eigenvectors:")
print(eigenvectors)
print()
# Matrix rank
rank_A = np.linalg.matrix_rank(A)
print(f"Matrix rank: {rank_A}")
print()
# Singular value decomposition
U, s, Vt = np.linalg.svd(A)
print("Singular value decomposition:")
print(f"U shape: {U.shape}")
print(f"Singular values: {s}")
print(f"Vt shape: {Vt.shape}")#Solving Linear Systems
import numpy as np
# Solve linear system Ax = b
# Example: 2x + y + z = 8
# x + 3y + 2z = 20
# x + 0y + 0z = 2
A = np.array([[2, 1, 1], [1, 3, 2], [1, 0, 0]])
b = np.array([8, 20, 2])
print("Coefficient matrix A:")
print(A)
print("Constant vector b:", b)
print()
# Solve using solve
x = np.linalg.solve(A, b)
print("Solution vector x:", x)
print()
# Verify solution
verification = A @ x
print("Verification A @ x:", verification)
print("Original b:", b)
print("Error:", np.abs(verification - b))
print()
# Least squares solution (for overdetermined systems)
A_over = np.array([[2, 1, 1], [1, 3, 2], [1, 0, 0], [0, 1, 1]])
b_over = np.array([8, 20, 2, 5])
print("Overdetermined system:")
print("A_over:")
print(A_over)
print("b_over:", b_over)
# Least squares solution
x_lstsq = np.linalg.lstsq(A_over, b_over, rcond=None)[0]
print("Least squares solution:", x_lstsq)
# Calculate residual
residual = A_over @ x_lstsq - b_over
print("Residual:", residual)
print(f"Sum of squared residuals: {np.sum(residual**2):.6f}")#Complex Number Operations
#Complex Basics
import numpy as np
# Create complex array
complex_array = np.array([1+2j, 3-4j, 5+0j, 0+6j])
print("Complex array:", complex_array)
print("Data type:", complex_array.dtype)
print()
# Complex properties
print("Real part np.real():", np.real(complex_array))
print("Imaginary part np.imag():", np.imag(complex_array))
print("Magnitude np.abs():", np.abs(complex_array))
print("Phase angle np.angle():", np.angle(complex_array))
print("Conjugate np.conj():", np.conj(complex_array))
print()
# Complex operations
z1 = 3 + 4j
z2 = 1 - 2j
print(f"z1 = {z1}")
print(f"z2 = {z2}")
print(f"z1 + z2 = {z1 + z2}")
print(f"z1 * z2 = {z1 * z2}")
print(f"z1 / z2 = {z1 / z2}")
print(f"|z1| = {abs(z1)}")
print(f"z1* = {np.conj(z1)}")#Broadcasting
#Understanding Broadcasting
import numpy as np
# Operations with different shaped arrays
array_2d = np.array([[1, 2, 3], [4, 5, 6]])
array_1d = np.array([10, 20, 30])
scalar = 100
print("2D array:")
print(array_2d)
print(f"Shape: {array_2d.shape}")
print()
print("1D array:", array_1d)
print(f"Shape: {array_1d.shape}")
print()
print(f"Scalar: {scalar}")
print()
# Broadcasting operations
print("2D array + 1D array:")
result_2d_1d = array_2d + array_1d
print(result_2d_1d)
print(f"Result shape: {result_2d_1d.shape}")
print()
print("2D array + scalar:")
result_2d_scalar = array_2d + scalar
print(result_2d_scalar)
print(f"Result shape: {result_2d_scalar.shape}")
print()
# Column vector broadcasting
column_vector = np.array([[1], [10]])
print("Column vector:")
print(column_vector)
print(f"Shape: {column_vector.shape}")
print("2D array * column vector:")
result_col = array_2d * column_vector
print(result_col)
print(f"Result shape: {result_col.shape}")#Practical Examples
#Example 1: Signal Processing
import numpy as np
# Generate signal
t = np.linspace(0, 1, 1000) # Time axis
frequency1 = 5 # 5Hz
frequency2 = 10 # 10Hz
# Combined signal
signal = np.sin(2 * np.pi * frequency1 * t) + 0.5 * np.sin(2 * np.pi * frequency2 * t)
noise = 0.2 * np.random.normal(0, 1, len(t))
noisy_signal = signal + noise
print(f"Signal length: {len(signal)}")
print(f"Signal statistics:")
print(f" Mean: {np.mean(noisy_signal):.4f}")
print(f" Standard deviation: {np.std(noisy_signal):.4f}")
print(f" Maximum: {np.max(noisy_signal):.4f}")
print(f" Minimum: {np.min(noisy_signal):.4f}")
print()
# Simple filtering (moving average)
window_size = 10
filtered_signal = np.convolve(noisy_signal, np.ones(window_size)/window_size, mode='same')
print("After filtering statistics:")
print(f" Standard deviation: {np.std(filtered_signal):.4f}")
print(f" SNR improvement: {np.std(noisy_signal)/np.std(filtered_signal):.2f}x")
print()
# Frequency domain analysis (simple power spectrum)
fft_result = np.fft.fft(signal)
power_spectrum = np.abs(fft_result)**2
frequencies = np.fft.fftfreq(len(signal), t[1] - t[0])
# Find main frequency components
positive_freq_idx = frequencies > 0
positive_freqs = frequencies[positive_freq_idx]
positive_power = power_spectrum[positive_freq_idx]
# Find frequency with maximum power
max_power_idx = np.argmax(positive_power)
dominant_frequency = positive_freqs[max_power_idx]
print(f"Dominant frequency component: {dominant_frequency:.1f} Hz")#Example 2: Financial Data Analysis
import numpy as np
# Simulate stock price data
np.random.seed(42)
days = 252 # Trading days per year
initial_price = 100
daily_returns = np.random.normal(0.001, 0.02, days) # Daily returns
prices = initial_price * np.cumprod(1 + daily_returns)
print(f"Stock price analysis ({days} trading days):")
print(f" Initial price: ${initial_price:.2f}")
print(f" Final price: ${prices[-1]:.2f}")
print(f" Total return: {(prices[-1]/initial_price - 1)*100:.2f}%")
print()
# Return statistics
print("Daily return statistics:")
print(f" Average daily return: {np.mean(daily_returns)*100:.3f}%")
print(f" Return std: {np.std(daily_returns)*100:.3f}%")
print(f" Largest gain: {np.max(daily_returns)*100:.2f}%")
print(f" Largest loss: {np.min(daily_returns)*100:.2f}%")
print()
# Risk metrics
# 1. Annualized volatility
annual_volatility = np.std(daily_returns) * np.sqrt(252)
print(f"Annualized volatility: {annual_volatility*100:.2f}%")
# 2. Maximum drawdown
running_max = np.maximum.accumulate(prices)
drawdown = (prices - running_max) / running_max
max_drawdown = np.min(drawdown)
print(f"Maximum drawdown: {max_drawdown*100:.2f}%")
# 3. Sharpe ratio (assuming 2% risk-free rate)
risk_free_rate = 0.02
excess_returns = daily_returns - risk_free_rate/252
sharpe_ratio = np.mean(excess_returns) / np.std(excess_returns) * np.sqrt(252)
print(f"Sharpe ratio: {sharpe_ratio:.3f}")
print()
# Moving averages
window_short = 20 # 20-day MA
window_long = 50 # 50-day MA
ma_short = np.convolve(prices, np.ones(window_short)/window_short, mode='valid')
ma_long = np.convolve(prices, np.ones(window_long)/window_long, mode='valid')
print(f"Technical indicators:")
print(f" Current price: ${prices[-1]:.2f}")
print(f" 20-day MA: ${ma_short[-1]:.2f}")
print(f" 50-day MA: ${ma_long[-1]:.2f}")
# Simple trading signal
if len(ma_short) > 0 and len(ma_long) > 0:
if ma_short[-1] > ma_long[-1]:
signal = "Buy signal"
else:
signal = "Sell signal"
print(f" Trading signal: {signal}")#Performance Optimization
#Vectorization vs Loops
import numpy as np
import time
# Compare vectorized and loop performance
size = 1000000
array1 = np.random.random(size)
array2 = np.random.random(size)
# Using loop
start_time = time.time()
result_loop = np.zeros(size)
for i in range(size):
result_loop[i] = array1[i] * array2[i] + np.sin(array1[i])
loop_time = time.time() - start_time
# Using vectorization
start_time = time.time()
result_vectorized = array1 * array2 + np.sin(array1)
vectorized_time = time.time() - start_time
print(f"Array size: {size:,}")
print(f"Loop method time: {loop_time:.4f} seconds")
print(f"Vectorized method time: {vectorized_time:.4f} seconds")
print(f"Performance improvement: {loop_time/vectorized_time:.1f}x")
print(f"Results equal: {np.allclose(result_loop, result_vectorized)}")#Chapter Summary
In this chapter, we learned:
- Basic arithmetic operations and mathematical functions
- Trigonometric and inverse trigonometric functions
- Statistical functions and descriptive statistics
- Linear algebra operations
- Complex number operations
- Broadcasting mechanism
- Practical examples (signal processing, financial analysis)
- Performance optimization techniques
#Next Steps
In the next chapter, we'll learn NumPy's advanced features, including random number generation, polynomials, Fourier transforms, and more.
#Exercises
- Create a function to calculate standardized scores (z-score) of an array
- Implement a simple moving average filter
- Write a function to calculate cosine similarity between two vectors
- Use linear algebra to solve a multiple linear regression problem
- Implement a simple image convolution function