Python Basics

Section Detail

The Python Ecosystem and Philosophy

The Philosophy of Python

Python is a high-level, interpreted programming language known for its readability and versatility. Often described as “executable pseudocode,” Python’s syntax allows developers to express concepts in fewer lines of code than might be possible in languages such as C++ or Java.

The Zen of Python

Python’s design is guided by a set of principles known as The Zen of Python (PEP 20). These principles emphasize simplicity, clarity, and beauty.

Beautiful is better than ugly.
Explicit is better than implicit.
Simple is better than complex.
Readability counts.

In the context of scientific computing, these principles are invaluable. Scientists and researchers need a language that stays out of their way, allowing them to focus on algorithms and data rather than memory management or boilerplate syntax.

Why Python for Science?

Python has become the de facto standard for data science, machine learning, and scientific research for several key reasons:

1. The “Glue Language” Property

Python is exceptionally good at interfaced with other languages. Most performance-critical scientific libraries (like NumPy, SciPy, and TensorFlow) are actually written in C, C++, or Fortran for speed, with Python providing a high-level “glue” interface that is easy to use.

2. Rich Library Ecosystem

Instead of reinventing the wheel, Python users leverage a massive ecosystem of specialized libraries:

NumPy: The foundation for numerical computing.
Pandas: Essential for data manipulation and analysis.
Matplotlib/Seaborn: For data visualization.
SciPy: For advanced scientific calculations (integration, optimization).
SymPy: For symbolic mathematics.
Scikit-Learn: The standard for classical machine learning.

3. Community and Documentation

The scientific Python community (SciPy stack) is one of the most robust in the world, ensuring that libraries are well-maintained and that help is always available.

The Execution Model

Python uses a bytecode-interpreted execution model. While slower than compiled languages for raw loop execution, its efficient C-based library backends often make scientific Python code nearly as fast as hand-coded C for vectorizable operations.

Code

skinparam componentStyle rectangle

package "Development" {
component "Python Source (.py)" as SRC
}

node "Execution Environment" {
component "CPython Interpreter" as INT
component "Bytecode Compiler" as COMP
component "PVM (Python Virtual Machine)" as PVM
}

node "Native Code" {
component "NumPy Core (C/Fortran)" as NP
component "System Libs" as SYS
}

SRC --> COMP
COMP --> INT : "Bytecode (.pyc)"
INT --> PVM
PVM <-> NP : "Vectorized Operations"
PVM <-> SYS

Getting Started

Let’s look at a simple Python script that demonstrates its clean syntax.

Interactive Lab


import math

def calculate_circle_area(radius):
  """Simple function to calculate circle area."""
  return math.pi * (radius ** 2)

# Calculate areas for a range of radii
radii = [1, 2, 3, 4, 5]
areas = [calculate_circle_area(r) for r in radii]

for r, a in zip(radii, areas):
  print(f"Radius: {r}, Area: {a:.2f}")

Expected output

Radius: 1, Area: 3.14
Radius: 2, Area: 12.57
Radius: 3, Area: 28.27
Radius: 4, Area: 50.27
Radius: 5, Area: 78.54

python

1

2import math

3

4def calculate_circle_area(radius):

5 """Simple function to calculate circle area."""

6 return math.pi * (radius ** 2)

7

8# Calculate areas for a range of radii

9radii = [1, 2, 3, 4, 5]

10areas = [calculate_circle_area(r) for r in radii]

11

12for r, a in zip(radii, areas):

13 print(f"Radius: {r}, Area: {a:.2f}")

14

In the following modules, we will dive deep into how to leverage this simplicity for complex scientific tasks.

NumPy Deep Dive

Section Detail

NumPy: The Foundation of Numerical Computing

The Importance of NumPy

NumPy (Numerical Python) is the fundamental package for scientific computing in Python. It provides a high-performance multidimensional array object and tools for working with these arrays. If you are doing any form of data science or scientific computing in Python, NumPy is the engine under the hood.

Why not use Python Lists?

Python lists are incredibly flexible—they can hold elements of different types, and they grow dynamically. However, this flexibility comes at a significant performance cost.

The Problem with Python Lists:

Memory Overhead: Each element in a Python list is a full-fledged object. A list of integers doesn’t just store the numbers; it stores pointers to integer objects, which each contain type information and reference counts.
Lack of Locality: Because list elements are pointers, they can be scattered across memory. This prevents hardware-level optimizations like CPU caching and pre-fetching.
Looping Speed: Iterating over a Python list in a loop is slow because the interpreter must check the type and perform dispatching for every single operation.

The NumPy Solution: `ndarray`

The core of NumPy is the ndarray (n-dimensional array) object.

Contiguous Memory: ndarrays store data in a single, contiguous block of memory.
Homogeneous Types: Every element in a NumPy array must be of the same type (e.g., all 64-bit floats).
Vectorized Operations: Operations on NumPy arrays are performed by compiled C/Fortran code, which can process entire arrays at once without Python loop overhead.

Creating Arrays

NumPy provides multiple ways to initialize arrays.

Interactive Lab


import numpy as np

# From a list
a = np.array([1, 2, 3, 4, 5])
print("From list:", a)

# Filled with zeros or ones
zeros = np.zeros((3, 3))
ones = np.ones((2, 4))
print("Zeros:\n", zeros)

# Ranges
arange = np.arange(0, 10, 2)
linspace = np.linspace(0, 1, 5) # 5 points between 0 and 1
print("Arange:", arange)
print("Linspace:", linspace)

# Random numbers
rand = np.random.rand(2, 2)
print("Random:\n", rand)

Expected output

