Challenges - Miscellaneous#

Challenges that require you to use a selection of the topics learned so far

Stick-Breaking Triangle Experiment#

You have a stick of total length 1.

  1. Drop the stick and let it break at a random point, chosen uniformly along its length.

  2. Take the longer of the two resulting pieces and drop it again, breaking it at a random point (uniformly along that piece).

  3. At the end, you have three pieces whose total length is still 1.

Write a Python program to estimate the probability that these three pieces can form a triangle.

Tip: To generate uniformly distributed random numbers, use random.uniform(a, b). Remember to import the random module.

Hint 1: Triangle Condition#

Three line segments with a total length of 1 can form a triangle if and only if no single piece is at least 0.5 units long.

This is equivalent to the triangle inequality.

Hint 2: Estimating Probabilities#

When an exact formula is difficult to find, we can approximate a probability by performing the same random experiment many times.

The estimated probability is:

number of successful trials / total number of trials

Try running your simulation with a large number of trials, such as 10.000 or 100.000, to get a more accurate estimate.

Expected Result#

If your simulation is working correctly, you should expect your estimated probability to be close to

\[2\ln(2) - 1 \approx 0.3863\]
[ ]:

from random import uniform

# write here

The Metropolis-Hastings Algorithm#

Context#

In many scientific and statistical applications, we need to draw random samples from a probability distribution \(p(x)\). However, standard sampling methods (like inverse transform sampling) fail when \(p(x)\) is high-dimensional or complex.

More crucially, in many real-world scenarios, such as Bayesian inference or statistical mechanics, we only know the target distribution up to an unknown normalization constant:

\[p(x) = \frac{\tilde{p}(x)}{Z}\]

where \(\tilde{p}(x)\) is an unnormalized density function that we can evaluate, but \(Z = \int \tilde{p}(x) \, dx\) is an intractable integral that is computationally impossible to calculate.

The Metropolis-Hastings (MH) algorithm is a Markov Chain Monte Carlo (MCMC) method designed specifically to solve this problem: it generates a sequence of samples drawn from \(p(x)\) without ever calculating :math:`Z`.

Basically, this algorithm made Bayesian statistics practical beyond conjugacy!

How the Algorithm Works#

