{ "cells": [ { "cell_type": "markdown", "id": "44440a4a", "metadata": {}, "source": [ "# Module 5 — Linear Algebra & Random\n" ] }, { "cell_type": "markdown", "id": "3ae79c2f", "metadata": {}, "source": [ "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/inigoochoa/comprehensive-numpy/blob/main/module_5_linalg_random.ipynb)\n" ] }, { "cell_type": "markdown", "id": "9201987a", "metadata": {}, "source": [ "## Welcome\n", "\n", "This module covers two pillars of scientific computing with NumPy: linear\n", "algebra (vectors, matrices, and the equations that relate them) and random\n", "number generation (sampling, simulation, and reproducibility). Together\n", "they're the backbone of everything from graphics to machine learning to\n", "Monte Carlo simulation — which is exactly where this module ends up.\n", "\n", "**Learning goals**\n", "- Compute dot products, norms, and cosine similarity between vectors\n", "- Multiply matrices with `@` and know how it differs from elementwise `*`\n", "- Solve linear systems with `np.linalg.solve` instead of inverting matrices by hand\n", "- Find eigenvalues/eigenvectors and verify them directly from their definition\n", "- Generate reproducible random numbers with the modern `Generator` API\n", "- Combine random sampling with vectorized boolean indexing to run a simulation\n", "\n", "**Setup** — run this first.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "a1c8ecb9", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "print(np.__version__)\n" ] }, { "cell_type": "markdown", "id": "30588e8a", "metadata": {}, "source": [ "## Lesson 1 — Vectors & Dot Products\n", "\n", "The **dot product** of two same-length vectors, `np.dot(a, b)` or simply\n", "`a @ b`, multiplies corresponding elements and sums the results. It shows\n", "up everywhere: it's the core of matrix multiplication (next lesson), and\n", "`np.linalg.norm(v)` — a vector's length/magnitude — is just\n", "`sqrt(v @ v)`.\n", "\n", "A useful derived quantity is **cosine similarity**, which measures the\n", "angle between two vectors regardless of their length:\n", "\n", "```\n", "cos(theta) = (a @ b) / (||a|| * ||b||)\n", "```\n", "\n", "It ranges from `-1` (opposite direction) to `1` (same direction), with `0`\n", "meaning the vectors are orthogonal.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d7f4dfdd", "metadata": {}, "outputs": [], "source": [ "a = np.array([1.0, 2.0, 3.0])\n", "b = np.array([4.0, 5.0, 6.0])\n", "\n", "print(\"dot: \", np.dot(a, b), a @ b)\n", "print(\"norm a:\", np.linalg.norm(a), np.sqrt(a @ a))\n", "\n", "cos_theta = (a @ b) / (np.linalg.norm(a) * np.linalg.norm(b))\n", "print(\"cos similarity:\", cos_theta)\n" ] }, { "cell_type": "markdown", "id": "d5659773", "metadata": {}, "source": [ "### 💪 Exercise 1.1: Unit vector\n", "\n", "Write `to_unit_vector(v)` that returns `v` scaled to have length (norm)\n", "exactly `1`, pointing in the same direction. Use `np.linalg.norm`.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d01df601", "metadata": {}, "outputs": [], "source": [ "def to_unit_vector(v):\n", " # YOUR CODE HERE\n", " raise NotImplementedError\n" ] }, { "cell_type": "code", "execution_count": null, "id": "e2f550bb", "metadata": {}, "outputs": [], "source": [ "result = to_unit_vector(np.array([3.0, 4.0]))\n", "print(\"to_unit_vector(np.array([3.0, 4.0])) ->\")\n", "print(result)\n", "\n", "expected = np.array([0.6, 0.8])\n", "_passed = np.allclose(np.linalg.norm(result), 1.0) and np.allclose(result, expected)\n", "\n", "if not _passed:\n", " print()\n", " print(\"Expected:\")\n", " print(expected)\n" ] }, { "cell_type": "markdown", "id": "6d781345", "metadata": {}, "source": [ "
\n", "💡 Hint (click to expand)\n", "\n", "Divide `v` by `np.linalg.norm(v)`.\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "bd951563", "metadata": {}, "source": [ "
\n", "✅ Solution (click to expand)\n", "\n", "```python\n", "def to_unit_vector(v):\n", " return v / np.linalg.norm(v)\n", "```\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "7be327a8", "metadata": {}, "source": [ "### 💪 Exercise 1.2: Cosine similarity\n", "\n", "Write `cosine_similarity(a, b)` that returns the cosine similarity between\n", "two 1-D vectors, as defined above.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "49a511ef", "metadata": {}, "outputs": [], "source": [ "def cosine_similarity(a, b):\n", " # YOUR CODE HERE\n", " raise NotImplementedError\n" ] }, { "cell_type": "code", "execution_count": null, "id": "13255fb2", "metadata": {}, "outputs": [], "source": [ "same_dir = cosine_similarity(np.array([1.0, 0.0]), np.array([5.0, 0.0]))\n", "print(\"cosine_similarity([1, 0], [5, 0]) ->\")\n", "print(same_dir)\n", "\n", "orthogonal = cosine_similarity(np.array([1.0, 0.0]), np.array([0.0, 1.0]))\n", "opposite = cosine_similarity(np.array([1.0, 0.0]), np.array([-1.0, 0.0]))\n", "print(\"cosine_similarity([1, 0], [0, 1]) ->\", orthogonal)\n", "print(\"cosine_similarity([1, 0], [-1, 0]) ->\", opposite)\n", "\n", "_passed = (\n", " np.allclose(same_dir, 1.0)\n", " and np.allclose(orthogonal, 0.0)\n", " and np.allclose(opposite, -1.0)\n", ")\n", "\n", "if not _passed:\n", " print()\n", " print(\"Expected: 1.0, 0.0, and -1.0 respectively\")\n" ] }, { "cell_type": "markdown", "id": "c843242b", "metadata": {}, "source": [ "
\n", "💡 Hint (click to expand)\n", "\n", "`(a @ b) / (np.linalg.norm(a) * np.linalg.norm(b))` — just transcribe the formula above.\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "c64283d8", "metadata": {}, "source": [ "
\n", "✅ Solution (click to expand)\n", "\n", "```python\n", "def cosine_similarity(a, b):\n", " return (a @ b) / (np.linalg.norm(a) * np.linalg.norm(b))\n", "```\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "341d017f", "metadata": {}, "source": [ "## Lesson 2 — Matrix Multiplication\n", "\n", "For 2-D arrays, `@` performs real **matrix multiplication**: given\n", "`A` of shape `(m, n)` and `B` of shape `(n, p)`, `A @ B` has shape\n", "`(m, p)`, where each entry is the dot product of a row of `A` and a\n", "column of `B`. The inner dimensions (`n`) must match.\n", "\n", "This is completely different from `A * B`, which is the **elementwise**\n", "(Hadamard) product and requires `A` and `B` to have the same shape (or be\n", "broadcastable) — there's no \"inner dimension\" concept at all.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "0f147678", "metadata": {}, "outputs": [], "source": [ "A = np.array([[1, 2], [3, 4], [5, 6]]) # shape (3, 2)\n", "B = np.array([[1, 0], [0, 1]]) # shape (2, 2)\n", "\n", "print(\"A @ B shape:\", (A @ B).shape) # (3, 2)\n", "print(A @ B)\n", "\n", "C = np.array([[1, 2], [3, 4], [5, 6]]) # same shape as A\n", "print(\"A * C (elementwise):\")\n", "print(A * C)\n" ] }, { "cell_type": "markdown", "id": "7b012889", "metadata": {}, "source": [ "### 💪 Exercise 2.1: Weighted row sums\n", "\n", "Write `weighted_sums(matrix, weights)` that takes a `(m, n)` matrix and a\n", "length-`n` `weights` vector, and returns a length-`m` vector where each\n", "entry is the weighted sum of the corresponding row. Use a single matrix\n", "multiplication — no loop.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c3333b19", "metadata": {}, "outputs": [], "source": [ "def weighted_sums(matrix, weights):\n", " # YOUR CODE HERE\n", " raise NotImplementedError\n" ] }, { "cell_type": "code", "execution_count": null, "id": "ba1dd10e", "metadata": {}, "outputs": [], "source": [ "m = np.array([[1, 2, 3], [4, 5, 6]])\n", "w = np.array([1, 0, 1])\n", "result = weighted_sums(m, w)\n", "print(\"weighted_sums(matrix, weights) ->\")\n", "print(result)\n", "\n", "expected = np.array([4, 10])\n", "_passed = np.allclose(result, expected)\n", "\n", "if not _passed:\n", " print()\n", " print(\"Expected:\")\n", " print(expected)\n" ] }, { "cell_type": "markdown", "id": "95580511", "metadata": {}, "source": [ "
\n", "💡 Hint (click to expand)\n", "\n", "`matrix @ weights` does exactly this: a `(m, n) @ (n,)` product gives a length-`m` vector.\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "0a3936df", "metadata": {}, "source": [ "
\n", "✅ Solution (click to expand)\n", "\n", "```python\n", "def weighted_sums(matrix, weights):\n", " return matrix @ weights\n", "```\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "78cf1207", "metadata": {}, "source": [ "### 💪 Exercise 2.2: Matmul vs. Hadamard\n", "\n", "Write `both_products(A, B)` that takes two same-shape square matrices and\n", "returns a tuple `(matmul_result, hadamard_result)`: the true matrix product\n", "`A @ B` and the elementwise product `A * B`.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "71ebe8bb", "metadata": {}, "outputs": [], "source": [ "def both_products(A, B):\n", " # YOUR CODE HERE\n", " raise NotImplementedError\n" ] }, { "cell_type": "code", "execution_count": null, "id": "fe09b01f", "metadata": {}, "outputs": [], "source": [ "A = np.array([[1, 2], [3, 4]])\n", "B = np.array([[5, 6], [7, 8]])\n", "matmul_result, hadamard_result = both_products(A, B)\n", "print(\"both_products(A, B) ->\")\n", "print(\"matmul: \", matmul_result.tolist())\n", "print(\"hadamard:\", hadamard_result.tolist())\n", "\n", "expected_matmul = A @ B\n", "expected_hadamard = A * B\n", "_passed = (\n", " np.array_equal(matmul_result, expected_matmul)\n", " and np.array_equal(hadamard_result, expected_hadamard)\n", " and not np.array_equal(matmul_result, hadamard_result)\n", ")\n", "\n", "if not _passed:\n", " print()\n", " print(\"Expected matmul: \", expected_matmul.tolist())\n", " print(\"Expected hadamard:\", expected_hadamard.tolist())\n" ] }, { "cell_type": "markdown", "id": "00c0983b", "metadata": {}, "source": [ "
\n", "💡 Hint (click to expand)\n", "\n", "Literally return `(A @ B, A * B)`.\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "493a6b64", "metadata": {}, "source": [ "
\n", "✅ Solution (click to expand)\n", "\n", "```python\n", "def both_products(A, B):\n", " return A @ B, A * B\n", "```\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "d76af64d", "metadata": {}, "source": [ "## Lesson 3 — Solving Linear Systems\n", "\n", "To solve `Ax = b` for `x`, use `np.linalg.solve(A, b)`. It is faster and\n", "numerically more stable than computing `np.linalg.inv(A) @ b` — computing\n", "a full inverse does unnecessary extra work and amplifies floating-point\n", "error. **Prefer `solve` over `inv` whenever you just need `x`.**\n", "\n", "`np.linalg.det(A)` computes the determinant. A matrix with `det(A) == 0`\n", "is **singular** — it has no unique inverse, and `np.linalg.solve` will\n", "raise a `LinAlgError` for it.\n", "\n", "That said, not every LAPACK build raises reliably for a singular matrix —\n", "some just hand back `nan`/`inf` values instead of erroring. `np.isfinite(x)`\n", "is the elementwise check for exactly this: it returns `False` wherever a\n", "value is `nan` or `inf`, and `True` everywhere else, so `np.all(np.isfinite(x))`\n", "tells you in one line whether an entire array/result is \"clean\" (no\n", "`nan`/`inf` anywhere).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "73387528", "metadata": {}, "outputs": [], "source": [ "A = np.array([[3.0, 1.0], [1.0, 2.0]])\n", "b = np.array([9.0, 8.0])\n", "\n", "x = np.linalg.solve(A, b)\n", "print(\"solution x:\", x)\n", "print(\"check A @ x == b:\", np.allclose(A @ x, b))\n", "print(\"det(A):\", np.linalg.det(A))\n", "\n", "singular = np.array([[1.0, 2.0], [2.0, 4.0]]) # second row = 2 * first row\n", "print(\"det(singular):\", np.linalg.det(singular))\n", "try:\n", " np.linalg.solve(singular, np.array([1.0, 2.0]))\n", "except np.linalg.LinAlgError as e:\n", " print(\"LinAlgError:\", e)\n" ] }, { "cell_type": "markdown", "id": "57f1ed4e", "metadata": {}, "source": [ "### 💪 Exercise 3.1: Solve or report singular\n", "\n", "Write `solve_or_none(A, b)` that returns the solution vector `x` to\n", "`Ax = b` using `np.linalg.solve`, or returns `None` if `A` is singular.\n", "Catch `LinAlgError` — but also treat a result containing `nan`/`inf` as\n", "singular: depending on the exact LAPACK build, `np.linalg.solve` doesn't\n", "always *raise* for a singular matrix, it sometimes just returns garbage\n", "values instead, so a robust solution checks both.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "cf38c828", "metadata": {}, "outputs": [], "source": [ "def solve_or_none(A, b):\n", " # YOUR CODE HERE\n", " raise NotImplementedError\n" ] }, { "cell_type": "code", "execution_count": null, "id": "3dfe2994", "metadata": {}, "outputs": [], "source": [ "A = np.array([[3.0, 1.0], [1.0, 2.0]])\n", "b = np.array([9.0, 8.0])\n", "result = solve_or_none(A, b)\n", "print(\"solve_or_none(A, b) ->\")\n", "print(result)\n", "\n", "ok_basic = result is not None and np.allclose(A @ result, b)\n", "\n", "singular = np.array([[1.0, 2.0], [2.0, 4.0]])\n", "singular_result = solve_or_none(singular, np.array([1.0, 2.0]))\n", "print()\n", "print(\"solve_or_none(singular_matrix, b) ->\", singular_result)\n", "\n", "_passed = ok_basic and singular_result is None\n", "\n", "if not _passed:\n", " print()\n", " print(\"Expected: a vector x with A @ x ≈ b for the first call, and None for the singular matrix\")\n" ] }, { "cell_type": "markdown", "id": "b8166759", "metadata": {}, "source": [ "
\n", "💡 Hint (click to expand)\n", "\n", "Wrap `np.linalg.solve(A, b)` in a `try/except np.linalg.LinAlgError` and return `None` in the `except` branch — then also check `np.all(np.isfinite(x))` on whatever it returns, since a singular matrix doesn't always raise.\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "af76dfa0", "metadata": {}, "source": [ "
\n", "✅ Solution (click to expand)\n", "\n", "```python\n", "def solve_or_none(A, b):\n", " try:\n", " x = np.linalg.solve(A, b)\n", " except np.linalg.LinAlgError:\n", " return None\n", " return x if np.all(np.isfinite(x)) else None\n", "```\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "eb6d9aaf", "metadata": {}, "source": [ "### 💪 Exercise 3.2: Is it singular?\n", "\n", "Write `is_singular(A, tol=1e-10)` that returns `True` if the square matrix\n", "`A` is singular (its determinant is within `tol` of zero), `False`\n", "otherwise.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "61fa15e4", "metadata": {}, "outputs": [], "source": [ "def is_singular(A, tol=1e-10):\n", " # YOUR CODE HERE\n", " raise NotImplementedError\n" ] }, { "cell_type": "code", "execution_count": null, "id": "cce51aa2", "metadata": {}, "outputs": [], "source": [ "result_1 = is_singular(np.array([[1.0, 2.0], [2.0, 4.0]]))\n", "result_2 = is_singular(np.array([[3.0, 1.0], [1.0, 2.0]]))\n", "print(\"is_singular(singular_matrix) ->\", result_1)\n", "print(\"is_singular(regular_matrix) ->\", result_2)\n", "\n", "_passed = bool(result_1) and not bool(result_2)\n", "\n", "if not _passed:\n", " print()\n", " print(\"Expected: True for the singular matrix, False for the regular one\")\n" ] }, { "cell_type": "markdown", "id": "6ad9fb6c", "metadata": {}, "source": [ "
\n", "💡 Hint (click to expand)\n", "\n", "Compute `np.linalg.det(A)` and compare `abs(det)` against `tol`.\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "a5376333", "metadata": {}, "source": [ "
\n", "✅ Solution (click to expand)\n", "\n", "```python\n", "def is_singular(A, tol=1e-10):\n", " return abs(np.linalg.det(A)) < tol\n", "```\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "0a646265", "metadata": {}, "source": [ "## Lesson 4 — Eigenvalues & Eigenvectors\n", "\n", "For a square matrix `A`, an **eigenvector** `v` is a (nonzero) vector whose\n", "direction doesn't change when `A` is applied to it — only its length\n", "scales, by a factor called the **eigenvalue** `λ`:\n", "\n", "```\n", "A @ v == λ * v\n", "```\n", "\n", "`np.linalg.eig(A)` returns `(eigenvalues, eigenvectors)`, where\n", "`eigenvectors[:, i]` is the eigenvector for `eigenvalues[i]`.\n", "\n", "For non-square matrices there's no eigendecomposition, but there is an\n", "analogous, even more general tool: the **Singular Value Decomposition**\n", "(`np.linalg.svd`), which factors *any* matrix into rotation/scaling\n", "components — useful for things like dimensionality reduction, but outside\n", "the scope of this lesson's exercises.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "fa940e55", "metadata": {}, "outputs": [], "source": [ "A = np.array([[2.0, 0.0], [0.0, 3.0]])\n", "eigenvalues, eigenvectors = np.linalg.eig(A)\n", "print(\"eigenvalues: \", eigenvalues)\n", "print(\"eigenvectors:\\n\", eigenvectors)\n", "\n", "v0 = eigenvectors[:, 0]\n", "lam0 = eigenvalues[0]\n", "print(\"A @ v0: \", A @ v0)\n", "print(\"lam0 * v0:\", lam0 * v0)\n", "print(\"match:\", np.allclose(A @ v0, lam0 * v0))\n" ] }, { "cell_type": "markdown", "id": "09ae1481", "metadata": {}, "source": [ "### 💪 Exercise 4.1: Verify an eigenpair\n", "\n", "Write `is_eigenpair(A, v, lam, tol=1e-8)` that returns `True` if\n", "`A @ v` is within `tol` of `lam * v` (elementwise, via `np.allclose`),\n", "`False` otherwise.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "44e7b247", "metadata": {}, "outputs": [], "source": [ "def is_eigenpair(A, v, lam, tol=1e-8):\n", " # YOUR CODE HERE\n", " raise NotImplementedError\n" ] }, { "cell_type": "code", "execution_count": null, "id": "80e6f431", "metadata": {}, "outputs": [], "source": [ "A = np.array([[2.0, 0.0], [0.0, 3.0]])\n", "eigenvalues, eigenvectors = np.linalg.eig(A)\n", "\n", "result_1 = is_eigenpair(A, eigenvectors[:, 0], eigenvalues[0])\n", "print(\"is_eigenpair(A, eigenvectors[:, 0], eigenvalues[0]) ->\")\n", "print(result_1)\n", "\n", "result_2 = is_eigenpair(A, eigenvectors[:, 1], eigenvalues[1])\n", "result_3 = is_eigenpair(A, eigenvectors[:, 0], eigenvalues[1])\n", "print(\"is_eigenpair(A, eigenvectors[:, 1], eigenvalues[1]) ->\", result_2)\n", "print(\"is_eigenpair(A, eigenvectors[:, 0], eigenvalues[1]) (mismatched) ->\", result_3)\n", "\n", "_passed = bool(result_1) and bool(result_2) and not result_3\n", "\n", "if not _passed:\n", " print()\n", " print(\"Expected: True, True, False\")\n" ] }, { "cell_type": "markdown", "id": "289daf71", "metadata": {}, "source": [ "
\n", "💡 Hint (click to expand)\n", "\n", "`np.allclose(A @ v, lam * v, atol=tol)`.\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "52264783", "metadata": {}, "source": [ "
\n", "✅ Solution (click to expand)\n", "\n", "```python\n", "def is_eigenpair(A, v, lam, tol=1e-8):\n", " return np.allclose(A @ v, lam * v, atol=tol)\n", "```\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "2264c9bc", "metadata": {}, "source": [ "### 💪 Exercise 4.2: Dominant eigenvalue\n", "\n", "Write `dominant_eigenvalue(A)` that returns the eigenvalue of `A` with the\n", "largest absolute value, as a real float (eigenvalues from `np.linalg.eig`\n", "come back as complex numbers even when the imaginary part is zero — use\n", "`.real` after picking the right one, or take `abs()` for comparison and\n", "return the original real part).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "0c3760c3", "metadata": {}, "outputs": [], "source": [ "def dominant_eigenvalue(A):\n", " # YOUR CODE HERE\n", " raise NotImplementedError\n" ] }, { "cell_type": "code", "execution_count": null, "id": "b05e3619", "metadata": {}, "outputs": [], "source": [ "A = np.array([[2.0, 0.0], [0.0, 3.0]])\n", "result = dominant_eigenvalue(A)\n", "print(\"dominant_eigenvalue(A) ->\")\n", "print(result)\n", "\n", "expected = 3.0\n", "_passed = np.allclose(result, expected) and (isinstance(result, float) or np.isrealobj(result))\n", "\n", "if not _passed:\n", " print()\n", " print(\"Expected:\")\n", " print(expected)\n" ] }, { "cell_type": "markdown", "id": "28214549", "metadata": {}, "source": [ "
\n", "💡 Hint (click to expand)\n", "\n", "Get `eigenvalues` from `np.linalg.eig(A)`, use `np.argmax(np.abs(eigenvalues))` to find the index of the largest-magnitude one, then take `.real`.\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "9acb8423", "metadata": {}, "source": [ "
\n", "✅ Solution (click to expand)\n", "\n", "```python\n", "def dominant_eigenvalue(A):\n", " eigenvalues, _ = np.linalg.eig(A)\n", " idx = np.argmax(np.abs(eigenvalues))\n", " return eigenvalues[idx].real\n", "```\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "9b942c2a", "metadata": {}, "source": [ "## Lesson 5 — Random Number Generation with the Generator API\n", "\n", "Modern NumPy code should create an explicit `Generator` instead of using\n", "the legacy global functions (`np.random.seed`, `np.random.rand`, ...):\n", "\n", "```python\n", "rng = np.random.default_rng(seed=42)\n", "```\n", "\n", "A `Generator` gives you **reproducibility** (same seed -> same sequence)\n", "*and* **independent streams** (two `default_rng()` calls never silently\n", "share or clobber global state, which matters once you have multiple\n", "pieces of code sampling randomness at once). The most common methods:\n", "\n", "| Method | Produces |\n", "|---|---|\n", "| `rng.random(size)` | floats uniform in `[0, 1)` |\n", "| `rng.integers(low, high, size)` | integers in `[low, high)` |\n", "| `rng.normal(loc, scale, size)` | samples from a normal distribution |\n", "| `rng.choice(a, size, replace=...)` | random samples from an array |\n", "\n", "For `rng.choice`, `replace=True` (the default) can pick the same element\n", "more than once; `replace=False` samples without repeats, like dealing\n", "cards from a deck instead of rolling the same die repeatedly.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "e3172063", "metadata": {}, "outputs": [], "source": [ "rng = np.random.default_rng(seed=42)\n", "print(rng.random(3))\n", "print(rng.integers(0, 10, size=5))\n", "print(rng.normal(loc=0.0, scale=1.0, size=3))\n", "\n", "# same seed -> same sequence, every time\n", "rng_a = np.random.default_rng(seed=7)\n", "rng_b = np.random.default_rng(seed=7)\n", "print(np.array_equal(rng_a.random(5), rng_b.random(5)))\n", "\n", "samples = np.random.default_rng(0).normal(size=10_000)\n", "plt.hist(samples, bins=30)\n", "plt.title(\"10,000 samples from a standard normal distribution\")\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "bfa052a2", "metadata": {}, "source": [ "### 💪 Exercise 5.1: Reproducible dice rolls\n", "\n", "Write `roll_dice(n_rolls, seed)` that returns an array of `n_rolls`\n", "simulated six-sided die rolls (integers `1`-`6` inclusive), using\n", "`np.random.default_rng(seed)`. The same `seed` must always produce the\n", "same sequence.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d78d04b5", "metadata": {}, "outputs": [], "source": [ "def roll_dice(n_rolls, seed):\n", " # YOUR CODE HERE\n", " raise NotImplementedError\n" ] }, { "cell_type": "code", "execution_count": null, "id": "71d3bdba", "metadata": {}, "outputs": [], "source": [ "rolls = roll_dice(1000, seed=1)\n", "print(\"roll_dice(1000, seed=1) -> (first 20 of 1000)\")\n", "print(rolls[:20])\n", "\n", "ok_shape = rolls.shape == (1000,)\n", "ok_range = rolls.min() >= 1 and rolls.max() <= 6\n", "ok_reproducible = np.array_equal(roll_dice(10, seed=5), roll_dice(10, seed=5))\n", "\n", "_passed = ok_shape and ok_range and ok_reproducible\n", "\n", "if not _passed:\n", " print()\n", " print(\"Expected: shape (1000,), every value in 1..6, and the same seed always giving the same rolls\")\n", " print(f\"(shape={rolls.shape}, min={rolls.min()}, max={rolls.max()}, reproducible={ok_reproducible})\")\n" ] }, { "cell_type": "markdown", "id": "18455893", "metadata": {}, "source": [ "
\n", "💡 Hint (click to expand)\n", "\n", "`np.random.default_rng(seed).integers(1, 7, size=n_rolls)` — remember `integers`'s upper bound is exclusive, so use `7` to include `6`.\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "67c37bbf", "metadata": {}, "source": [ "
\n", "✅ Solution (click to expand)\n", "\n", "```python\n", "def roll_dice(n_rolls, seed):\n", " rng = np.random.default_rng(seed)\n", " return rng.integers(1, 7, size=n_rolls)\n", "```\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "c27695c7", "metadata": {}, "source": [ "### 💪 Exercise 5.2: Shuffle without replacement\n", "\n", "Write `random_sample_no_replace(arr, k, seed)` that returns `k` elements\n", "randomly chosen from `arr` *without* replacement (no duplicates), using\n", "`rng.choice` with the right arguments.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "b1804db9", "metadata": {}, "outputs": [], "source": [ "def random_sample_no_replace(arr, k, seed):\n", " # YOUR CODE HERE\n", " raise NotImplementedError\n" ] }, { "cell_type": "code", "execution_count": null, "id": "83d13ad9", "metadata": {}, "outputs": [], "source": [ "arr = np.arange(20)\n", "result = random_sample_no_replace(arr, 5, seed=3)\n", "print(\"random_sample_no_replace(np.arange(20), 5, seed=3) ->\")\n", "print(result)\n", "\n", "ok_shape = result.shape == (5,)\n", "ok_unique = len(set(result.tolist())) == 5\n", "ok_subset = np.all(np.isin(result, arr))\n", "\n", "_passed = ok_shape and ok_unique and ok_subset\n", "\n", "if not _passed:\n", " print()\n", " print(\"Expected: 5 unique values, all drawn from the original array\")\n", " print(f\"(shape={result.shape}, unique={ok_unique}, subset={ok_subset})\")\n" ] }, { "cell_type": "markdown", "id": "f78b9f8c", "metadata": {}, "source": [ "
\n", "💡 Hint (click to expand)\n", "\n", "`np.random.default_rng(seed).choice(arr, size=k, replace=False)`.\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "ccd8eb57", "metadata": {}, "source": [ "
\n", "✅ Solution (click to expand)\n", "\n", "```python\n", "def random_sample_no_replace(arr, k, seed):\n", " rng = np.random.default_rng(seed)\n", " return rng.choice(arr, size=k, replace=False)\n", "```\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "cd2a35e2", "metadata": {}, "source": [ "## Lesson 6 — Simulation Case Study: Monte Carlo Estimate of π\n", "\n", "Here's a classic way to estimate π without any trigonometry: scatter\n", "random points uniformly inside a `2x2` square centered at the origin, and\n", "check what fraction land inside the unit circle inscribed in that square.\n", "\n", "The area of the circle is `π * r² = π` (since `r = 1`), and the area of\n", "the square is `2 * 2 = 4`. So:\n", "\n", "```\n", "fraction inside circle ≈ π / 4 => π ≈ 4 * fraction inside circle\n", "```\n", "\n", "The whole \"check what fraction land inside the circle\" step is just\n", "**boolean indexing** from Module 2 — no loop needed: generate `x` and `y`\n", "coordinate arrays, build the mask `x**2 + y**2 <= 1`, and take its mean\n", "(since a boolean array's mean IS the fraction of `True` values).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "296cd71c", "metadata": {}, "outputs": [], "source": [ "rng = np.random.default_rng(seed=0)\n", "n = 100_000\n", "x = rng.uniform(-1, 1, size=n)\n", "y = rng.uniform(-1, 1, size=n)\n", "\n", "inside = x**2 + y**2 <= 1\n", "pi_estimate = 4 * inside.mean()\n", "print(\"pi estimate:\", pi_estimate)\n", "print(\"true pi: \", np.pi)\n", "\n", "plt.figure(figsize=(4, 4))\n", "plt.scatter(x[inside][:2000], y[inside][:2000], s=1, label=\"inside\")\n", "plt.scatter(x[~inside][:2000], y[~inside][:2000], s=1, label=\"outside\")\n", "plt.gca().set_aspect(\"equal\")\n", "plt.legend()\n", "plt.title(f\"Monte Carlo π ≈ {pi_estimate:.4f}\")\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "11204a4d", "metadata": {}, "source": [ "### 💪 Exercise 6.1: estimate_pi(n_samples, seed)\n", "\n", "Write `estimate_pi(n_samples, seed)` that runs the Monte Carlo procedure\n", "above: sample `n_samples` uniform `(x, y)` points in `[-1, 1] x [-1, 1]`\n", "using `np.random.default_rng(seed)`, and return `4 * (fraction of points\n", "with x**2 + y**2 <= 1)` as a float. Use vectorized operations only — no\n", "Python loop over the samples.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "53dee652", "metadata": {}, "outputs": [], "source": [ "def estimate_pi(n_samples, seed):\n", " # YOUR CODE HERE\n", " raise NotImplementedError\n" ] }, { "cell_type": "code", "execution_count": null, "id": "e841aa54", "metadata": {}, "outputs": [], "source": [ "# With 200,000 samples and this fixed seed, the estimate is reliably\n", "# within 0.02 of the true value — verified when writing this check.\n", "result = estimate_pi(200_000, seed=42)\n", "print(\"estimate_pi(200_000, seed=42) ->\")\n", "print(result)\n", "print(\"true pi:\", np.pi)\n", "\n", "ok_close = abs(result - np.pi) < 0.02\n", "ok_reproducible = estimate_pi(1000, seed=1) == estimate_pi(1000, seed=1)\n", "\n", "_passed = ok_close and ok_reproducible\n", "\n", "if not _passed:\n", " print()\n", " print(f\"Expected: a value within 0.02 of {np.pi}, and the same seed always giving the same estimate\")\n" ] }, { "cell_type": "markdown", "id": "bae61d10", "metadata": {}, "source": [ "
\n", "💡 Hint (click to expand)\n", "\n", "Reuse the exact pattern from the example above: `rng.uniform(-1, 1, size=n_samples)` for `x` and `y`, build the `inside` mask, return `4 * inside.mean()`.\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "c874e456", "metadata": {}, "source": [ "
\n", "✅ Solution (click to expand)\n", "\n", "```python\n", "def estimate_pi(n_samples, seed):\n", " rng = np.random.default_rng(seed)\n", " x = rng.uniform(-1, 1, size=n_samples)\n", " y = rng.uniform(-1, 1, size=n_samples)\n", " inside = x**2 + y**2 <= 1\n", " return 4 * inside.mean()\n", "```\n", "\n", "
\n" ] }, { "cell_type": "markdown", "id": "cd4c20f9", "metadata": {}, "source": [ "## 🎯 Key Takeaways\n", "\n", "- `np.dot`/`@` on vectors gives a scalar; `np.linalg.norm` gives length;\n", " cosine similarity measures direction regardless of magnitude.\n", "- `@` on 2-D arrays is true matrix multiplication (`(m,n) @ (n,p) ->\n", " (m,p)`) — completely different from elementwise `*`, which needs\n", " matching/broadcastable shapes.\n", "- Solve `Ax = b` with `np.linalg.solve`, not by inverting `A` — it's\n", " faster and more numerically stable. A zero determinant means `A` is\n", " singular and has no unique solution.\n", "- `np.linalg.eig` gives eigenvalues/eigenvectors satisfying `A @ v == λ *\n", " v`; `np.linalg.svd` generalizes the idea to non-square matrices.\n", "- Use `np.random.default_rng(seed)` for reproducible, independent random\n", " streams — never the legacy global `np.random.seed`/`rand` functions.\n", "- Random sampling + vectorized boolean indexing (Module 2) is enough to\n", " build real simulations, like estimating π with Monte Carlo.\n", "\n", "Next up: **Module 6 — Performance, I/O & Mini-Projects**, where you'll\n", "learn how memory layout affects speed, how to save/load arrays, and put\n", "everything from this course together in a few capstone exercises.\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.x" } }, "nbformat": 4, "nbformat_minor": 5 }