From list: [1 2 3 4 5]
Zeros:
 [[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]
Arange: [0 2 4 6 8]
Linspace: [0.   0.25 0.5  0.75 1.  ]
Random:
 [[0.45 0.12]
 [0.89 0.34]]

python

1

2import numpy as np

3

4# From a list

5a = np.array([1, 2, 3, 4, 5])

6print("From list:", a)

7

8# Filled with zeros or ones

9zeros = np.zeros((3, 3))

10ones = np.ones((2, 4))

11print("Zeros:\n", zeros)

12

13# Ranges

14arange = np.arange(0, 10, 2)

15linspace = np.linspace(0, 1, 5) # 5 points between 0 and 1

16print("Arange:", arange)

17print("Linspace:", linspace)

18

19# Random numbers

20rand = np.random.rand(2, 2)

21print("Random:\n", rand)

22

Array Attributes

Every ndarray has properties that describe its shape, size, and data type.

shape: A tuple indicating the size of each dimension (e.g., (rows, cols)).
ndim: The number of dimensions (axes).
size: The total number of elements.
dtype: The data type of the elements (e.g., int32, float64).

Interactive Lab


import numpy as np

arr = np.array([[1, 2, 3], [4, 5, 6]])

print(f"Shape: {arr.shape}")
print(f"Dimensions: {arr.ndim}")
print(f"Size: {arr.size}")
print(f"Data Type: {arr.dtype}")

Expected output

Shape: (2, 3)
Dimensions: 2
Size: 6
Data Type: int64

python

1

2import numpy as np

3

4arr = np.array([[1, 2, 3], [4, 5, 6]])

5

6print(f"Shape: {arr.shape}")

7print(f"Dimensions: {arr.ndim}")

8print(f"Size: {arr.size}")

9print(f"Data Type: {arr.dtype}")

10

Vectorization: The Secret Sauce

Vectorization is the process of replacing explicit Python loops with array expressions. This is the primary way to achieve high performance in NumPy.

Consider adding two large vectors:

Interactive Lab


import numpy as np
import time

size = 1000000
a = list(range(size))
b = list(range(size))

# Python Loop approach
start = time.time()
c = [a[i] + b[i] for i in range(size)]
end = time.time()
print(f"Python list addition time: {end - start:.4f}s")

# NumPy approach
na = np.arange(size)
nb = np.arange(size)
start = time.time()
nc = na + nb
end = time.time()
print(f"NumPy vector addition time: {end - start:.4f}s")

Expected output

Python list addition time: 0.0824s
NumPy vector addition time: 0.0012s

python

1

2import numpy as np

3import time

4

5size = 1000000

6a = list(range(size))

7b = list(range(size))

8

9# Python Loop approach

10start = time.time()

11c = [a[i] + b[i] for i in range(size)]

12end = time.time()

13print(f"Python list addition time: {end - start:.4f}s")

14

15# NumPy approach

16na = np.arange(size)

17nb = np.arange(size)

18start = time.time()

19nc = na + nb

20end = time.time()

21print(f"NumPy vector addition time: {end - start:.4f}s")

22

As you can see, the NumPy version is orders of magnitude faster. It’s not just “shorter code”—it’s fundamentally different execution.

Memory Layout and Slicing

NumPy slicing is unique because it creates views of the data rather than copies. This is extremely efficient for large datasets but requires caution: changing a slice changes the original array.

Interactive Lab


import numpy as np

arr = np.array([0, 1, 2, 3, 4, 5])
s = arr[1:4]
s[0] = 99

print("Original array after modifying slice:", arr)

Expected output

Original array after modifying slice: [ 0 99  2  3  4  5]

python

1

2import numpy as np

3

4arr = np.array([0, 1, 2, 3, 4, 5])

5s = arr[1:4]

6s[0] = 99

7

8print("Original array after modifying slice:", arr)

9

In the next module, we will explore advanced indexing and multi-dimensional array manipulations.

Section Detail

Advanced NumPy: Indexing, Slicing, and Reshaping

Multi-Dimensional Indexing

In the previous module, we touched on basic slicing. However, NumPy’s true power comes from its ability to manipulate multi-dimensional arrays with surgical precision.

Comma-Separated Indexing

Unlike Python lists, where you access nested elements using list[i][j], NumPy allows you to use a single set of brackets with comma-separated indices: arr[i, j]. This is not just syntactic sugar; it is more efficient and allows for complex expressions.

Interactive Lab


import numpy as np

arr = np.array([[10, 20, 30], [40, 50, 60], [70, 80, 90]])

print("Element at [1, 2]:", arr[1, 2]) # Row 1, Col 2
print("All rows, first column:", arr[:, 0])
print("First two rows, last two columns:\n", arr[:2, 1:])

Expected output

Element at [1, 2]: 60
All rows, first column: [10 40 70]
First two rows, last two columns:
 [[20 30]
 [50 60]]

python

1

2import numpy as np

3

4arr = np.array([[10, 20, 30], [40, 50, 60], [70, 80, 90]])

5

6print("Element at [1, 2]:", arr[1, 2]) # Row 1, Col 2

7print("All rows, first column:", arr[:, 0])

8print("First two rows, last two columns:\n", arr[:2, 1:])

9

Integer Array Indexing (Fancy Indexing)

Fancy indexing is a term for passing arrays of indices to access multiple array elements at once.

Interactive Lab


import numpy as np

a = np.arange(12).reshape((3, 4))
print("Array:\n", a)

# Select elements at [0,0], [1,1], [2,0]
rows = np.array([0, 1, 2])
cols = np.array([0, 1, 0])
print("\nSelected elements:", a[rows, cols])

Expected output

Array:
 [[ 0  1  2  3]
  [ 4  5  6  7]
  [ 8  9 10 11]]

Selected elements: [0 5 8]

python

1

2import numpy as np

3

4a = np.arange(12).reshape((3, 4))

5print("Array:\n", a)

6

7# Select elements at [0,0], [1,1], [2,0]

8rows = np.array([0, 1, 2])

9cols = np.array([0, 1, 0])

10print("\nSelected elements:", a[rows, cols])

11

Boolean Indexing (Masking)

This is perhaps the most useful feature for data processing. You can index an array using another array of booleans of the same shape.

Interactive Lab


import numpy as np

arr = np.array([1, 2, 3, 4, 5, 6])

# Create a mask where elements are greater than 3
mask = arr > 3
print("Mask:", mask)

# Apply the mask
print("Filtered Array:", arr[mask])

# Complex conditions
mask_complex = (arr > 2) & (arr < 6)
print("Filtered (2 < x < 6):", arr[mask_complex])

Expected output

Mask: [False False False  True  True  True]
Filtered Array: [4 5 6]
Filtered (2 < x < 6): [3 4 5]

python

1

2import numpy as np

3

4arr = np.array([1, 2, 3, 4, 5, 6])

5

6# Create a mask where elements are greater than 3

7mask = arr > 3

8print("Mask:", mask)

9

10# Apply the mask

11print("Filtered Array:", arr[mask])

12

13# Complex conditions

14mask_complex = (arr > 2) & (arr < 6)

15print("Filtered (2 < x < 6):", arr[mask_complex])

16

Shape Manipulation

Frequently, data arrives in a shape that doesn’t match the input requirements of your algorithm (e.g., deep learning models often require a batch dimension).

Reshape

reshape() gives a new shape to an array without changing its data.

Interactive Lab


import numpy as np

# 1D array of 0-11
a = np.arange(12)
print("1D array:", a)

# Reshape to 3x4
a_2d = a.reshape((3, 4))
print("\n3x4 array:\n", a_2d)

# Flatten back to 1D
print("\nFlattened:", a_2d.flatten())

Expected output

1D array: [ 0  1  2  3  4  5  6  7  8  9 10 11]

3x4 array:
 [[ 0  1  2  3]
  [ 4  5  6  7]
  [ 8  9 10 11]]

Flattened: [ 0  1  2  3  4  5  6  7  8  9 10 11]

python

1

2import numpy as np

3

4# 1D array of 0-11

5a = np.arange(12)

6print("1D array:", a)

7

8# Reshape to 3x4

9a_2d = a.reshape((3, 4))

10print("\n3x4 array:\n", a_2d)

11

12# Flatten back to 1D

13print("\nFlattened:", a_2d.flatten())

14

Stacking and Splitting

Combining multiple arrays into one or splitting one into many.

Interactive Lab


import numpy as np

x = np.array([1, 2, 3])
y = np.array([4, 5, 6])

# Vertical Stack
print("V-Stack:\n", np.vstack((x, y)))

# Horizontal Stack
print("H-Stack:", np.hstack((x, y)))

Expected output

V-Stack:
 [[1 2 3]
 [4 5 6]]
H-Stack: [1 2 3 4 5 6]

python

1

2import numpy as np

3

4x = np.array([1, 2, 3])

5y = np.array([4, 5, 6])

6

7# Vertical Stack

8print("V-Stack:\n", np.vstack((x, y)))

9

10# Horizontal Stack

11print("H-Stack:", np.hstack((x, y)))

12

Broadcasting: A Preview

If you try to add a scalar to an array, NumPy implicitly expands the scalar to match the array’s shape. This is called broadcasting.

Interactive Lab


import numpy as np

arr = np.array([[1, 2, 3], [4, 5, 6]])
print("Array + 10:\n", arr + 10)

Expected output

Array + 10:
 [[11 12 13]
 [14 15 16]]

python

1

2import numpy as np

3

4arr = np.array([[1, 2, 3], [4, 5, 6]])

5print("Array + 10:\n", arr + 10)

6

In the next module, we will explore the formal rules of broadcasting and how it allows for high-performance operations without unnecessary memory replication.

Section Detail

The Rules of NumPy Broadcasting

Efficient Operations with Broadcasting

Broadcasting is a powerful mechanism that allows NumPy to perform arithmetic operations on arrays of different shapes. Usually, for element-wise operations, the arrays must be the same size. Broadcasting allows “stretching” the smaller array across the larger one to make them compatible.

The Mental Model

Imagine you have a $3 \times 3$ matrix and you want to add a $1 \times 3$ vector to every row. Using a loop is slow. Instead, NumPy “broadcasts” the vector.

Code

skinparam componentStyle rectangle

component "Row 1" as M1
component "Row 2" as M2
component "Row 3" as M3

component "Vector (1x3)" as V1

component "Result Row 1" as R1
component "Result Row 2" as R2
component "Result Row 3" as R3

M1 --> R1 : "+ V1"
M2 --> R2 : "+ V1"
M3 --> R3 : "+ V1"

note right of V1: "Vector values: [10, 20, 30]\nBroadcast across each matrix row"
note left of R1: "[11, 22, 33]"
note left of R2: "[14, 25, 36]"
note left of R3: "[17, 28, 39]"

The Four Rules of Broadcasting

When operating on two arrays, NumPy compares their shapes element-wise. It starts with the trailing (rightmost) dimensions and works its way left. Two dimensions are compatible when:

They are equal.
One of them is 1.

Example 1: Scalar and Array

A scalar has effectively an infinite number of dimensions of size 1.

Array A: (3, 3)
Array B: (1, ) (Scalar) Result: (3, 3)

Example 2: Vector and Matrix

Array A: (3, 3)
Array B: (3, ) -> broadcast becomes (1, 3) Comparing trailing dimensions: 3 and 3 (Equal). Comparing next: 3 and 1 (One of them is 1). Result: (3, 3)

Example 3: Incompatible Shapes

Array A: (3, 3)
Array B: (2, ) Trailing dimensions: 3 and 2 (Not equal, neither is 1). Result: ValueError!

Broadcasting in Action

Let’s test these rules with code.

Interactive Lab


import numpy as np

# Rule 1: Matrix + Scalar
m = np.ones((3, 3))
print("Matrix + 5:\n", m + 5)

# Rule 2: Matrix + Row Vector
v_row = np.array([1, 2, 3])
print("\nMatrix + Row Vector (v.shape=(3,)):\n", m + v_row)

# Rule 3: Matrix + Column Vector
v_col = np.array([[1], [2], [3]]) # Shape (3, 1)
print("\nMatrix + Column Vector (v.shape=(3, 1)):\n", m + v_col)

Expected output

Matrix + 5:
 [[6. 6. 6.]
 [6. 6. 6.]
 [6. 6. 6.]]

Matrix + Row Vector (v.shape=(3,)):
 [[2. 3. 4.]
 [2. 3. 4.]
 [2. 3. 4.]]

Matrix + Column Vector (v.shape=(3, 1)):
 [[2. 2. 2.]
 [3. 3. 3.]
 [4. 4. 4.]]

python

1

2import numpy as np

3

4# Rule 1: Matrix + Scalar

5m = np.ones((3, 3))

6print("Matrix + 5:\n", m + 5)

7

8# Rule 2: Matrix + Row Vector

9v_row = np.array([1, 2, 3])

10print("\nMatrix + Row Vector (v.shape=(3,)):\n", m + v_row)

11

12# Rule 3: Matrix + Column Vector

13v_col = np.array([[1], [2], [3]]) # Shape (3, 1)

14print("\nMatrix + Column Vector (v.shape=(3, 1)):\n", m + v_col)

15

Centering an Array

Broadcasting is vital for data preprocessing. For instance, to “center” a dataset (subtract the mean of each feature), you can use broadcasting.

Interactive Lab


import numpy as np

# Random data: 10 samples, 3 features
data = np.random.randn(10, 3)

# Calculate mean along the columns (feature mean)
# Shape will be (3,)
mean = data.mean(axis=0)

# Subtract mean from data
# (10, 3) - (3,) -> (3,) is broadcast to (10, 3)
centered_data = data - mean

print("Original shape:", data.shape)
print("Mean shape:", mean.shape)
print("Centered mean (should be ~0):", centered_data.mean(axis=0))

Expected output

Original shape: (10, 3)
Mean shape: (3,)
Centered mean (should be ~0): [ 0.e+00 -1.e-17 -3.e-17]

python

1

2import numpy as np

3

4# Random data: 10 samples, 3 features

5data = np.random.randn(10, 3)

6

7# Calculate mean along the columns (feature mean)

8# Shape will be (3,)

9mean = data.mean(axis=0)

10

11# Subtract mean from data

12# (10, 3) - (3,) -> (3,) is broadcast to (10, 3)

13centered_data = data - mean

14

15print("Original shape:", data.shape)

16print("Mean shape:", mean.shape)

17print("Centered mean (should be ~0):", centered_data.mean(axis=0))

18

Performance Benefits

Broadcasting is computationally efficient because it does not actually replicate the data in memory. The “stretching” is conceptual and handled at the C-level, minimizing memory bandwidth usage.

In the next module, we will apply these techniques to Linear Algebra operations like matrix multiplication and decompositions.

SymPy Deep Dive

Section Detail

SymPy: Symbolic Mathematics in Python

Introduction to Symbolic Computation

Most numerical libraries (like NumPy) work with floating-point numbers. While fast, floating-point numbers are always approximations. For example, $\sqrt{2}$ in NumPy is represented as 1.4142135623730951. In symbolic mathematics, we want to keep it as $\sqrt{2}$ to maintain exactness throughout our derivations.

SymPy is a Python library for symbolic mathematics. It aims to become a full-featured Computer Algebra System (CAS) while keeping the code as simple as possible in order to be comprehensible and easily extensible.

Symbols and Expressions

In SymPy, we must explicitly define variables as symbolic objects.

Interactive Lab


import sympy as sp

# Define symbols
x, y = sp.symbols('x y')

# Create an expression
expr = x**2 + 2*x + 1

print("Expression:", expr)
print("Substituted (x=2):", expr.subs(x, 2))

Expected output

Expression: x**2 + 2*x + 1
Substituted (x=2): 9

python

1

2import sympy as sp

3

4# Define symbols

5x, y = sp.symbols('x y')

6

7# Create an expression

8expr = x**2 + 2*x + 1

9

10print("Expression:", expr)

11print("Substituted (x=2):", expr.subs(x, 2))

12

Structural Manipulations

One of the most powerful features of SymPy is the ability to simplify, expand, and factor expressions automatically.

Expansion and Factoring

Interactive Lab


import sympy as sp
x = sp.symbols('x')

poly = (x + 1)**3
expanded = sp.expand(poly)
factored = sp.factor(expanded)

print(f"Original: {poly}")
print(f"Expanded: {expanded}")
print(f"Factored: {factored}")

Expected output

Original: (x + 1)**3
Expanded: x**3 + 3*x**2 + 3*x + 1
Factored: (x + 1)**3

python

1

2import sympy as sp

3x = sp.symbols('x')

4

5poly = (x + 1)**3

6expanded = sp.expand(poly)

7factored = sp.factor(expanded)

8

9print(f"Original: {poly}")

10print(f"Expanded: {expanded}")

11print(f"Factored: {factored}")

12

Simplification

SymPy has a general-purpose simplify() function that attempts to find the most compact form of an expression.

Interactive Lab


import sympy as sp
x = sp.symbols('x')

expr = (x**2 + x) / x
simplified = sp.simplify(expr)

print(f"Expression: {expr}")
print(f"Simplified: {simplified}")

# Trigonometric simplification
trig_expr = sp.sin(x)**2 + sp.cos(x)**2
print(f"sin^2 + cos^2 = {sp.simplify(trig_expr)}")

Expected output

Expression: (x**2 + x)/x
Simplified: x + 1
sin^2 + cos^2 = 1

python

1

2import sympy as sp

3x = sp.symbols('x')

4

5expr = (x**2 + x) / x

6simplified = sp.simplify(expr)

7

8print(f"Expression: {expr}")

9print(f"Simplified: {simplified}")

10

11# Trigonometric simplification

12trig_expr = sp.sin(x)**2 + sp.cos(x)**2

13print(f"sin^2 + cos^2 = {sp.simplify(trig_expr)}")

14

Calculus with SymPy

SymPy can perform various calculus operations exactly.

Differentiation

Interactive Lab


import sympy as sp
x = sp.symbols('x')

f = sp.sin(x) * sp.exp(x)
derivative = sp.diff(f, x)

print(f"f(x) = {f}")
print(f"f'(x) = {derivative}")

Expected output

f(x) = exp(x)*sin(x)
f'(x) = exp(x)*sin(x) + exp(x)*cos(x)

python

1

2import sympy as sp

3x = sp.symbols('x')

4

5f = sp.sin(x) * sp.exp(x)

6derivative = sp.diff(f, x)

7

8print(f"f(x) = {f}")

9print(f"f'(x) = {derivative}")

10

Integration

SymPy can handle both definite and indefinite integrals.

Interactive Lab


import sympy as sp
x = sp.symbols('x')

# Indefinite integral
indef = sp.integrate(sp.cos(x), x)
print(f"Integral of cos(x): {indef}")

# Definite integral from 0 to pi/2
def_int = sp.integrate(sp.exp(-x), (x, 0, sp.oo)) # oo is infinity
print(f"Integral of e^-x from 0 to infinity: {def_int}")

Expected output

Integral of cos(x): sin(x)
Integral of e^-x from 0 to infinity: 1

python

1

2import sympy as sp

3x = sp.symbols('x')

4

5# Indefinite integral

6indef = sp.integrate(sp.cos(x), x)

7print(f"Integral of cos(x): {indef}")

8

9# Definite integral from 0 to pi/2

10def_int = sp.integrate(sp.exp(-x), (x, 0, sp.oo)) # oo is infinity

11print(f"Integral of e^-x from 0 to infinity: {def_int}")

12

Limits

Interactive Lab


import sympy as sp
x = sp.symbols('x')

limit_val = sp.limit(sp.sin(x)/x, x, 0)
print(f"Limit of sin(x)/x as x -> 0: {limit_val}")

Expected output

Limit of sin(x)/x as x -> 0: 1

python

1

2import sympy as sp

3x = sp.symbols('x')

4

5limit_val = sp.limit(sp.sin(x)/x, x, 0)

6print(f"Limit of sin(x)/x as x -> 0: {limit_val}")

7

Equation Solving

SymPy can solve equations and systems of equations symbolically.

Interactive Lab


import sympy as sp
x, y = sp.symbols('x y')

# Solve x^2 - 1 = 0
solutions = sp.solve(x**2 - 1, x)
print(f"Solutions for x^2 - 1 = 0: {solutions}")

# Solve system of linear equations
# x + y = 5
# x - y = 1
sol_system = sp.solve([x + y - 5, x - y - 1], [x, y])
print(f"System solution: {sol_system}")

Expected output

Solutions for x^2 - 1 = 0: [-1, 1]
System solution: {x: 3, y: 2}

python

1

2import sympy as sp

3x, y = sp.symbols('x y')

4

5# Solve x^2 - 1 = 0

6solutions = sp.solve(x**2 - 1, x)

7print(f"Solutions for x^2 - 1 = 0: {solutions}")

8

9# Solve system of linear equations

10# x + y = 5

11# x - y = 1

12sol_system = sp.solve([x + y - 5, x - y - 1], [x, y])

13print(f"System solution: {sol_system}")

14

Pretty Printing and LaTeX

SymPy can output expressions in various formats, including LaTeX, which is beautiful for documentation and papers.

Interactive Lab


import sympy as sp
x = sp.symbols('x')
expr = sp.sqrt(x**2 + 1) / sp.sin(x)

print("Standard print:", expr)
print("LaTeX:", sp.latex(expr))

Expected output

Standard print: sqrt(x**2 + 1)/sin(x)
LaTeX: \\frac{\\sqrt{x^{2} + 1}}{\\sin{\\left(x \\right)}}

python

1

2import sympy as sp

3x = sp.symbols('x')

4expr = sp.sqrt(x**2 + 1) / sp.sin(x)

5

6print("Standard print:", expr)

7print("LaTeX:", sp.latex(expr))

8

Bridging SymPy and NumPy: `lambdify`

Often you derive an expression symbolically in SymPy but need to evaluate it numerically for thousands of points using NumPy. The lambdify function creates a fast numerical function from a SymPy expression.

Interactive Lab


import sympy as sp
import numpy as np

x = sp.symbols('x')
expr = sp.sin(x)**2

# Convert to a numpy-friendly function
f = sp.lambdify(x, expr, 'numpy')

# Now evaluate on a numpy array
data = np.linspace(0, 10, 5)
results = f(data)
print("Numerical results from symbolic expr:", results)

Expected output

Numerical results from symbolic expr: [0.         0.35816891 0.91953576 0.4265434  0.29562444]

python

1

2import sympy as sp

3import numpy as np

4

5x = sp.symbols('x')

6expr = sp.sin(x)**2

7

8# Convert to a numpy-friendly function

9f = sp.lambdify(x, expr, 'numpy')

10

11# Now evaluate on a numpy array

12data = np.linspace(0, 10, 5)

13results = f(data)

14print("Numerical results from symbolic expr:", results)

15

In the next modules, we will dive deeper into Matrix algebra and structural mechanics applications.

Section Detail

Advanced SymPy: Calculus, Differential Equations, and Optimization

Beyond Simple Algebra

While simplifying polynomials is useful, SymPy’s real power for scientists lies in its ability to solve complex calculus problems and differential equations—tasks that usually require expensive software like Mathematica or Maple.

Taylor Series Expansions

In physics and engineering, we often approximate complex functions using Taylor series. SymPy can generate these to any arbitrary order.

Interactive Lab


import sympy as sp
x = sp.symbols('x')

# Taylor series of sin(x) around 0 up to 10th order
series = sp.series(sp.sin(x), x, 0, 10)
print(f"Taylor series of sin(x):\n{series}")

# Remove the 'O' term for numerical evaluation
func_approx = series.removeO()
print(f"\nPolynomial version: {func_approx}")

Expected output

Taylor series of sin(x):
x - x**3/6 + x**5/120 - x**7/5040 + x**9/362880 + O(x**10)

Polynomial version: x**9/362880 - x**7/5040 + x**5/120 - x**3/6 + x

python

1

2import sympy as sp

3x = sp.symbols('x')

4

5# Taylor series of sin(x) around 0 up to 10th order

6series = sp.series(sp.sin(x), x, 0, 10)

7print(f"Taylor series of sin(x):\n{series}")

8

9# Remove the 'O' term for numerical evaluation

10func_approx = series.removeO()

11print(f"\nPolynomial version: {func_approx}")

12

Solving Ordinary Differential Equations (ODEs)

SymPy can solve many classes of ODEs analytically. This is extremely useful for verifying numerical solvers or finding exact solutions for simple physical models.

Example: The Harmonic Oscillator

The equation is $f''(x) + \omega^2 f(x) = 0$ .

Interactive Lab


import sympy as sp
x = sp.symbols('x')
f = sp.Function('f')
omega = sp.symbols('omega', positive=True)

# Define the ODE
diffeq = sp.Eq(f(x).diff(x, x) + omega**2 * f(x), 0)
print(f"ODE: {diffeq}")

# Solve the ODE
solution = sp.dsolve(diffeq, f(x))
print(f"General Solution: {solution}")

Expected output

ODE: Eq(omega**2*f(x) + Derivative(f(x), (x, 2)), 0)
General Solution: Eq(f(x), C1*sin(omega*x) + C2*cos(omega*x))

python

1

2import sympy as sp

3x = sp.symbols('x')

4f = sp.Function('f')

5omega = sp.symbols('omega', positive=True)

6

7# Define the ODE

8diffeq = sp.Eq(f(x).diff(x, x) + omega**2 * f(x), 0)

9print(f"ODE: {diffeq}")

10

11# Solve the ODE

12solution = sp.dsolve(diffeq, f(x))

13print(f"General Solution: {solution}")

14

Multivariable Calculus

SymPy handles gradients, Jacobians, and Hessians with ease.

Interactive Lab


import sympy as sp
x, y, z = sp.symbols('x y z')

f = x**2 + y**2 + z**2

# Partial derivatives
df_dx = sp.diff(f, x)
print(f"df/dx = {df_dx}")

# Gradient vector
grad = [sp.diff(f, var) for var in (x, y, z)]
print(f"Gradient: {grad}")

Expected output

df/dx = 2*x
Gradient: [2*x, 2*y, 2*z]

python

1

2import sympy as sp

3x, y, z = sp.symbols('x y z')

4

5f = x**2 + y**2 + z**2

6

7# Partial derivatives

8df_dx = sp.diff(f, x)

9print(f"df/dx = {df_dx}")

10

11# Gradient vector

12grad = [sp.diff(f, var) for var in (x, y, z)]

13print(f"Gradient: {grad}")

14

Symbolic Linear Algebra

Calculating the eigenvalues or inverse of a matrix with symbols is a common requirement in theoretical research.

Interactive Lab


import sympy as sp
a, b, c, d = sp.symbols('a b c d')

M = sp.Matrix([[a, b], [c, d]])
print("Matrix M:\n", M)

# Determinant
print("\nDeterminant:", M.det())

# Inverse
print("\nInverse:\n", M.inv())

# Eigenvalues (returns a dict of value: multiplicity)
print("\nEigenvalues:", M.eigenvals())

Expected output

Matrix M:
 Matrix([[a, b], [c, d]])

Determinant: a*d - b*c

Inverse:
 Matrix([[d/(a*d - b*c), -b/(a*d - b*c)], [-c/(a*d - b*c), a/(a*d - b*c)]])

Eigenvalues: {a/2 + d/2 - sqrt(a**2 - 2*a*d + 4*b*c + d**2)/2: 1, a/2 + d/2 + sqrt(a**2 - 2*a*d + 4*b*c + d**2)/2: 1}

python

1

2import sympy as sp

3a, b, c, d = sp.symbols('a b c d')

4

5M = sp.Matrix([[a, b], [c, d]])

6print("Matrix M:\n", M)

7

8# Determinant

9print("\nDeterminant:", M.det())

10

11# Inverse

12print("\nInverse:\n", M.inv())

13

14# Eigenvalues (returns a dict of value: multiplicity)

15print("\nEigenvalues:", M.eigenvals())

16

Global Optimization (Symbolic)

If a function is simple enough, we can find its global minimum by solving where the derivative is zero.

Interactive Lab


import sympy as sp
x = sp.symbols('x')

f = x**4 - 4*x**2 + 3

# Find critical points
critical_points = sp.solve(sp.diff(f, x), x)
print(f"Critical points: {critical_points}")

# Determine which is a minimum using second derivative
for point in critical_points:
  second_deriv = sp.diff(f, x, 2).subs(x, point)
  nature = "minimum" if second_deriv > 0 else "maximum"
  print(f"x = {point} is a {nature} (f''={second_deriv})")

Expected output

Critical points: [0, -sqrt(2), sqrt(2)]
x = 0 is a maximum (f''=-8)
x = -sqrt(2) is a minimum (f''=16)
x = sqrt(2) is a minimum (f''=16)

python

1

2import sympy as sp

3x = sp.symbols('x')

4

5f = x**4 - 4*x**2 + 3

6

7# Find critical points

8critical_points = sp.solve(sp.diff(f, x), x)

9print(f"Critical points: {critical_points}")

10

11# Determine which is a minimum using second derivative

12for point in critical_points:

13 second_deriv = sp.diff(f, x, 2).subs(x, point)

14 nature = "minimum" if second_deriv > 0 else "maximum"

15 print(f"x = {point} is a {nature} (f''={second_deriv})")

16

In the next section, we will integrate these symbolic results into numerical workflows using the SciPy stack.

Section Detail

SymPy: Matrices, Eigenvalues, and Linear Algebra

Symbolic Linear Algebra

While NumPy is the king of numerical linear algebra, SymPy allows us to perform linear algebra with variables. This is essential for deriving formulas, checking proofs, and solving systems where some parameters are not yet known.

Creating Symbolic Matrices

In SymPy, matrices are created using the Matrix class. Unlike NumPy arrays, these matrices are mutable by default and can contain any SymPy expression.

Interactive Lab


import sympy as sp
x, y = sp.symbols('x y')

# Create a 2x2 matrix
M = sp.Matrix([[x, y], [y, x]])
print("Symbolic Matrix M:\n")
sp.pprint(M)

# Basic operations
print("\nMatrix Squared:\n")
sp.pprint(M**2)

Expected output

Symbolic Matrix M:

⎡x  y⎤
⎢    ⎥
⎣y  x⎦

Matrix Squared:

⎡ 2    2            ⎤
⎢x  + y      2⋅x⋅y  ⎥
⎢                   ⎥
⎣  2⋅x⋅y     2    2 ⎦
            x  + y

python

1

2import sympy as sp

3x, y = sp.symbols('x y')

4

5# Create a 2x2 matrix

6M = sp.Matrix([[x, y], [y, x]])

7print("Symbolic Matrix M:\n")

8sp.pprint(M)

9

10# Basic operations

11print("\nMatrix Squared:\n")

12sp.pprint(M**2)

13

Determinants and Inversion

Calculating a determinant symbolically is a common task in stability analysis and geometry.

Interactive Lab


import sympy as sp
a, b, c, d = sp.symbols('a b c d')

# The general 2x2 matrix
A = sp.Matrix([[a, b], [c, d]])

det_A = A.det()
print(f"Determinant of A: {det_A}")

inv_A = A.inv()
print("\nInverse of A:")
sp.pprint(inv_A)

Expected output

Determinant of A: a*d - b*c

Inverse of A:
⎡   d          -b   ⎤
⎢───────     ───────⎥
⎢ad - bc     ad - bc⎥
⎢                   ⎥
⎢  -c          a    ⎥
⎢───────     ───────⎥
⎣ad - bc     ad - bc⎦

python

1

2import sympy as sp

3a, b, c, d = sp.symbols('a b c d')

4

5# The general 2x2 matrix

6A = sp.Matrix([[a, b], [c, d]])

7

8det_A = A.det()

9print(f"Determinant of A: {det_A}")

10

11inv_A = A.inv()

12print("\nInverse of A:")

13sp.pprint(inv_A)

14

Eigenvalues and Eigenvectors

One of the most powerful features of SymPy is calculating exact eigenvalues. Numerical methods might lose precision for nearly singular matrices or complex clusters; SymPy finds the exact roots of the characteristic polynomial.

Interactive Lab


import sympy as sp
lam = sp.symbols('lambda')
M = sp.Matrix([[3, -2], [4, -1]])

# Characteristic polynomial
poly = M.charpoly(lam)
print(f"Characteristic Polynomial: {poly.as_expr()}")

# Eigenvalues
eigenvals = M.eigenvals()
print(f"Eigenvalues: {eigenvals}") # Returns {value: multiplicity}

# Eigenvectors
eigenvects = M.eigenvects()
print("\nEigenvectors (Value, Multiplicity, Basis):")
for v in eigenvects:
  print(v)

Expected output

Characteristic Polynomial: lambda**2 - 2*lambda + 5
Eigenvalues: {1 - 2*I: 1, 1 + 2*I: 1}

Eigenvectors (Value, Multiplicity, Basis):
(1 - 2*I, 1, [Matrix([
[1/2 + I/2],
[        1]])])
(1 + 2*I, 1, [Matrix([
[1/2 - I/2],
[        1]])])

python

1

2import sympy as sp

3lam = sp.symbols('lambda')

4M = sp.Matrix([[3, -2], [4, -1]])

5

6# Characteristic polynomial

7poly = M.charpoly(lam)

8print(f"Characteristic Polynomial: {poly.as_expr()}")

9

10# Eigenvalues

11eigenvals = M.eigenvals()

12print(f"Eigenvalues: {eigenvals}") # Returns {value: multiplicity}

13

14# Eigenvectors

15eigenvects = M.eigenvects()

16print("\nEigenvectors (Value, Multiplicity, Basis):")

17for v in eigenvects:

18 print(v)

19

Matrix Decompositions

SymPy supports several decompositions, including LU, QR, and Diagonalization.

Interactive Lab


import sympy as sp
M = sp.Matrix([[1, 2], [2, 1]])

# Diagonalization: M = P * D * P^-1
P, D = M.diagonalize()

print("Diagonal Matrix D:")
sp.pprint(D)
print("\nTransformation Matrix P:")
sp.pprint(P)

# Verify
assert P * D * P.inv() == M
print("\nVerification successful: P * D * P^-1 == M")

Expected output

Diagonal Matrix D:
⎡-1  0⎤
⎢     ⎥
⎣ 0  3⎦

Transformation Matrix P:
⎡-1  1⎤
⎢     ⎥
⎣ 1  1⎦

Verification successful: P * D * P^-1 == M

python

1

2import sympy as sp

3M = sp.Matrix([[1, 2], [2, 1]])

4

5# Diagonalization: M = P * D * P^-1

6P, D = M.diagonalize()

7

8print("Diagonal Matrix D:")

9sp.pprint(D)

10print("\nTransformation Matrix P:")

11sp.pprint(P)

12

13# Verify

14assert P * D * P.inv() == M

15print("\nVerification successful: P * D * P^-1 == M")

16

Solving Linear Systems ( $Ax = b$ )

You can solve a system of linear equations by passing a matrix and a column vector to the LUsolve method or by using solve_linear_system.

Interactive Lab


import sympy as sp
x1, x2 = sp.symbols('x1 x2')

# System:
# x1 + 2*x2 = 5
# 3*x1 + 4*x2 = 11

A = sp.Matrix([[1, 2], [3, 4]])
b = sp.Matrix([5, 11])

sol = A.LUsolve(b)
print(f"Solution vector: {sol}")

Expected output

Solution vector: Matrix([[1], [2]])

python

1

2import sympy as sp

3x1, x2 = sp.symbols('x1 x2')

4

5# System:

6# x1 + 2*x2 = 5

7# 3*x1 + 4*x2 = 11

8

9A = sp.Matrix([[1, 2], [3, 4]])

10b = sp.Matrix([5, 11])

11

12sol = A.LUsolve(b)

13print(f"Solution vector: {sol}")

14

Applications in Engineering

Symbolic matrices are often used to define Stiffness Matrices in structural analysis or Jacobian Matrices in robot kinematics. Because the variables are preserved, we can calculate the Jacobian once and then substitute specific joint angles thousands of times during a simulation.

Interactive Lab


import sympy as sp
theta = sp.symbols('theta')

# Rotation matrix around Z-axis
R = sp.Matrix([
  [sp.cos(theta), -sp.sin(theta), 0],
  [sp.sin(theta),  sp.cos(theta), 0],
  [0, 0, 1]
])

# Derivative of rotation matrix with respect to theta
dR_dtheta = R.diff(theta)

print("Rotation Matrix Derivative:")
sp.pprint(dR_dtheta)

Expected output

Rotation Matrix Derivative:
⎡-sin(theta)  -cos(theta)  0⎤
⎢                           ⎥
⎢ cos(theta)  -sin(theta)  0⎥
⎢                           ⎥
⎣     0            0       0⎦

python

1

2import sympy as sp

3theta = sp.symbols('theta')

4

5# Rotation matrix around Z-axis

6R = sp.Matrix([

7 [sp.cos(theta), -sp.sin(theta), 0],

8 [sp.sin(theta), sp.cos(theta), 0],

9 [0, 0, 1]

10])

11

12# Derivative of rotation matrix with respect to theta

13dR_dtheta = R.diff(theta)

14

15print("Rotation Matrix Derivative:")

16sp.pprint(dR_dtheta)

17

In the next module, we will conclude our deep dive into the scientific stack by looking at how to bridge these symbolic results back into the numerical world of SciPy.

Pandas Deep Dive

Section Detail

Pandas: High-Performance Data Structures

Why Pandas?

While NumPy provides the computational horsepower for numerical arrays, it lacks the high-level features needed for real-world data analysis, such as:

Handling missing data (NaN).
Working with labeled axes (instead of just integer indices).
Merging and joining datasets (SQL-like operations).
Time-series functionality.

Pandas builds on top of NumPy to provide these features. It is the primary tool for data “munging” or “wrangling”—the process of cleaning and transforming raw data into a format suitable for analysis.

Core Data Structures: Series and DataFrame

There are two primary data structures in Pandas:

Series: A 1D array-like object with an associated index.
DataFrame: A 2D table-like structure with labeled rows and columns.

The Series Object

Think of a Series as a cross between a NumPy array and a Python dictionary.

Interactive Lab


import pandas as pd
import numpy as np

# Creating a Series from a list
data = pd.Series([0.25, 0.5, 0.75, 1.0], index=['a', 'b', 'c', 'd'])
print("Series:\n", data)

# Accessing by label
print("\nValue at index 'b':", data['b'])

# Series behaves like a numpy array
print("\nMean of Series:", data.mean())

Expected output

Series:
a    0.25
b    0.50
c    0.75
d    1.00
dtype: float64

Value at index 'b': 0.5

Mean of Series: 0.625

python

1

2import pandas as pd

3import numpy as np

4

5# Creating a Series from a list

6data = pd.Series([0.25, 0.5, 0.75, 1.0], index=['a', 'b', 'c', 'd'])

7print("Series:\n", data)

8

9# Accessing by label

10print("\nValue at index 'b':", data['b'])

11

12# Series behaves like a numpy array

13print("\nMean of Series:", data.mean())

14

The DataFrame Object

The DataFrame is the most important structure in Pandas. It represents a table of data, similar to a spreadsheet or a SQL table. Each column in a DataFrame is itself a Series.

Interactive Lab


import pandas as pd

data = {
  'Name': ['Alice', 'Bob', 'Charlie', 'David'],
  'Age': [25, 30, 35, 40],
  'City': ['New York', 'London', 'Paris', 'Tokyo']
}

df = pd.DataFrame(data)
print("DataFrame:\n", df)

# Check basic info
print("\nColumns:", df.columns)
print("Index:", df.index)

Expected output

DataFrame:
       Name  Age      City
0     Alice   25  New York
1       Bob   30    London
2   Charlie   35     Paris
3     David   40     Tokyo

Columns: Index(['Name', 'Age', 'City'], dtype='object')
Index: RangeIndex(start=0, stop=4, step=1)

python

1

2import pandas as pd

3

4data = {

5 'Name': ['Alice', 'Bob', 'Charlie', 'David'],

6 'Age': [25, 30, 35, 40],

7 'City': ['New York', 'London', 'Paris', 'Tokyo']

8}

9

10df = pd.DataFrame(data)

11print("DataFrame:\n", df)

12

13# Check basic info

14print("\nColumns:", df.columns)

15print("Index:", df.index)

16

Data Selection and Indexing

Selection in Pandas can be confusing because there are multiple ways to do it.

1. `loc`: Label-based selection

loc is used to select data using the labels of rows and columns.

2. `iloc`: Integer-based selection

iloc is used to select data using the 0-based integer position.

Interactive Lab


import pandas as pd

df = pd.DataFrame({
  'A': [1, 2, 3],
  'B': [4, 5, 6]
}, index=['row1', 'row2', 'row3'])

print("Selection with loc (labels):\n", df.loc['row2', 'A'])
print("\nSelection with iloc (positions):\n", df.iloc[1, 0])

Expected output

Selection with loc (labels):
 2

Selection with iloc (positions):
 2

python

1

2import pandas as pd

3

4df = pd.DataFrame({

5 'A': [1, 2, 3],

6 'B': [4, 5, 6]

7}, index=['row1', 'row2', 'row3'])

8

9print("Selection with loc (labels):\n", df.loc['row2', 'A'])

10print("\nSelection with iloc (positions):\n", df.iloc[1, 0])

11

Importing Data

In practice, you rarely create DataFrames by hand. You load them from files (CSV, Excel, JSON, SQL).

# Reading a CSV file
df = pd.read_csv('data.csv')

# Writing to an Excel file
df.to_excel('output.xlsx')

Handling Missing Data

Real-world data is messy. Pandas uses NaN (Not a Number) to represent missing values and provides robust tools to handle them.

Interactive Lab


import pandas as pd
import numpy as np

df = pd.DataFrame({
  'A': [1, 2, np.nan, 4],
  'B': [5, np.nan, np.nan, 8],
  'C': [1, 2, 3, 4]
})

print("Original DF with NaNs:\n", df)

# Fill NaNs with a value
print("\nFilled NaNs:\n", df.fillna(value=0))

# Drop rows with any NaN
print("\nDropped NaNs:\n", df.dropna())

Expected output

Original DF with NaNs:
      A    B  C
0  1.0  5.0  1
1  2.0  NaN  2
2  NaN  NaN  3
3  4.0  8.0  4

Filled NaNs:
      A    B  C
0  1.0  5.0  1
1  2.0  0.0  2
2  0.0  0.0  3
3  4.0  8.0  4

Dropped NaNs:
      A    B  C
0  1.0  5.0  1
3  4.0  8.0  4

python

1

2import pandas as pd

3import numpy as np

4

5df = pd.DataFrame({

6 'A': [1, 2, np.nan, 4],

7 'B': [5, np.nan, np.nan, 8],

8 'C': [1, 2, 3, 4]

9})

10

11print("Original DF with NaNs:\n", df)

12

13# Fill NaNs with a value

14print("\nFilled NaNs:\n", df.fillna(value=0))

15

16# Drop rows with any NaN

17print("\nDropped NaNs:\n", df.dropna())

18

GroupBy and Aggregation

Similar to GROUP BY in SQL, Pandas allows you to split data into groups and apply functions to each group separately.

Interactive Lab


import pandas as pd

df = pd.DataFrame({
  'Company': ['GOOG', 'GOOG', 'MSFT', 'MSFT', 'FB', 'FB'],
  'Person': ['Sam', 'Charlie', 'Amy', 'Vanessa', 'Carl', 'Sarah'],
  'Sales': [200, 120, 340, 124, 243, 350]
})

by_comp = df.groupby("Company")
print("Mean sales per company:\n", by_comp['Sales'].mean())

Expected output

Mean sales per company:
Company
FB      296.5
GOOG    160.0
MSFT    232.0
Name: Sales, dtype: float64

python

1

2import pandas as pd

3

4df = pd.DataFrame({

5 'Company': ['GOOG', 'GOOG', 'MSFT', 'MSFT', 'FB', 'FB'],

6 'Person': ['Sam', 'Charlie', 'Amy', 'Vanessa', 'Carl', 'Sarah'],

7 'Sales': [200, 120, 340, 124, 243, 350]

8})

9

10by_comp = df.groupby("Company")

11print("Mean sales per company:\n", by_comp['Sales'].mean())

12

In the following modules, we will explore advanced joining techniques and time-series analysis.

Section Detail

Pandas: Data Cleaning and Wrangling

The Reality of Data

In textbooks, data is clean. In reality, data is missing, duplicated, inconsistent, and scattered across multiple files. Data scientists spend roughly 80% of their time cleaning and wrangling data. Pandas is designed to make this “grunt work” efficient.

Duplicates and Noise

Common cleaning tasks include removing duplicate rows and inconsistent strings.

Interactive Lab


import pandas as pd

df = pd.DataFrame({
  'k1': ['one', 'two'] * 3 + ['two'],
  'k2': [1, 1, 2, 3, 3, 4, 4]
})

print("Original DF:\n", df)

# Identify duplicates
print("\nIs Duplicate:\n", df.duplicated())

# Remove duplicates
print("\nCleaned DF:\n", df.drop_duplicates())

Expected output

Original DF:
     k1  k2
 0  one   1
 1  two   1
 2  one   2
 3  two   3
 4  one   3
 5  two   4
 6  two   4

Is Duplicate:
 0    False
 1    False
 2    False
 3    False
 4    False
 5    False
 6     True
 dtype: bool

Cleaned DF:
     k1  k2
 0  one   1
 1  two   1
 2  one   2
 3  two   3
 4  one   3
 5  two   4

python

1

2import pandas as pd

3

4df = pd.DataFrame({

5 'k1': ['one', 'two'] * 3 + ['two'],

6 'k2': [1, 1, 2, 3, 3, 4, 4]

7})

8

9print("Original DF:\n", df)

10

11# Identify duplicates

12print("\nIs Duplicate:\n", df.duplicated())

13

14# Remove duplicates

15print("\nCleaned DF:\n", df.drop_duplicates())

16

Transforming Data: `map` and `apply`

Sometimes you need to apply a custom function to every element or row in a DataFrame.

map(): Used on a Series to substitute values based on a dictionary or function.
apply(): Used on a DataFrame to apply a function along an axis (rows or columns).

Interactive Lab


import pandas as pd

df = pd.DataFrame({
  'food': ['bacon', 'pulled pork', 'bacon', 'Pastrami', 'corned beef'],
  'ounces': [4, 3, 12, 6, 7.5]
})

# Mapping food to animal source
meat_to_animal = {
  'bacon': 'pig',
  'pulled pork': 'pig',
  'pastrami': 'cow',
  'corned beef': 'cow'
}

# Normalize string and map
df['animal'] = df['food'].str.lower().map(meat_to_animal)
print("Transformation results:\n", df)

Expected output

Transformation results:
           food  ounces animal
 0        bacon     4.0    pig
 1  pulled pork     3.0    pig
 2        bacon    12.0    pig
 3     Pastrami     6.0    cow
 4  corned beef     7.5    cow

python

1

2import pandas as pd

3

4df = pd.DataFrame({

5 'food': ['bacon', 'pulled pork', 'bacon', 'Pastrami', 'corned beef'],

6 'ounces': [4, 3, 12, 6, 7.5]

7})

8

9# Mapping food to animal source

10meat_to_animal = {

11 'bacon': 'pig',

12 'pulled pork': 'pig',

13 'pastrami': 'cow',

14 'corned beef': 'cow'

15}

16

17# Normalize string and map

18df['animal'] = df['food'].str.lower().map(meat_to_animal)

19print("Transformation results:\n", df)

20

Combining Datasets: Merge, Join, and Concatenate

In a relational database model, data is split across tables. Pandas allows you to bring them together.

Merge (SQL Join)

Merge connects rows in DataFrames based on one or more keys.

Interactive Lab


import pandas as pd

df1 = pd.DataFrame({'key': ['b', 'b', 'a', 'c', 'a', 'a', 'b'], 'data1': range(7)})
df2 = pd.DataFrame({'key': ['a', 'b', 'd'], 'data2': range(3)})

# Inner Join (default)
merged = pd.merge(df1, df2, on='key')
print("Merged Data (Inner):\n", merged)

# Outer Join
outer = pd.merge(df1, df2, on='key', how='outer')
print("\nMerged Data (Outer):\n", outer)

Expected output

Merged Data (Inner):
   key  data1  data2
 0   b      0      1
 1   b      1      1
 2   b      6      1
 3   a      2      0
 4   a      4      0
 5   a      5      0

Merged Data (Outer):
   key  data1  data2
 0   a    2.0    0.0
 1   a    4.0    0.0
 2   a    5.0    0.0
 3   b    0.0    1.0
 4   b    1.0    1.0
 5   b    6.0    1.0
 6   c    3.0    NaN
 7   d    NaN    2.0

python

1

2import pandas as pd

3

4df1 = pd.DataFrame({'key': ['b', 'b', 'a', 'c', 'a', 'a', 'b'], 'data1': range(7)})

5df2 = pd.DataFrame({'key': ['a', 'b', 'd'], 'data2': range(3)})

6

7# Inner Join (default)

8merged = pd.merge(df1, df2, on='key')

9print("Merged Data (Inner):\n", merged)

10

11# Outer Join

12outer = pd.merge(df1, df2, on='key', how='outer')

13print("\nMerged Data (Outer):\n", outer)

14

Concatenate

Concatenate “stacks” DataFrames on top of each other or side-by-side.

Interactive Lab


import pandas as pd
import numpy as np

s1 = pd.Series([0, 1], index=['a', 'b'])
s2 = pd.Series([2, 3, 4], index=['c', 'd', 'e'])

print("Concatenated Series:\n", pd.concat([s1, s2]))

Expected output

Concatenated Series:
a    0
b    1
c    2
d    3
e    4
dtype: int64

python

1

2import pandas as pd

3import numpy as np

4

5s1 = pd.Series([0, 1], index=['a', 'b'])

6s2 = pd.Series([2, 3, 4], index=['c', 'd', 'e'])

7

8print("Concatenated Series:\n", pd.concat([s1, s2]))

9

Reshaping and Pivoting

Changing the layout of a DataFrame (Long to Wide format or vice versa) is essential for visualization and certain machine learning models.

Interactive Lab


import pandas as pd

df = pd.DataFrame({
  'date': ['2021-01-01', '2021-01-01', '2021-01-02'],
  'variable': ['temp', 'humidity', 'temp'],
  'value': [22.5, 60, 23.0]
})

# Pivot to wide format
pivoted = df.pivot(index='date', columns='variable', values='value')
print("Pivoted (Wide) Data:\n", pivoted)

Expected output

Pivoted (Wide) Data:
 variable    humidity  temp
 date                      
 2021-01-01      60.0  22.5
 2021-01-02       NaN  23.0

python

1

2import pandas as pd

3

4df = pd.DataFrame({

5 'date': ['2021-01-01', '2021-01-01', '2021-01-02'],

6 'variable': ['temp', 'humidity', 'temp'],

7 'value': [22.5, 60, 23.0]

8})

9

10# Pivot to wide format

11pivoted = df.pivot(index='date', columns='variable', values='value')

12print("Pivoted (Wide) Data:\n", pivoted)

13

Advanced String Manipulation

Pandas provides a .str accessor for performing vectorized string operations.

Interactive Lab


import pandas as pd

names = pd.Series([' Alice ', 'bOB', ' Charlie  ', 'd_avid'])

# Strip whitespace and capitalize
clean_names = names.str.strip().str.capitalize()
print("Clean names:\n", clean_names)

# Regex operations
print("\nContains 'a'?", clean_names.str.contains('a'))

Expected output

Clean names:
 0      Alice
 1        Bob
 2    Charlie
 3     D_avid
 dtype: object

Contains 'a'? 0    False
 1    False
 2     True
 3     True
 dtype: bool

python

1

2import pandas as pd

3

4names = pd.Series([' Alice ', 'bOB', ' Charlie ', 'd_avid'])

5

6# Strip whitespace and capitalize

7clean_names = names.str.strip().str.capitalize()

8print("Clean names:\n", clean_names)

9

10# Regex operations

11print("\nContains 'a'?", clean_names.str.contains('a'))

12

In the next module, we’ll dive into Time-Series analysis, one of Pandas’ strongest features.

Section Detail

Pandas: Time Series Analysis

Why Time Series in Pandas?

Pandas was originally developed in a financial context (at AQR Capital Management), which explains why it has world-class support for time series. Temporal data often requires specialized operations like resampling (converting frequency), time zone conversion, and rolling windows.

The `DatetimeIndex`

The key to time series in Pandas is having a DatetimeIndex.

Interactive Lab


import pandas as pd
import numpy as np

# Create a range of dates
dates = pd.date_range('2023-01-01', periods=6, freq='D')
print("Dates Index:\n", dates)

# Create a Series with the date index
ts = pd.Series(np.random.randn(6), index=dates)
print("\nTime Series:\n", ts)

Expected output

Dates Index:
 DatetimeIndex(['2023-01-01', '2023-01-02', '2023-01-03', '2023-01-04',
               '2023-01-05', '2023-01-06'],
              dtype='datetime64[ns]', freq='D')

Time Series:
 2023-01-01   -1.11
 2023-01-02    0.45
 2023-01-03    1.23
 2023-01-04   -0.56
 2023-01-05    0.89
 2023-01-06   -0.12
 Freq: D, dtype: float64

python

1

2import pandas as pd

3import numpy as np

4

5# Create a range of dates

6dates = pd.date_range('2023-01-01', periods=6, freq='D')

7print("Dates Index:\n", dates)

8

9# Create a Series with the date index

10ts = pd.Series(np.random.randn(6), index=dates)

11print("\nTime Series:\n", ts)

12

Resampling: Upsampling and Downsampling

Resampling is the process of changing the frequency of your time series observations.

Downsampling: Aggregating data (e.g., daily to monthly).
Upsampling: Increasing frequency (e.g., monthly to daily), often requiring interpolation.

Interactive Lab


import pandas as pd
import numpy as np

# Hourly data for 3 days
rng = pd.date_range('1/1/2023', periods=72, freq='H')
ts = pd.Series(np.random.randn(len(rng)), index=rng)

# Downsample to Daily frequency and get the mean
daily_summary = ts.resample('D').mean()
print("Daily Mean:\n", daily_summary)

Expected output

Daily Mean:
 2023-01-01   -0.01
 2023-01-02    0.04
 2023-01-03   -0.12
 Freq: D, dtype: float64

python

1

2import pandas as pd

3import numpy as np

4

5# Hourly data for 3 days

6rng = pd.date_range('1/1/2023', periods=72, freq='H')

7ts = pd.Series(np.random.randn(len(rng)), index=rng)

8

9# Downsample to Daily frequency and get the mean

10daily_summary = ts.resample('D').mean()

11print("Daily Mean:\n", daily_summary)

12

Moving Windows (Rolling Operations)

Rolling operations allow you to calculate statistics over a sliding window of time. This is common for smoothing noisy data or calculating moving averages in finance.

Interactive Lab


import pandas as pd
import numpy as np

# Create 100 days of data
ts = pd.Series(np.random.randn(100), index=pd.date_range('1/1/2023', periods=100))

# Calculate 7-day rolling mean
rolling_mean = ts.rolling(window=7).mean()

print("Original Data (first 10):\n", ts.head(10))
print("\n7-day Moving Average (first 10):\n", rolling_mean.head(10))

Expected output

Original Data (first 10):
 2023-01-01    0.12
 2023-01-02   -0.45
 2023-01-03    0.89
 2023-01-04   -0.34
 2023-01-05    1.56
 2023-01-06    0.11
 2023-01-07   -0.78
 2023-01-08    0.23
 2023-01-09    1.11
 2023-01-10   -0.56
 Freq: D, dtype: float64

7-day Moving Average (first 10):
 2023-01-01         NaN
 2023-01-02         NaN
 2023-01-03         NaN
 2023-01-04         NaN
 2023-01-05         NaN
 2023-01-06         NaN
 2023-01-07    0.3157
 2023-01-08    0.3342
 2023-01-09    0.4568
 2023-01-10    0.3456
 Freq: D, dtype: float64

python

1

2import pandas as pd

3import numpy as np

4

5# Create 100 days of data

6ts = pd.Series(np.random.randn(100), index=pd.date_range('1/1/2023', periods=100))

7

8# Calculate 7-day rolling mean

9rolling_mean = ts.rolling(window=7).mean()

10

11print("Original Data (first 10):\n", ts.head(10))

12print("\n7-day Moving Average (first 10):\n", rolling_mean.head(10))

13

Handling Time Zones

Global data often requires reconciling different time zones.

Interactive Lab


import pandas as pd

ts = pd.Series([1, 2, 3], index=pd.date_range('2023-01-01', periods=3, freq='D'))

# Localize to UTC
ts_utc = ts.tz_localize('UTC')
print("UTC Series:\n", ts_utc)

# Convert to US Eastern time
ts_eastern = ts_utc.tz_convert('US/Eastern')
print("\nEastern Time Series:\n", ts_eastern)

Expected output

UTC Series:
 2023-01-01 00:00:00+00:00    1
 2023-01-02 00:00:00+00:00    2
 2023-01-03 00:00:00+00:00    3
 Freq: D, dtype: int64

Eastern Time Series:
 2022-12-31 19:00:00-05:00    1
 2023-01-01 19:00:00-05:00    2
 2023-01-02 19:00:00-05:00    3
 Freq: D, dtype: int64

python

1

2import pandas as pd

3

4ts = pd.Series([1, 2, 3], index=pd.date_range('2023-01-01', periods=3, freq='D'))

5

6# Localize to UTC

7ts_utc = ts.tz_localize('UTC')

8print("UTC Series:\n", ts_utc)

9

10# Convert to US Eastern time

11ts_eastern = ts_utc.tz_convert('US/Eastern')

12print("\nEastern Time Series:\n", ts_eastern)

13

Shifting and Lagging

In time series modeling, you often want to shift data forward or backward in time (e.g., to calculate percentage changes).

Interactive Lab


import pandas as pd

ts = pd.Series([10, 20, 30, 40], index=pd.date_range('2023-01-01', periods=4))

# Shift forward by 1
shifted = ts.shift(1)
print("Shifted (Lag 1):\n", shifted)

# Calculate percentage change
pct_change = (ts - ts.shift(1)) / ts.shift(1)
print("\nPercent Change:\n", pct_change)

Expected output

Shifted (Lag 1):
 2023-01-01     NaN
 2023-01-02    10.0
 2023-01-03    20.0
 2023-01-04    30.0
 Freq: D, dtype: float64

Percent Change:
 2023-01-01     NaN
 2023-01-02    1.00
 2023-01-03    0.50
 2023-01-04    0.33
 Freq: D, dtype: float64

python

1

2import pandas as pd

3

4ts = pd.Series([10, 20, 30, 40], index=pd.date_range('2023-01-01', periods=4))

5

6# Shift forward by 1

7shifted = ts.shift(1)

8print("Shifted (Lag 1):\n", shifted)

9

10# Calculate percentage change

11pct_change = (ts - ts.shift(1)) / ts.shift(1)

12print("\nPercent Change:\n", pct_change)

13

Period Containers

While Timestamp represents a point in time, Period represents a duration (a day, a month, a year).

Interactive Lab


import pandas as pd

p = pd.Period('2023-01', freq='M')
print(f"Period: {p}")
print(f"Next month: {p + 1}")

# Converting from Timestamp to Period
ts = pd.Timestamp('2023-01-15')
period = ts.to_period('M')
print(f"Timestamp to Monthly Period: {period}")

Expected output

Period: 2023-01
Next month: 2023-02
Timestamp to Monthly Period: 2023-01

python

1

2import pandas as pd

3

4p = pd.Period('2023-01', freq='M')

5print(f"Period: {p}")

6print(f"Next month: {p + 1}")

7

8# Converting from Timestamp to Period

9ts = pd.Timestamp('2023-01-15')

10period = ts.to_period('M')

11print(f"Timestamp to Monthly Period: {period}")

12

In the next section, we will integrate these data manipulation skills with Scikit-Learn to build predictive models.

Scikit-Learn Deep Dive

Section Detail

Scikit-Learn: Machine Learning Principles

The Scikit-Learn Ecosystem

Scikit-Learn (frequently abbreviated as sklearn) is the primary library for classical machine learning in Python. Built on top of NumPy, SciPy, and Matplotlib, it provides simple and efficient tools for predictive data analysis.

Key Concepts and the Estimator API

The brilliance of Scikit-Learn lies in its consistent API. Whether you are performing linear regression, support vector machines, or random forests, the workflow is almost identical.

Code

skinparam componentStyle rectangle

package "Data Preparation" {
[Features (X)]
[Target (y)]
}

node "Scikit-Learn Workflow" {
component "Instantiate Model" as MK
component "fit(X, y)" as FIT
component "predict(X_new)" as PRED
}

[Features (X)] --> MK
[Target (y)] --> MK
MK --> FIT : "Training"
FIT --> PRED : "Inference"

The Three Main Steps:

Instantiate: Choose your model class and set hyperparameters.
Fit: Train the model on your data using the fit() method.
Predict: Apply the trained model to new data using predict() or transform().

Data Representation in Scikit-Learn

Scikit-Learn expects data in a specific format:

X (Feature Matrix): A 2D array or DataFrame of shape [n_samples, n_features].
y (Target Vector): A 1D array or Series of length n_samples.

A Simple Example: Linear Regression

Let’s see how we can predict a continuous value using Scikit-Learn.

Interactive Lab


from sklearn.linear_model import LinearRegression
import numpy as np

# 1. Create Data
X = np.array([[1], [2], [3], [4]]) # Feature matrix
y = np.array([2, 4, 6, 8])         # Target vector (y = 2*x)

# 2. Instantiate Model
model = LinearRegression()

# 3. Fit Model
model.fit(X, y)

# 4. Predict
X_new = np.array([[5], [10]])
predictions = model.predict(X_new)

print(f"Predictions for 5 and 10: {predictions}")
print(f"Coefficient (Slope): {model.coef_[0]}")
print(f"Intercept: {model.intercept_}")

Expected output

Predictions for 5 and 10: [10. 20.]
Coefficient (Slope): 2.0
Intercept: 0.0

python

1

2from sklearn.linear_model import LinearRegression

3import numpy as np

4

5# 1. Create Data

6X = np.array([[1], [2], [3], [4]]) # Feature matrix

7y = np.array([2, 4, 6, 8]) # Target vector (y = 2*x)

8

9# 2. Instantiate Model

10model = LinearRegression()

11

12# 3. Fit Model

13model.fit(X, y)

14

15# 4. Predict

16X_new = np.array([[5], [10]])

17predictions = model.predict(X_new)

18

19print(f"Predictions for 5 and 10: {predictions}")

20print(f"Coefficient (Slope): {model.coef_[0]}")

21print(f"Intercept: {model.intercept_}")

22

The Machine Learning Workflow

A real-world project involves more than just fitting a model. It requires rigorous evaluation.

Train/Test Splitting

We must never evaluate a model on the same data it was trained on. Scikit-Learn provides train_test_split to handle this.

Interactive Lab


from sklearn.model_selection import train_test_split
from sklearn.datasets import load_iris

# Load sample data
iris = load_iris()
X, y = iris.data, iris.target

# Split into 80% training, 20% testing
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

print(f"Training set size: {X_train.shape[0]}")
print(f"Testing set size: {X_test.shape[0]}")

Expected output

Training set size: 120
Testing set size: 30

python

1

2from sklearn.model_selection import train_test_split

3from sklearn.datasets import load_iris

4

5# Load sample data

6iris = load_iris()

7X, y = iris.data, iris.target

8

9# Split into 80% training, 20% testing

10X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

11

12print(f"Training set size: {X_train.shape[0]}")

13print(f"Testing set size: {X_test.shape[0]}")

14

Preprocessing: Scaling and Encoding

Machine learning models are sensitive to the scale of features. For example, a model might give more weight to a feature ranging from 0 to 1000 than to one ranging from 0 to 1.

Interactive Lab


from sklearn.preprocessing import StandardScaler
import numpy as np

data = np.array([[10, 0.001], [20, 0.002], [30, 0.003]])
scaler = StandardScaler()

scaled_data = scaler.fit_transform(data)
print("Original Data:\n", data)
print("\nScaled Data (mean=0, std=1):\n", scaled_data)

Expected output

Original Data:
 [[1.e+01 1.e-03]
 [2.e+01 2.e-03]
 [3.e+01 3.e-03]]

Scaled Data (mean=0, std=1):
 [[-1.22474487 -1.22474487]
 [ 0.          0.        ]
 [ 1.22474487  1.22474487]]

python

1

2from sklearn.preprocessing import StandardScaler

3import numpy as np

4

5data = np.array([[10, 0.001], [20, 0.002], [30, 0.003]])

6scaler = StandardScaler()

7

8scaled_data = scaler.fit_transform(data)

9print("Original Data:\n", data)

10print("\nScaled Data (mean=0, std=1):\n", scaled_data)

11

Evaluation Metrics

How do we know if our model is any good? Scikit-Learn offers a suite of metrics.

For Regression: Mean Squared Error (MSE), R-squared.
For Classification: Accuracy, Precision, Recall, F1-Score.

Interactive Lab


from sklearn.metrics import accuracy_score

# Simulated ground truth and predictions
y_true = [0, 1, 2, 0, 1]
y_pred = [0, 2, 1, 0, 1]

acc = accuracy_score(y_true, y_pred)
print(f"Model Accuracy: {acc * 100:.1f}%")

Expected output

Model Accuracy: 60.0%

python

1

2from sklearn.metrics import accuracy_score

3

4# Simulated ground truth and predictions

5y_true = [0, 1, 2, 0, 1]

6y_pred = [0, 2, 1, 0, 1]

7

8acc = accuracy_score(y_true, y_pred)

9print(f"Model Accuracy: {acc * 100:.1f}%")

10

In the next module, we will explore supervised learning models like Decision Trees and Random Forests in much greater detail.

Section Detail

Scikit-Learn: Supervised Learning Deep Dive

Supervised Learning Categories

Supervised learning is divided into two main tasks:

Regression: Predicting a continuous numerical value (e.g., house prices).
Classification: Predicting a discrete category or label (e.g., spam vs. not spam).

Classification: Beyond the Basics

We already saw the Estimator API. Let’s look at more complex classifiers.

Support Vector Machines (SVM)

SVMs are powerful models that attempt to find the hyperplane that best separates classes with the maximum margin.

Interactive Lab


from sklearn.svm import SVC
from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split

# Generate synthetic data
X, y = make_classification(n_samples=100, n_features=2, n_redundant=0, random_state=42)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)