Starting from an initial state \(x_0\), the algorithm iteratively constructs a sequence of states \(x_0, x_1, x_2, \dots\) using a chosen proposal distribution \(q(x' \mid x)\) (e.g., a Gaussian distribution centered at the current state \(x\)).

For each step \(t = 1, 2, \dots, N\):

  1. Propose a candidate state: Draw a candidate \(x'\) from the proposal distribution \(q(x' \mid x_{t-1})\).

  2. Calculate the acceptance probability :math:`alpha`:

    \[\alpha = \min\left(1, \, \frac{\tilde{p}(x') \, q(x_{t-1} \mid x')}{\tilde{p}(x_{t-1}) \, q(x' \mid x_{t-1})}\right)\]

    (Note: If the proposal distribution is symmetric, :math:`q(x’ mid x) = q(x mid x’)`, the ratio simplifies to :math:`frac{tilde{p}(x’)}{tilde{p}(x_{t-1})}`)

  3. Accept or Reject:

    • Generate a random number \(u \sim \text{Uniform}(0, 1)\).

    • If \(u \le \alpha\), accept the candidate: set \(x_t = x'\).

    • Otherwise, reject the candidate: set \(x_t = x_{t-1}\).

Over time, the stationary distribution of this generated Markov chain converges to \(p(x)\).

Exercise: Implementation#

Write a Python function metropolis_hastings that implements the Metropolis-Hastings algorithm for a 1D target distribution.

Function Requirements#

  1. Parameters:

    • target_density: A function f(x) proportional to the target PDF \(\tilde{p}(x)\).

    • x0: Initial state (float).

    • n_samples: Number of iterations to run (int).

    • proposal_std: Standard deviation \(\sigma\) of a symmetric Gaussian proposal \(q(x' \mid x) = \mathcal{N}(x, \sigma^2)\) (default: 1.0).

  2. Return Value:

    • An array of length n_samples containing the sampled chain \(x_0, x_1, \dots, x_{N-1}\).

    • An integer representing the total count of accepted proposals (to evaluate the acceptance rate).

  3. Task:

    • Implement the function below.

    • Test your function by sampling from a target distribution \(\tilde{p}(x) = \exp(-x^2) \cdot (1 + \sin(3x)^2)\) using \(x_0 = 0.0\) and \(n\_samples = 20.000\).

[ ]:

import numpy as np
import matplotlib.pyplot as plt


# this is an example of an unnormalized pdf (i.e., its integral is not 1),
# you can change this to be any unnormalized pdf you want (but in that case if you want to
# plot your empirical result with the true distribution p
# you also need to change the normalisation constant in true_pdf function -> 2.65857)
def unnormalized_target(x):
    return np.exp(-x**2) * (1 + np.sin(3 * x)**2)

def true_pdf(x):
    return unnormalized_target(x) / 2.65857



def metropolis_hastings(target_density, x0, n_samples, proposal_std=1.0):
    raise Exception("TODO Implement me!")


samples, n_accepted = metropolis_hastings(
    target_density=unnormalized_target,
    x0=0.0,
    n_samples=200_000,
    proposal_std=1.0
)

# this is for plotting purposes, do not touch

x_grid = np.linspace(-4, 4, 1000)
y_true = true_pdf(x_grid)

plt.figure(figsize=(10, 5), dpi=100)
plt.hist(
    samples, # this should be returned by metropolis_hastings function!!
    bins=70,
    density=True,
    alpha=0.6,
    color='skyblue',
    edgecolor='navy',
    label='MH MCMC Samples'
)
plt.plot(
    x_grid,
    y_true,
    color='crimson',
    linewidth=2.5,
    label='True Normalized PDF $p(x)$'
)

plt.xlabel('$x$', fontsize=11)
plt.ylabel('Density', fontsize=11)
plt.legend(fontsize=11)
plt.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

\(k\)-Nearest Neighbors (\(k\)-NN) from Scratch#

Context#

The :math:`k`-Nearest Neighbors (:math:`k`-NN) algorithm is one of the foundational algorithms in supervised machine learning. Unlike parametric models (such as linear regression or logistic regression) that learn a set of parameters \(\boldsymbol{\theta}\) during a training phase, \(k\)-NN is a non-parametric, instance-based (or “lazy”) learning algorithm.

It makes no explicit assumptions about the functional form of the underlying data distribution and defers all computations until prediction time.

How the Algorithm Works#

Given a training dataset of \(N\) samples \(\mathbf{X}_{\text{train}} \in \mathbb{R}^{N \times d}\) with labels \(\mathbf{y}_{\text{train}} \in \{0, 1, \dots, C-1\}\), and a new query point \(\mathbf{x}_* \in \mathbb{R}^d\):

  1. Distance Calculation: Compute the distance between \(\mathbf{x}_*\) and every training point \(\mathbf{x}_i \in \mathbf{X}_{\text{train}}\). The standard choice is the Euclidean distance (\(L_2\) norm):

    \[d(\mathbf{x}_*, \mathbf{x}_i) = \Vert\mathbf{x}_* - \mathbf{x}_i\Vert_2 = \sqrt{\sum_{j=1}^{d} (x_{*,j} - x_{i,j})^2}\]

  2. Neighbor Identification: Sort the calculated distances in ascending order and select the indices corresponding to the \(k\) smallest distances, denoted by the neighborhood set \(\mathcal{N}_k(\mathbf{x}_*)\).

  3. Majority Voting (Classification): Assign \(\mathbf{x}_*\) the label that appears most frequently among its \(k\) nearest neighbors:

    \[\hat{y}_* = \arg\max_{c} \sum_{i \in \mathcal{N}_k(\mathbf{x}_*)} \mathbb{I}(y_i = c)\]

    where \(\mathbb{I}(\cdot)\) is the indicator function (evaluates to \(1\) if true, \(0\) otherwise).

Exercise: Implementation#

Implement clean, modular Python functions to perform \(k\)-NN classification.

Requirements#

  1. Function Structure:

    • compute_distances(X_train, X_test): Computes the \(M \times N\) Euclidean distance matrix between test queries \(\mathbf{X}_{\text{test}} \in \mathbb{R}^{M \times d}\) and training points \(\mathbf{X}_{\text{train}} \in \mathbb{R}^{N \times d}\).

    • knn_predict(X_train, y_train, X_test, k=3): Accepts training features, training labels, test queries, and the hyperparameter \(k\). Returns a NumPy array of predicted class labels for X_test.

  2. Task:

    • Complete the implementation of compute_distances and knn_predict below.

    • Evaluate your functions on the synthetic multi-class dataset provided in the test block.

    • Compute the accuracy score and compare it against scikit-learn’s KNeighborsClassifier.

Dataset Provided#

A synthetic dataset representing a multi-class classification problem is generated for you in the testing code below. You will use this fictional dataset to test your functions

[ ]:

import numpy as np

# X_train, X_test, and y_train are numpy arrays.
# in python, you can iterate over numpy arrays in the same way as lists,
# so you can use regular for loops. However, you can also use numpy
# functions and operations if you want!

def compute_distances(X_train, X_test):
    raise Exception("TODO: Implement me")


def knn_predict(X_train, y_train, X_test, k=3):
       raise Exception("TODO: Implement me")

[ ]:

k_value = 5
n_classes = 2
n_samples = 500

# this below is for testing your code, do not touch it

# if you dont have sklearn installed, type this in terminal:
# pip install -U scikit-learn

from sklearn.datasets import make_classification
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

X, y = make_classification(
    n_samples=n_samples,
    n_features=4,
    n_informative=4,
    n_redundant=0,
    n_classes=n_classes,
    random_state=42
)

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=0
)



print(f"Dataset summary: {X_train.shape[0]} train samples, {X_test.shape[0]} test samples.")



y_pred_custom = knn_predict(X_train, y_train, X_test, k=k_value)
custom_acc = accuracy_score(y_test, y_pred_custom)
print(f"\nAccuracy (k={k_value}): {custom_acc:.2%}")