# Use Radial Basis Function (RBF) kernel
model = SVC(kernel='rbf', C=1.0)
model.fit(X_train, y_train)

print(f"SVM Test Accuracy: {model.score(X_test, y_test):.2f}")

Expected output

SVM Test Accuracy: 0.95

python

1

2from sklearn.svm import SVC

3from sklearn.datasets import make_classification

4from sklearn.model_selection import train_test_split

5

6# Generate synthetic data

7X, y = make_classification(n_samples=100, n_features=2, n_redundant=0, random_state=42)

8X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)

9

10# Use Radial Basis Function (RBF) kernel

11model = SVC(kernel='rbf', C=1.0)

12model.fit(X_train, y_train)

13

14print(f"SVM Test Accuracy: {model.score(X_test, y_test):.2f}")

15

Decision Trees and Random Forests

Decision Trees mimic human decision-making by splitting data based on feature thresholds. Random Forests are “ensembles” of many decision trees, which reduces overfitting.

Interactive Lab


from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import load_digits

digits = load_digits()
X, y = digits.data, digits.target

model = RandomForestClassifier(n_estimators=100)
model.fit(X, y)

print(f"Random Forest Accuracy on Digits: {model.score(X, y):.2f}")

Expected output

Random Forest Accuracy on Digits: 1.00

python

1

2from sklearn.ensemble import RandomForestClassifier

3from sklearn.datasets import load_digits

4

5digits = load_digits()

6X, y = digits.data, digits.target

7

8model = RandomForestClassifier(n_estimators=100)

9model.fit(X, y)

10

11print(f"Random Forest Accuracy on Digits: {model.score(X, y):.2f}")

12

Regression: Complexity and Regularization

Simple linear regression often underfits complex data. We can use techniques like Ridge and Lasso regression to prevent overfitting by penalizing large coefficients.

Interactive Lab


from sklearn.linear_model import Ridge, Lasso
import numpy as np

# Synthetic noisy data
X = np.random.rand(100, 10)
y = 2*X[:, 0] + 3*X[:, 1] + np.random.randn(100) * 0.1

ridge = Ridge(alpha=1.0)
ridge.fit(X, y)

lasso = Lasso(alpha=0.1)
lasso.fit(X, y)

print(f"Ridge Coefficients: {ridge.coef_[:2]}")
print(f"Lasso Coefficients: {lasso.coef_[:2]}")

Expected output

Ridge Coefficients: [1.98 2.95]
Lasso Coefficients: [1.85 2.82]

python

1

2from sklearn.linear_model import Ridge, Lasso

3import numpy as np

4

5# Synthetic noisy data

6X = np.random.rand(100, 10)

7y = 2*X[:, 0] + 3*X[:, 1] + np.random.randn(100) * 0.1

8

9ridge = Ridge(alpha=1.0)

10ridge.fit(X, y)

11

12lasso = Lasso(alpha=0.1)

13lasso.fit(X, y)

14

15print(f"Ridge Coefficients: {ridge.coef_[:2]}")

16print(f"Lasso Coefficients: {lasso.coef_[:2]}")

17

Hyperparameter Tuning

How do we choose the best alpha for Ridge or the best n_estimators for a Random Forest? We use Grid Search.

Interactive Lab


from sklearn.model_selection import GridSearchCV
from sklearn.ensemble import RandomForestClassifier

param_grid = {
  'n_estimators': [10, 50, 100],
  'max_depth': [None, 5, 10]
}

grid = GridSearchCV(RandomForestClassifier(), param_grid, cv=5)
# grid.fit(X, y) # Uncomment to run (takes time)
# print(f"Best params: {grid.best_params_}")
print("GridSearchCV defined with cross-validation.")

Expected output

GridSearchCV defined with cross-validation.

python

1

2from sklearn.model_selection import GridSearchCV

3from sklearn.ensemble import RandomForestClassifier

4

5param_grid = {

6 'n_estimators': [10, 50, 100],

7 'max_depth': [None, 5, 10]

8}

9

10grid = GridSearchCV(RandomForestClassifier(), param_grid, cv=5)

11# grid.fit(X, y) # Uncomment to run (takes time)

12# print(f"Best params: {grid.best_params_}")

13print("GridSearchCV defined with cross-validation.")

14

Pipelines: Chaining Transformations

A pipeline combines a series of preprocessing steps and a final estimator into one object. This prevents data leakage during cross-validation.

Interactive Lab


from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC

pipeline = Pipeline([
  ('scaler', StandardScaler()),
  ('svc', SVC())
])

# Fit and predict just like a single estimator
# pipeline.fit(X_train, y_train)
print("Pipeline created: [Scaler -> SVC]")

Expected output

Pipeline created: [Scaler -> SVC]

python

1

2from sklearn.pipeline import Pipeline

3from sklearn.preprocessing import StandardScaler

4from sklearn.svm import SVC

5

6pipeline = Pipeline([

7 ('scaler', StandardScaler()),

8 ('svc', SVC())

9])

10

11# Fit and predict just like a single estimator

12# pipeline.fit(X_train, y_train)

13print("Pipeline created: [Scaler -> SVC]")

14

In the next module, we’ll explore Unsupervised Learning techniques like Clustering and Dimensionality Reduction.

Section Detail

Scikit-Learn: Unsupervised Learning

Learning Without Labels

Unsupervised learning is used when we have features ( $X$ ) but no target ( $y$ ). The goal is to find inherent patterns or structures within the data.

1. Clustering: Grouping Similar Samples

Clustering algorithms attempt to partition the data into groups (clusters) where samples in the same group are more similar to each other than to those in other groups.

K-Means Clustering

The most popular clustering algorithm. it partitions data into $K$ clusters by minimizing the distance between points and their cluster centroids.

Interactive Lab


from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt

# Create 3 distinct blobs of data
X, _ = make_blobs(n_samples=300, centers=3, cluster_std=0.60, random_state=0)

# Instantiate and fit
kmeans = KMeans(n_clusters=3)
kmeans.fit(X)

# Predicted labels
y_kmeans = kmeans.predict(X)
print(f"Centroids:\n{kmeans.cluster_centers_}")

Expected output

Centroids:
[[ 0.94973332  4.41906906]
 [ 2.01527128  0.92212359]
 [-1.58438467  2.85831692]]

python

1

2from sklearn.cluster import KMeans

3from sklearn.datasets import make_blobs

4import matplotlib.pyplot as plt

5

6# Create 3 distinct blobs of data

7X, _ = make_blobs(n_samples=300, centers=3, cluster_std=0.60, random_state=0)

8

9# Instantiate and fit

10kmeans = KMeans(n_clusters=3)

11kmeans.fit(X)

12

13# Predicted labels

14y_kmeans = kmeans.predict(X)

15print(f"Centroids:\n{kmeans.cluster_centers_}")

16

2. Dimensionality Reduction

High-dimensional data (hundreds or thousands of features) is difficult to visualize and can lead to the “curse of dimensionality.” Dimensionality reduction seeks to represent data in a lower-dimensional space while preserving as much information as possible.

Principal Component Analysis (PCA)

PCA finds the “principal components”—the orthogonal axes along which the data varies the most.

Interactive Lab


from sklearn.decomposition import PCA
from sklearn.datasets import load_iris

iris = load_iris()
X = iris.data # 4 features

# Reduce from 4D to 2D
pca = PCA(n_components=2)
X_reduced = pca.fit_transform(X)

print(f"Original shape: {X.shape}")
print(f"Reduced shape: {X_reduced.shape}")
print(f"Explained variance ratio: {pca.explained_variance_ratio_}")

Expected output

Original shape: (150, 4)
Reduced shape: (150, 2)
Explained variance ratio: [0.92461872 0.05306648]

python

1

2from sklearn.decomposition import PCA

3from sklearn.datasets import load_iris

4

5iris = load_iris()

6X = iris.data # 4 features

7

8# Reduce from 4D to 2D

9pca = PCA(n_components=2)

10X_reduced = pca.fit_transform(X)

11

12print(f"Original shape: {X.shape}")

13print(f"Reduced shape: {X_reduced.shape}")

14print(f"Explained variance ratio: {pca.explained_variance_ratio_}")

15

3. Anomaly Detection

Identifying outliers or unusual patterns that do not conform to expected behavior.

Isolation Forest

Isolation Forest works by isolating anomalies using trees. Anomalies are easier to isolate and thus have shorter path lengths in the trees.

Interactive Lab


from sklearn.ensemble import IsolationForest
import numpy as np

# Generate normal data
X_train = 0.3 * np.random.randn(100, 2)
X_train = np.r_[X_train + 2, X_train - 2]

# Generate some abnormal novel observations
X_outliers = np.random.uniform(low=-4, high=4, size=(20, 2))

clf = IsolationForest(contamination=0.1)
clf.fit(X_train)

y_pred_outliers = clf.predict(X_outliers)
print(f"Predictions for outliers (-1 is anomaly):\n{y_pred_outliers}")

Expected output

Predictions for outliers (-1 is anomaly):
[-1 -1 -1  1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1  1]

python

1

2from sklearn.ensemble import IsolationForest

3import numpy as np

4

5# Generate normal data

6X_train = 0.3 * np.random.randn(100, 2)

7X_train = np.r_[X_train + 2, X_train - 2]

8

9# Generate some abnormal novel observations

10X_outliers = np.random.uniform(low=-4, high=4, size=(20, 2))

11

12clf = IsolationForest(contamination=0.1)

13clf.fit(X_train)

14

15y_pred_outliers = clf.predict(X_outliers)

16print(f"Predictions for outliers (-1 is anomaly):\n{y_pred_outliers}")

17

Summary of Unsupervised Techniques

Code

skinparam componentStyle rectangle

package "Unsupervised Tasks" {
component "Clustering" as CLUST
component "Dimensionality Reduction" as DIM
component "Density Estimation" as DENS
}

node "Algorithms" {
component "K-Means / DBSCAN / Agglomerative" as CALG
component "PCA / t-SNE / UMAP" as DALG
component "GMM / Kernel Density" as DEALG
}

CLUST --> CALG
DIM --> DALG
DENS --> DEALG

In the next section, we will explore specialized scientific libraries like SciPy and Matplotlib to round out our scientific Python toolkit.

Advanced Scientific Stack

Section Detail

SciPy: Essential Scientific Algorithms

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.]

python

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

python

1

2from scipy.integrate import quad

3import numpy as np

4

5# Integrate exp(-x^2) from 0 to infinity (The Gaussian integral)

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

python

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# Create linear and cubic interpolation functions

9f_linear = interp1d(x, y)

10f_cubic = interp1d(x, y, kind='cubic')

11

12xnew = np.linspace(0, 10, num=41, endpoint=True)

13# Evaluate at new points

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.

python

1

2from scipy import signal

3import numpy as np

4

5# Create a noisy signal

6t = np.linspace(0, 1, 100)

7sig = np.sin(2*np.pi*10*t) + np.random.randn(100)*0.1

8

9# Apply a Butterworth filter

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

python

1

2from scipy import stats

3import numpy as np

4

5# Normal distribution

6data = np.random.normal(0, 1, 1000)

7

8# Calculate parameters

9mean, std = stats.norm.fit(data)

10print(f"Fitted Mean: {mean:.2f}, Std: {std:.2f}")

11

12# Perform a T-test

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.

Data Visualization Deep Dive

Section Detail

Matplotlib: The Figure and Axes Architecture

Beyond “Simple” Plotting

Many beginners start with plt.plot() and quickly get frustrated when they can’t control the details of their figure. To master Matplotlib, you must understand its Object-Oriented (OO) Architecture.

The Hierarchy of a Plot

A Matplotlib plot is not just a bunch of lines; it’s a hierarchy of objects called Artists.

Code

skinparam componentStyle rectangle

package "Figure" as Fig {
component "Axes" as AX1
}

component "Axis" as Axis
component "Title" as Title
component "Legend" as Legend
component "Ticks" as Ticks
component "Grid" as Grid
component "Lines" as Lines

AX1 *-- Axis
AX1 *-- Title
AX1 *-- Legend
AX1 *-- Ticks
AX1 *-- Grid
AX1 *-- Lines

note bottom of AX1: "A Figure can contain\nmultiple Axes (subplots)."

The Figure vs. The Axes

Figure: The top-level container. Think of it as the physical page or empty window.
Axes: A “subplot.” It is the area where the data is actually plotted. Most operations (setting titles, labels, grids) happen on an Axes object.

Creating Subplots the Pro Way

The most common way to start a plot is with plt.subplots().

OO Plotting Example


import matplotlib.pyplot as plt
import numpy as np

# Create figure and axes
fig, ax = plt.subplots(figsize=(10, 5))

# Generate data
x = np.linspace(0, 10, 100)
y = np.sin(x)

# Plotting on the Axes object
ax.plot(x, y, label='Sine Wave', color='#2563eb', linewidth=2.5)

# Setting attributes on the Axes object
ax.set_title('Hierarchical Plot Construction', fontsize=16, fontweight='bold')
ax.set_xlabel('Time (s)')
ax.set_ylabel('Amplitude')
ax.grid(True, linestyle='--', alpha=0.6)
ax.legend()

plt.show()

python

1

2import matplotlib.pyplot as plt

3import numpy as np

4

5# Create figure and axes

6fig, ax = plt.subplots(figsize=(10, 5))

7

8# Generate data

9x = np.linspace(0, 10, 100)

10y = np.sin(x)

11

12# Plotting on the Axes object

13ax.plot(x, y, label='Sine Wave', color='#2563eb', linewidth=2.5)

14

15# Setting attributes on the Axes object

16ax.set_title('Hierarchical Plot Construction', fontsize=16, fontweight='bold')

17ax.set_xlabel('Time (s)')

18ax.set_ylabel('Amplitude')

19ax.grid(True, linestyle='--', alpha=0.6)

20ax.legend()

21

22plt.show()

23

Anatomy of a Figure

Each element you see—the title, the tick labels, the lines—is a separate Artist object that can be customized.

Artist	Purpose
Line2D	Represents the data lines.
Text	Used for titles, labels, and annotations.
Axis	Manages the scale, limits, and tick locations.
Patch	Used for boxes, circles, and bars.

Multiple Subplots

The power of the OO API becomes clear when you have multiple plots in a single figure.

Multiple Subplots


import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(0, 10, 100)

# Create a 2x1 grid of plots
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 8), sharex=True)

ax1.plot(x, np.sin(x), 'r')
ax1.set_title('Sine')

ax2.plot(x, np.cos(x), 'b')
ax2.set_title('Cosine')

# Automatically adjust spacing
plt.tight_layout()
plt.show()

python

1

2import matplotlib.pyplot as plt

3import numpy as np

4

5x = np.linspace(0, 10, 100)

6

7# Create a 2x1 grid of plots

8fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 8), sharex=True)

9

10ax1.plot(x, np.sin(x), 'r')

11ax1.set_title('Sine')

12

13ax2.plot(x, np.cos(x), 'b')

14ax2.set_title('Cosine')

15

16# Automatically adjust spacing

17plt.tight_layout()

18plt.show()

19

In the next lesson, we will explore how to style these components to create publication-quality figures.

Section Detail

Matplotlib: Styling and Customization

Plotting for Publication

Scientists don’t just produce plots; they produce evidence. A good plot should be clear, accessible, and aesthetically pleasing. In this lesson, we explore how to move beyond default settings.

Global Styles with `plt.style`

Matplotlib comes with several pre-defined style sheets. You can view them all using plt.style.available.

Switching Styles


import matplotlib.pyplot as plt
import numpy as np

# Use a clean, modern style
plt.style.use('ggplot')

x = np.random.randn(100)
y = np.random.randn(100)

fig, ax = plt.subplots()
ax.scatter(x, y, alpha=0.5, color='purple')
ax.set_title('Styled Scatter Plot')

plt.show()

python

1

2import matplotlib.pyplot as plt

3import numpy as np

4

5# Use a clean, modern style

6plt.style.use('ggplot')

7

8x = np.random.randn(100)

9y = np.random.randn(100)

10

11fig, ax = plt.subplots()

12ax.scatter(x, y, alpha=0.5, color='purple')

13ax.set_title('Styled Scatter Plot')

14

15plt.show()

16

Working with Colormaps

Choosing the right colormap is critical. Avoid “jet” and other non-perceptually uniform maps. Use viridis, magma, or inferno for continuous data.

Colormap Comparison


import matplotlib.pyplot as plt
import numpy as np

x = np.linspace(0, 5, 100)
y = np.linspace(0, 5, 100)
X, Y = np.meshgrid(x, y)
Z = np.sin(X)**10 + np.cos(10 + Y*X) * np.cos(X)

fig, ax = plt.subplots()
im = ax.imshow(Z, origin='lower', extent=[0, 5, 0, 5], cmap='viridis')
fig.colorbar(im, label='Intensity')
ax.set_title('Perceptually Uniform Colormap')

plt.show()

python

1

2import matplotlib.pyplot as plt

3import numpy as np

4

5x = np.linspace(0, 5, 100)

6y = np.linspace(0, 5, 100)

7X, Y = np.meshgrid(x, y)

8Z = np.sin(X)**10 + np.cos(10 + Y*X) * np.cos(X)

9

10fig, ax = plt.subplots()

11im = ax.imshow(Z, origin='lower', extent=[0, 5, 0, 5], cmap='viridis')

12fig.colorbar(im, label='Intensity')

13ax.set_title('Perceptually Uniform Colormap')

14

15plt.show()

16

Annotations: Highlighting Key Data

Sometimes you need to point directly at a feature in your data.

Using Annotations


import matplotlib.pyplot as plt
import numpy as np

fig, ax = plt.subplots()

t = np.arange(0.0, 5.0, 0.01)
s = np.cos(2*np.pi*t)
ax.plot(t, s, lw=2)

# Point to the peak
ax.annotate('Local Maximum', xy=(2, 1), xytext=(3, 1.5),
          arrowprops=dict(facecolor='black', shrink=0.05))

ax.set_ylim(-2, 2)
plt.show()

python

1

2import matplotlib.pyplot as plt

3import numpy as np

4

5fig, ax = plt.subplots()

6

7t = np.arange(0.0, 5.0, 0.01)

8s = np.cos(2*np.pi*t)

9ax.plot(t, s, lw=2)

10

11# Point to the peak

12ax.annotate('Local Maximum', xy=(2, 1), xytext=(3, 1.5),

13 arrowprops=dict(facecolor='black', shrink=0.05))

14

15ax.set_ylim(-2, 2)

16plt.show()

17

Complex Layouts with GridSpec

While plt.subplots() handles simple grids, GridSpec allows for plots that span multiple rows or columns.

GridSpec for Complex Dashboards


import matplotlib.pyplot as plt

fig = plt.figure(figsize=(8, 6))
gs = fig.add_gridspec(3, 3)

ax1 = fig.add_subplot(gs[0, :])
ax1.set_title('Top row, full width')

ax2 = fig.add_subplot(gs[1:, :2])
ax2.set_title('Bottom left, big')

ax3 = fig.add_subplot(gs[1, 2])
ax3.set_title('Small 1')

ax4 = fig.add_subplot(gs[2, 2])
ax4.set_title('Small 2')

plt.tight_layout()
plt.show()

python

1

2import matplotlib.pyplot as plt

3

4fig = plt.figure(figsize=(8, 6))

5gs = fig.add_gridspec(3, 3)

6

7ax1 = fig.add_subplot(gs[0, :])

8ax1.set_title('Top row, full width')

9

10ax2 = fig.add_subplot(gs[1:, :2])

11ax2.set_title('Bottom left, big')

12

13ax3 = fig.add_subplot(gs[1, 2])

14ax3.set_title('Small 1')

15

16ax4 = fig.add_subplot(gs[2, 2])

17ax4.set_title('Small 2')

18

19plt.tight_layout()

20plt.show()

21

In the next lesson, we will see how Seaborn makes these complex statistical visualizations much easier to achieve.

Section Detail

Seaborn: Statistical Data Visualization

High-Level Statistical Plotting

If Matplotlib is the “foundational” layer, Seaborn is the “sophisticated” layer. Built on top of Matplotlib, it integrates deeply with Pandas and automates many complex statistical visualization tasks.

Why Seaborn?

DataFrame Integration: You can pass column names directly as strings.
Beautiful Defaults: Better color palettes and styles out of the box.
Statistical Aggregation: It automatically calculates confidence intervals and means for you.

Exploring Relationships with `relplot`

relplot is the primary entry point for visualizing the relationship between two variables.

Multivariate Regression with Seaborn


import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

# Create synthetic tips-like data
df = pd.DataFrame({
  'total_bill': np.random.uniform(10, 50, 100),
  'tip': np.random.uniform(1, 10, 100),
  'smoker': np.random.choice(['Yes', 'No'], 100),
  'day': np.random.choice(['Thur', 'Fri', 'Sat', 'Sun'], 100)
})

# Add a correlation
df['tip'] = df['total_bill'] * 0.15 + np.random.normal(0, 1, 100)

sns.set_theme(style="ticks")

# Plot relationship by smoker status and day
g = sns.relplot(
  data=df, x="total_bill", y="tip",
  col="day", hue="smoker", style="smoker",
  kind="scatter"
)

plt.show()

python

1

2import seaborn as sns

3import matplotlib.pyplot as plt

4import pandas as pd

5import numpy as np

6

7# Create synthetic tips-like data

8df = pd.DataFrame({

9 'total_bill': np.random.uniform(10, 50, 100),

10 'tip': np.random.uniform(1, 10, 100),

11 'smoker': np.random.choice(['Yes', 'No'], 100),

12 'day': np.random.choice(['Thur', 'Fri', 'Sat', 'Sun'], 100)

13})

14

15# Add a correlation

16df['tip'] = df['total_bill'] * 0.15 + np.random.normal(0, 1, 100)

17

18sns.set_theme(style="ticks")

19

20# Plot relationship by smoker status and day

21g = sns.relplot(

22 data=df, x="total_bill", y="tip",

23 col="day", hue="smoker", style="smoker",

24 kind="scatter"

25)

26

27plt.show()

28

Categorical Data and Distributions

When dealing with categories, catplot and displot are your best friends.

Boxen Plots and Distributions


import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

data = pd.DataFrame({
  'Value': np.concatenate([np.random.normal(0, 1, 100), np.random.normal(2, 0.5, 100)]),
  'Group': ['A']*100 + ['B']*100
})

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# Violin plot
sns.violinplot(data=data, x="Group", y="Value", ax=ax1, palette="muted")
ax1.set_title("Distribution by Group (Violin)")

# ECDF (Empirical Cumulative Distribution Function)
sns.ecdfplot(data=data, x="Value", hue="Group", ax=ax2)
ax2.set_title("CDF Comparison")

plt.show()

python

1

2import seaborn as sns

3import matplotlib.pyplot as plt

4import pandas as pd

5import numpy as np

6

7data = pd.DataFrame({

8 'Value': np.concatenate([np.random.normal(0, 1, 100), np.random.normal(2, 0.5, 100)]),

9 'Group': ['A']*100 + ['B']*100

10})

11

12fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

13

14# Violin plot

15sns.violinplot(data=data, x="Group", y="Value", ax=ax1, palette="muted")

16ax1.set_title("Distribution by Group (Violin)")

17

18# ECDF (Empirical Cumulative Distribution Function)

19sns.ecdfplot(data=data, x="Value", hue="Group", ax=ax2)

20ax2.set_title("CDF Comparison")

21

22plt.show()

23

Regression Models with `lmplot`

Seaborn can perform linear regression right inside the visualization.

Linear Regression Visualization


import seaborn as sns
import matplotlib.pyplot as plt

df = sns.load_dataset("anscombe").query("dataset == 'I'")

sns.lmplot(x="x", y="y", data=df, ci=95, scatter_kws={"s": 80})
plt.title("Statistical Regression Analysis")
plt.show()

python

1

2import seaborn as sns

3import matplotlib.pyplot as plt

4

5df = sns.load_dataset("anscombe").query("dataset == 'I'")

6

7sns.lmplot(x="x", y="y", data=df, ci=95, scatter_kws={"s": 80})

8plt.title("Statistical Regression Analysis")

9plt.show()

10

Heatmaps and Matrix Plots

Heatmaps are essential for visualizing correlation matrices.

Correlation Heatmap


import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

# Correlation matrix for synthetic features
data = np.random.rand(10, 10)
cols = [f'Feat_{i}' for i in range(10)]
df = pd.DataFrame(data, columns=cols)
corr = df.corr()

plt.figure(figsize=(10, 8))
sns.heatmap(corr, annot=True, cmap='coolwarm', fmt=".2f")
plt.title("Feature Correlation Matrix")
plt.show()

python

1

2import seaborn as sns

3import matplotlib.pyplot as plt

4import numpy as np

5import pandas as pd

6

7# Correlation matrix for synthetic features

8data = np.random.rand(10, 10)

9cols = [f'Feat_{i}' for i in range(10)]

10df = pd.DataFrame(data, columns=cols)

11corr = df.corr()

12

13plt.figure(figsize=(10, 8))

14sns.heatmap(corr, annot=True, cmap='coolwarm', fmt=".2f")

15plt.title("Feature Correlation Matrix")

16plt.show()

17

In the final lesson of this module, we will discuss the principles of scientific storytelling and how to choose the right chart for your data.

Section Detail

Scientific Storytelling: The Art of the Visual

Visualization as Communication

Data visualization is the language of science. A poorly designed plot can obscure the truth, while a great one can reveal deep insights that numbers alone cannot convey.

The Principles of Tufte

Edward Tufte, a pioneer in data visualization, emphasized the concept of Data-Ink Ratio. The goal is to maximize the ink used for the data and minimize the ink used for everything else (decorative borders, unnecessary grids, etc.).

Tufte’s Rules:

Above all else, show the data.
Maximize the data-ink ratio.
Erase non-data-ink.
Erase redundant data-ink.
Revise and edit.

Avoiding Common Pitfalls

1. The Truncated Y-Axis

Starting a Y-axis at a non-zero value can exaggerate small differences. While sometimes necessary, it should be done with caution and clear labeling.

2. Overplotting

Too many points in a scatter plot can hide the density. Use alpha transparency or hexbins to solve this.

Handling Overplotting


import matplotlib.pyplot as plt
import numpy as np

# Generate thousands of points
n = 10000
x = np.random.standard_normal(n)
y = np.random.standard_normal(n)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# Bad: Solid points hide density
ax1.scatter(x, y, s=5, alpha=1)
ax1.set_title("Opaque (Hidden Density)")

# Good: Hexbin shows density
hb = ax2.hexbin(x, y, gridsize=40, cmap='inferno')
fig.colorbar(hb, ax=ax2, label='Counts')
ax2.set_title("Hexbin (Clear Density)")

plt.show()

python

1

2import matplotlib.pyplot as plt

3import numpy as np

4

5# Generate thousands of points

6n = 10000

7x = np.random.standard_normal(n)

8y = np.random.standard_normal(n)

9

10fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

11

12# Bad: Solid points hide density

13ax1.scatter(x, y, s=5, alpha=1)

14ax1.set_title("Opaque (Hidden Density)")

15

16# Good: Hexbin shows density

17hb = ax2.hexbin(x, y, gridsize=40, cmap='inferno')

18fig.colorbar(hb, ax=ax2, label='Counts')

19ax2.set_title("Hexbin (Clear Density)")

20

21plt.show()

22

Choosing the Right Plot

Data Type	Best Plot	Why?
Trend over Time	Line Plot	Emphasizes continuity and sequence.
Comparing Categories	Bar Plot / Box Plot	Clear separation and ranking of groups.
Relationships (2 vars)	Scatter Plot	Shows correlation or lack thereof.
Distribution	Histogram / KDE	Shows the shape and spread of data.
Composition	Stacked Area	Shows how parts change relative to the whole.

Narrative Visualization

A good figure should walk the reader through a story. Use labels and arrows to highlight the “turning points” in your data.

Narrative Plotting


import matplotlib.pyplot as plt
import numpy as np

time = np.linspace(0, 10, 100)
signal = np.exp(-0.5 * time) * np.sin(2 * np.pi * time)

fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(time, signal, color='darkblue', label='Sensor A Output')

# Add narrative elements
ax.axhline(0, color='black', alpha=0.3)
ax.annotate('System Damping Begins', xy=(2.3, 0.2), xytext=(5, 0.5),
          arrowprops=dict(facecolor='black', shrink=0.05),
          fontsize=12, fontweight='bold')

ax.fill_between(time, signal, 0, alpha=0.1, color='blue')

ax.set_title("Signal Decay Analysis", loc='left', fontsize=18)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

plt.show()

python

1

2import matplotlib.pyplot as plt

3import numpy as np

4

5time = np.linspace(0, 10, 100)

6signal = np.exp(-0.5 * time) * np.sin(2 * np.pi * time)

7

8fig, ax = plt.subplots(figsize=(10, 5))

9ax.plot(time, signal, color='darkblue', label='Sensor A Output')

10

11# Add narrative elements

12ax.axhline(0, color='black', alpha=0.3)

13ax.annotate('System Damping Begins', xy=(2.3, 0.2), xytext=(5, 0.5),

14 arrowprops=dict(facecolor='black', shrink=0.05),

15 fontsize=12, fontweight='bold')

16

17ax.fill_between(time, signal, 0, alpha=0.1, color='blue')

18

19ax.set_title("Signal Decay Analysis", loc='left', fontsize=18)

20ax.spines['top'].set_visible(False)

21ax.spines['right'].set_visible(False)

22

23plt.show()

24

Conclusion of the module

You have now moved from plotting basic lines to constructing complex, data-driven narratives. By combining the power of Matplotlib and Seaborn with the principles of Tufte, you can communicate your scientific findings with clarity and impact.

Advanced Scientific Stack

Section Detail

Capstone: Integrating the Full Scientific Stack

The Power of Integration

Throughout this course, we have looked at libraries in isolation. However, the true strength of the Python scientific ecosystem is how seamlessly these tools work together.

In this final project, we will:

Generate synthetic data using NumPy.
Structure and analyze it using Pandas.
Build a predictive model using Scikit-Learn.
Visualize the results using Matplotlib and Seaborn.

The Scenario: Predicting Sensor Failures

Imagine we have a chemical reactor with two sensors: Temperature and Pressure. We want to predict if the reactor is in a “Stable” or “Unstable” state based on these readings.

End-to-End Scientific Workflow


import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import confusion_matrix

# 1. Generate Synthetic Data (NumPy)
np.random.seed(42)
n_samples = 200
temp = np.random.normal(300, 20, n_samples)
pressure = np.random.normal(50, 10, n_samples)

# Define stability condition: unstable if temp > 330 or pressure > 65
stability = ((temp > 330) | (pressure > 65)).astype(int)

# 2. Structure Data (Pandas)
df = pd.DataFrame({
  'Temperature': temp,
  'Pressure': pressure,
  'Unstable': stability
})

print("First 5 samples of the dataset:")
print(df.head())

# 3. Model Training (Scikit-Learn)
X = df[['Temperature', 'Pressure']]
y = df['Unstable']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25)

model = RandomForestClassifier(n_estimators=50)
model.fit(X_train, y_train)

# 4. Comprehensive Visualization (Matplotlib/Seaborn)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

# Plot A: The raw data distribution
sns.scatterplot(data=df, x='Temperature', y='Pressure', hue='Unstable', ax=ax1, palette='coolwarm')
ax1.set_title("Reactor States: Temperature vs Pressure")

# Plot B: The Model Performance (Confusion Matrix)
y_pred = model.predict(X_test)
cm = confusion_matrix(y_test, y_pred)
sns.heatmap(cm, annot=True, fmt='d', cmap='Blues', ax=ax2)
ax2.set_xlabel('Predicted')
ax2.set_ylabel('Actual')
ax2.set_title("Model Accuracy: Confusion Matrix")

plt.tight_layout()
plt.show()

score = model.score(X_test, y_test)
print(f"\nModel Accuracy on Test Set: {score*100:.1f}%")

python

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import confusion_matrix

# 1. Generate Synthetic Data (NumPy)
np.random.seed(42)
n_samples = 200
temp = np.random.normal(300, 20, n_samples)
pressure = np.random.normal(50, 10, n_samples)

# Define stability condition: unstable if temp > 330 or pressure > 65
stability = ((temp > 330) | (pressure > 65)).astype(int)

# 2. Structure Data (Pandas)
df = pd.DataFrame({
  'Temperature': temp,
  'Pressure': pressure,
  'Unstable': stability
})

print("First 5 samples of the dataset:")
print(df.head())

# 3. Model Training (Scikit-Learn)
X = df[['Temperature', 'Pressure']]
y = df['Unstable']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25)

model = RandomForestClassifier(n_estimators=50)
model.fit(X_train, y_train)

# 4. Comprehensive Visualization (Matplotlib/Seaborn)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

# Plot A: The raw data distribution
sns.scatterplot(data=df, x='Temperature', y='Pressure', hue='Unstable', ax=ax1, palette='coolwarm')
ax1.set_title("Reactor States: Temperature vs Pressure")

# Plot B: The Model Performance (Confusion Matrix)
y_pred = model.predict(X_test)
cm = confusion_matrix(y_test, y_pred)
sns.heatmap(cm, annot=True, fmt='d', cmap='Blues', ax=ax2)
ax2.set_xlabel('Predicted')
ax2.set_ylabel('Actual')
ax2.set_title("Model Accuracy: Confusion Matrix")

plt.tight_layout()
plt.show()

score = model.score(X_test, y_test)
print(f"\nModel Accuracy on Test Set: {score*100:.1f}%")

1

2import numpy as np

3import pandas as pd

4import matplotlib.pyplot as plt

5import seaborn as sns

6from sklearn.model_selection import train_test_split

7from sklearn.ensemble import RandomForestClassifier

8from sklearn.metrics import confusion_matrix

9

10# 1. Generate Synthetic Data (NumPy)

11np.random.seed(42)

12n_samples = 200

13temp = np.random.normal(300, 20, n_samples)

14pressure = np.random.normal(50, 10, n_samples)

15

16# Define stability condition: unstable if temp > 330 or pressure > 65

17stability = ((temp > 330) | (pressure > 65)).astype(int)

18

19# 2. Structure Data (Pandas)

20df = pd.DataFrame({

21 'Temperature': temp,

22 'Pressure': pressure,

23 'Unstable': stability

24})

25

26print("First 5 samples of the dataset:")

27print(df.head())

28

29# 3. Model Training (Scikit-Learn)

30X = df[['Temperature', 'Pressure']]

31y = df['Unstable']

32

33X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25)

34

35model = RandomForestClassifier(n_estimators=50)

36model.fit(X_train, y_train)

37

38# 4. Comprehensive Visualization (Matplotlib/Seaborn)

39fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

40

41# Plot A: The raw data distribution

42sns.scatterplot(data=df, x='Temperature', y='Pressure', hue='Unstable', ax=ax1, palette='coolwarm')

43ax1.set_title("Reactor States: Temperature vs Pressure")

44

45# Plot B: The Model Performance (Confusion Matrix)

46y_pred = model.predict(X_test)

47cm = confusion_matrix(y_test, y_pred)

48sns.heatmap(cm, annot=True, fmt='d', cmap='Blues', ax=ax2)

49ax2.set_xlabel('Predicted')

50ax2.set_ylabel('Actual')

51ax2.set_title("Model Accuracy: Confusion Matrix")

52

53plt.tight_layout()

54plt.show()

55

56score = model.score(X_test, y_test)

57print(f"\nModel Accuracy on Test Set: {score*100:.1f}%")

58

Course Conclusion

Congratulations! You have mastered the core components of the Python scientific stack. You are now equipped to:

Perform high-performance numerical computing with NumPy.
Derive exact mathematical formulas with SymPy.
Clean and analyze massive datasets with Pandas.
Solve complex scientific problems with SciPy.
Build machine learning pipelines with Scikit-Learn.
Communicate your insights through professional-grade visualizations with Matplotlib and Seaborn.

The journey doesn’t end here. The ecosystem is constantly evolving with libraries like PyTorch for Deep Learning, Dask for Parallel Computing, and Plotly for Interactive Dashboards.

Happy coding!

Python for Scientific Computing

Contents

Python Basics

NumPy Deep Dive

SymPy Deep Dive

Pandas Deep Dive

Scikit-Learn Deep Dive

Advanced Scientific Stack

Data Visualization Deep Dive

Advanced Scientific Stack

Python Basics

The Philosophy of Python

The Zen of Python

Why Python for Science?

1. The “Glue Language” Property

2. Rich Library Ecosystem

3. Community and Documentation

The Execution Model

Getting Started

NumPy Deep Dive

The Importance of NumPy

Why not use Python Lists?

The Problem with Python Lists:

The NumPy Solution: ndarray

Creating Arrays

Array Attributes

Vectorization: The Secret Sauce

Memory Layout and Slicing

Multi-Dimensional Indexing

Comma-Separated Indexing

Integer Array Indexing (Fancy Indexing)

Boolean Indexing (Masking)

Shape Manipulation

Reshape

Stacking and Splitting

Broadcasting: A Preview

Efficient Operations with Broadcasting

The Mental Model

The Four Rules of Broadcasting

Example 1: Scalar and Array

Example 2: Vector and Matrix

Example 3: Incompatible Shapes

Broadcasting in Action

Centering an Array

Performance Benefits

SymPy Deep Dive

Introduction to Symbolic Computation

Symbols and Expressions

Structural Manipulations

Expansion and Factoring

Simplification

Calculus with SymPy

Differentiation

Integration

Limits

Equation Solving

Pretty Printing and LaTeX

Bridging SymPy and NumPy: lambdify

Beyond Simple Algebra

Taylor Series Expansions

Solving Ordinary Differential Equations (ODEs)

Example: The Harmonic Oscillator

Multivariable Calculus

Symbolic Linear Algebra

Global Optimization (Symbolic)

Symbolic Linear Algebra

Creating Symbolic Matrices

Determinants and Inversion

Eigenvalues and Eigenvectors

Matrix Decompositions

Solving Linear Systems (Ax=bAx = b)

Applications in Engineering

Pandas Deep Dive

Why Pandas?

Core Data Structures: Series and DataFrame

The Series Object

The DataFrame Object

Data Selection and Indexing

1. loc: Label-based selection

2. iloc: Integer-based selection

The NumPy Solution: `ndarray`

Bridging SymPy and NumPy: `lambdify`

Solving Linear Systems ( $Ax = b$ )

1. `loc`: Label-based selection

2. `iloc`: Integer-based selection

Transforming Data: `map` and `apply`

The `DatetimeIndex`

General Integration (`quad`)