Problem Set 1 — Advanced Data Analytics

HEC Lausanne, autumn 2026 · Topics: Python and NumPy, Monte Carlo simulation, probability, numerical differentiation, feature encoding, linear regression, constrained optimization

Based on course materials by Prof. Simon Scheidegger.

This problem set practises the content of the first lectures. Exercise 8 uses linear regression and is meant to be done after the week-3 lecture; everything else you can start right away.

How to work through this notebook. Cells marked TODO are yours to complete — code cells contain ... where your code goes, markdown cells say (your answer here). Everything runs top to bottom with no internet access and no absolute paths; the data ships in data/. When you are done, Kernel → Restart & Run All must execute cleanly to the end.

Nothing is submitted. These problem sets are practice; the solution is discussed in the exercise session. Because nobody grades them, the notebook has to tell you whether you got it right — that is what the "Certify it" parts are for: checks you run on your own code, with the expected outcome stated so you can see for yourself.

AI assistants. You may — and realistically will — work with an AI assistant (ChatGPT, Claude, Copilot). That is fine. Assistants are very good at producing code that runs and looks right; the skill that will carry your project is telling whether it is right. Treat every "Certify it" part as the question you should be asking the assistant's code — and your own.

import math
import time

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.optimize import minimize
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import PolynomialFeatures

# One random generator for the whole notebook, seeded so that results are
# reproducible: rerunning the notebook gives the same numbers.
rng = np.random.default_rng(0)

Exercise 1 — Monte Carlo estimation of $\pi$

Monte Carlo methods are a broad class of algorithms that rely on repeated random sampling to obtain numerical results. The classic first example is estimating $\pi$.

Draw $N$ random points $(x, y)$ uniformly in the unit square $[0,1]^2$. The quarter disc $x^2 + y^2 \le 1$ has area $\pi/4$, the square has area $1$, so for large $N$

$$\frac{\text{points inside the quarter disc}}{N} \;\approx\; \frac{\pi}{4}, \qquad\text{i.e.}\qquad \pi \approx 4 \cdot \frac{\#\{x^2 + y^2 \le 1\}}{N}.$$

The picture shows the idea with $2000$ points. Your code does not need to draw anything: generate the pairs, count how many satisfy $x^2 + y^2 \le 1$, multiply by four. The more samples, the more accurate the estimate.

Task. Implement estimate_pi(N) and print a table showing how the accuracy of the estimate improves with the number of samples, for $N = 10^2, 10^3, \dots, 10^6$.

def estimate_pi(N):
    # TODO: draw N points in [0,1]^2 with rng.uniform, count those with
    # x^2 + y^2 <= 1, return 4 * count / N.
    ...


# TODO: print a table  N | estimate | |estimate - pi|  for N = 1e2 ... 1e6
...

Certify it

Each sample is a Bernoulli trial with success probability $p = \pi/4$, so the estimate $4\hat p$ has standard error $\;\mathrm{SE}(N) = 4\sqrt{p(1-p)/N}$. A correct implementation is therefore within $4\,\mathrm{SE}(N)$ of $\pi$ essentially always (probability of a miss $\approx 6\cdot10^{-5}$ per row), and its error shrinks like $1/\sqrt{N}$ — ten times the samples buys you only $\sqrt{10} \approx 3.2$ times the accuracy.

Check: for every $N$ in your table, assert abs(estimate - pi) < 4 * SE(N).

# TODO: compute SE(N) = 4*sqrt(p*(1-p)/N) with p = pi/4 and assert that
# each estimate in your table is within 4*SE(N) of pi.
...

Exercise 2 — Flipping a coin

Write a program that simulates flipping a coin $N$ times. Print head or tail for each flip and count the number of heads. (Hint: draw r = rng.uniform() and call it head if r <= 0.5, or draw an integer from $\{1, 2\}$ with rng.integers(1, 3) and call it head if it is 1.)

Task. Run it for $N = 20$ (printing every flip), then for $N = 10^6$ (printing only the count and the fraction of heads).

def flip_coins(N, verbose=False):
    # TODO: simulate N flips; if verbose, print "head"/"tail" for each;
    # return the number of heads.
    ...


# TODO: run for N=20 with printing, then N=1_000_000 without.
...

Certify it

The number of heads in $N$ fair flips has mean $N/2$ and standard deviation $\sqrt{N}/2$. Check: for $N = 10^6$, assert that the fraction of heads is within $4 \cdot 0.5/\sqrt{N} = 0.002$ of $0.5$.

# TODO: assert |fraction of heads - 0.5| < 4 * 0.5 / sqrt(N)
...

Exercise 3 — Probabilities of rolling dice

1. You throw a die. What is the probability of getting a 6?
2. You throw a die four times in a row. What is the probability of getting a 6 all four times?
3. Suppose you have thrown the die three times with 6 coming up every time. What is the probability of getting a 6 in the fourth throw?
4. Suppose you have thrown the die 100 times and got a 6 every single time. What do you think about the probability of getting a 6 in the next throw?

Task. First answer the four questions from a theoretical (or common-sense) point of view in the markdown cell below. Then write functions that simulate cases 1, 2 and 3 with $N = 10^6$ trials each and compare with your theoretical answers.

TODO (your answer here)

N = 1_000_000


def prob_six(N):
    # TODO: fraction of N single throws that are a 6
    ...


def prob_four_sixes(N):
    # TODO: fraction of N four-throw sequences that are all 6
    ...


def prob_six_after_three_sixes(N):
    # TODO: among N four-throw sequences, restrict to those whose first
    # three throws are 6; fraction of those whose fourth throw is a 6
    ...


# TODO: print the three simulated probabilities next to the exact values
...

Certify it

A simulated frequency of an event with probability $p$ over $n$ trials has standard error $\sqrt{p(1-p)/n}$. Check: each simulated value is within $4$ standard errors of its exact value — for case 3, $n$ is the number of qualifying sequences (about $N/216 \approx 4600$), not $N$.

# TODO: assert |simulated - exact| < 4*sqrt(p(1-p)/n) for the three cases
...

Exercise 4 — Throwing two dice

Write a program that throws two dice a large number of times. Record the sum of the eyes each time and count how often each possible sum $s = 2, 3, \dots, 12$ appears — a dictionary with the sum as key and the count as value is convenient. Divide the counts by the number of trials to get the frequency of each sum.

Then compute the exact probabilities: set up all $6 \times 6$ possible outcomes and count how many of them have sum $s$.

Task. Print the frequencies and the exact probabilities side by side for $S = 10^6$ throws.

S = 1_000_000

# TODO: simulate S throws of two dice; build a dict {sum: frequency}
...

# TODO: enumerate the 36 outcomes; build a dict {sum: exact probability}
...

# TODO: print  s | frequency | exact  for s = 2..12
...

Certify it

Check: the largest gap between frequency and exact probability over $s = 2, \dots, 12$ is below $4\sqrt{p_{\max}(1-p_{\max})/S}$ with $p_{\max} = 1/6$ (the sum $7$), i.e. below about $0.0015$ for $S = 10^6$. Also check that both dictionaries sum to $1$.

# TODO: assert max_s |frequency[s] - exact[s]| < 4*sqrt((1/6)(5/6)/S)
#       and that sum(frequency.values()) == sum(exact.values()) == 1
...

Exercise 5 — Waves in a laboratory tank

In a laboratory experiment, waves are generated by the impact of a model slide into a wave tank (the experiment models a future tsunami in a fjord, triggered by rocks falling into the water). At a fixed location the surface elevation $\eta$ is measured at discrete times with an ultrasound gauge, giving a time series of vertical positions in metres: $\eta(t_0), \eta(t_1), \dots, \eta(t_n)$. There are 300 observations per second, so the time step between neighbouring measurements is $h = 1/300$ s.

Task. Write code that

1. defines the time step h in a parameter cell (a notebook has no command line; a clearly marked parameter cell is the notebook equivalent);
2. reads the $\eta$ values from data/gauge.dat into an array eta;
3. plots eta against time;
4. computes the surface velocity by the central difference $$v_i \approx \frac{\eta_{i+1} - \eta_{i-1}}{2h}, \qquad i = 1, \dots, n-1,$$ and plots it against time in a separate figure;
5. computes the surface acceleration by $$a_i \approx \frac{\eta_{i+1} - 2\eta_i + \eta_{i-1}}{h^2}, \qquad i = 1, \dots, n-1,$$ and plots it against time in a separate figure.

# --- parameters ---
h = 1 / 300  # seconds between measurements

# TODO: read data/gauge.dat into a NumPy array eta (np.loadtxt is enough)
...

# TODO: time axis t = i*h; plot eta vs t
...

# TODO: velocity and acceleration at the interior points i = 1..n-2, and
# their plots (note: they are defined at t[1:-1])
...

Certify it

NumPy ships reference implementations of the same stencils: np.gradient(eta, h) uses exactly the central difference at interior points, and np.diff(eta, 2) / h**2 is exactly the second difference. Check: your velocity equals np.gradient(eta, h)[1:-1] and your acceleration equals np.diff(eta, 2) / h**2, both to machine precision (np.allclose). Off-by-one slices and misplaced factors of $2$ or $h$ fail this instantly.

# TODO: assert np.allclose against np.gradient(eta, h)[1:-1] and
#       np.diff(eta, 2) / h**2
...

Exercise 6 — Differentiating noisy signals

The purpose of this exercise is to see what numerical differentiation does to a time series that contains measurement noise. The insight will help you interpret the velocity and acceleration plots from Exercise 5.

1. Compute the clean signal $$\bar\eta_i = A \sin\!\left(\frac{2\pi}{T} t_i\right), \qquad t_i = i\,\frac{T}{40}, \quad i = 0, \dots, 200,$$ with $A = 1$ and $T = 2\pi$. Plot $\bar\eta$ against $t$ and store the values in an array etabar.
2. Compute a noisy signal $\eta_i = \bar\eta_i + E_i$, where the $E_i$ are drawn from a normal distribution with mean $0$ and standard deviation $\sigma = 0.04A$. Plot $\eta$ as circles in the same figure as $\bar\eta$. Keep the $E_i$ in an array E for later.
3. Compute the first derivative of $\bar\eta$ by the central difference $$\frac{\bar\eta_{i+1} - \bar\eta_{i-1}}{2h}, \qquad i = 1, \dots, n-1,$$ store it in detabar and plot it. What is $h$ here? Read it off the definition of $t_i$.
4. Apply the same formula to the noise, $\;(E_{i+1} - E_{i-1})/2h$, store it in dE, and compute the mean and standard deviation of dE.
5. Plot detabar and detabar + dE. Use the standard deviations to explain the qualitative features of the two curves.
6. The second derivative of a signal is $$\frac{\eta_{i+1} - 2\eta_i + \eta_{i-1}}{h^2}, \qquad i = 1, \dots, n-1.$$ Apply it to etabar (→ d2etabar) and to E (→ d2E). Plot d2etabar and d2etabar + d2E. Compute the standard deviation of d2E and compare it with those of dE and E. Discuss the plot in the light of these standard deviations.

A = 1.0
T = 2 * math.pi
i = np.arange(201)
t = i * T / 40
h = ...  # TODO: the spacing between consecutive t_i
sigma = 0.04 * A


def first_derivative(arr, h):
    # TODO: central difference (arr[i+1] - arr[i-1]) / (2h) at the interior points
    ...


def second_derivative(arr, h):
    # TODO: (arr[i+1] - 2 arr[i] + arr[i-1]) / h^2 at the interior points
    ...


# TODO steps 1-2: etabar, E, eta, and the plot
...

# TODO steps 3-5: detabar, dE, their statistics, the plot
...

# TODO step 6: d2etabar, d2E, the plot, the standard deviations
...

TODO (your answer here)

Certify it

Two checks, one for the stencils and one for the noise statistics.

C1 — stencils against calculus. On the clean signal the derivatives are known exactly: $\bar\eta' = A\frac{2\pi}{T}\cos(\frac{2\pi}{T}t)$ and $\bar\eta'' = -A(\frac{2\pi}{T})^2\sin(\frac{2\pi}{T}t)$. Central differences have truncation error $O(h^2)$; with $h \approx 0.157$ that is a few $10^{-3}$. Check: the max absolute error of detabar and of d2etabar against the exact derivatives at $t_1, \dots, t_{n-1}$ is below $0.01$. A wrong $h$, a missing factor $2$, or a /2*h precedence slip all produce errors of order $1$ or larger.

C2 — noise amplification against theory. Check: std(dE) / (sigma/(sqrt(2) h)) and std(d2E) / (sqrt(6) sigma/h**2) are both between $0.8$ and $1.2$ (about 200 samples, so a 20 % band is generous).

# TODO C1: exact derivatives at t[1:-1]; assert max abs error < 0.01 for both
...

# TODO C2: assert 0.8 < std(dE)/theory1 < 1.2 and 0.8 < std(d2E)/theory2 < 1.2
...

Exercise 7 — Kinds of features

Consider the following list of household appliances.

Category	Price	kWh/year	Energy efficiency	Weight	Brand	Sales
Freezer	699	210	A	99	Mühle	27.96
Washing machine	399	302	B	50	Boss	99.75
Freezer	779	475	C	75	Nerd	33.558
Washing machine	989	705	E	125	Prolet	3.996
Tableware dryer	199	790	D	42	Boss	2.786
Tableware dryer	479	810	D	37	Mühle	21.076

(a) For each of the seven features, say whether it is nominal, ordinal, or numerical.

(b) Encode the nominal and ordinal features so that the table can be used in a linear regression. Build the table as a pandas.DataFrame, apply the encoding in code, and display the result.

TODO (your answer here)

appliances = pd.DataFrame({
    "category":   ["Freezer", "Washing machine", "Freezer", "Washing machine", "Tableware dryer", "Tableware dryer"],
    "price":      [699, 399, 779, 989, 199, 479],
    "kwh_year":   [210, 302, 475, 705, 790, 810],
    "efficiency": ["A", "B", "C", "E", "D", "D"],
    "weight":     [99, 50, 75, 125, 42, 37],
    "brand":      ["Mühle", "Boss", "Nerd", "Prolet", "Boss", "Mühle"],
    "sales":      [27.96, 99.75, 33.558, 3.996, 2.786, 21.076],
})

# TODO: map efficiency A..E -> 1..5; one-hot encode category and brand
# (pd.get_dummies with drop_first=True); display the encoded frame.
...

Certify it

Check: every column of the encoded table is numeric (encoded.dtypes contains no object/string columns), and the number of dummy columns is (number of categories − 1) + (number of brands − 1) $= 2 + 3 = 5$.

# TODO: assert all columns are numeric and there are exactly 5 dummy columns
...

Exercise 8 — Linear regression on housing data (after the week-3 lecture)

data/california_housing.csv contains the 1990 California census: one row per block group (20,640 rows), with these columns

Column	Meaning
MedInc	median income in the block group (tens of thousands of USD)
HouseAge	median house age
AveRooms, AveBedrms	average rooms / bedrooms per household
Population	block-group population
AveOccup	average household size
Latitude, Longitude	location
MedHouseVal	target: median house value (hundreds of thousands of USD)

(a) Look at the description of the columns and explain, for each, whether it is nominal, ordinal, or numerical. Is there any column you would hesitate to feed to a linear model as-is?

(b) Use the numerical features to fit linear models for MedHouseVal with scikit-learn. Choose the features yourself — a correlation heatmap and a few scatter plots are the tools. Split the data into a training and a test set (80/20, fixed seed) and report the regression coefficients, $R^2$ and the MSE on the test set. Do not only fit a linear model: also fit a model with second-order polynomial features and compare.

(c) Implement the linear regression yourself, using the closed form derived in the lecture, $$\mathbf{w} = (\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{X}^\top \mathbf{y},$$ with a column of ones in $\mathbf{X}$ for the intercept. Compare the parameters, $R^2$ and MSE with those from (b).

TODO (your answer here)

housing = pd.read_csv("data/california_housing.csv")
housing.describe().T

# TODO: correlation heatmap (sns.heatmap(housing.corr(), ...)); pick features
...

# TODO: train/test split (test_size=0.2, random_state=0); fit LinearRegression;
# print coefficients, test R^2, test MSE
...

# TODO: same with PolynomialFeatures(degree=2) on your features; compare
...

TODO (your answer here)

# TODO: build X_train_1 = [1, X_train] and solve the normal equations for w;
# predict on the test set; compute R^2 and MSE; compare with (b).
# Hint: np.linalg.solve(X.T @ X, X.T @ y) computes (X^T X)^{-1} X^T y
# without explicitly forming the inverse — same formula, better numerics.
...

Certify it

Check: your closed-form parameters agree with scikit-learn's to $10^{-6}$ (absolute), and $R^2$ and MSE agree to $10^{-9}$. If the intercept is off but the slopes agree, the column of ones is missing; if everything is off, the formula is applied to the wrong matrix or the wrong split.

# TODO: assert the closed-form parameters, R^2 and MSE match sklearn
...

Exercise 9 — Constrained optimization

Solve the following problem with a Python optimization toolkit of your choice (scipy.optimize.minimize with method="trust-constr" is a good default):

$$\begin{aligned} \min_{x}\quad & x_1 x_4 (x_1 + x_2 + x_3) + x_3 && \text{objective}\\ \text{s.t.}\quad & x_1 x_2 x_3 x_4 \ge 25 && \text{inequality constraint}\\ & x_1^2 + x_2^2 + x_3^2 + x_4^2 = 40 && \text{equality constraint}\\ & 1 \le x_1, x_2, x_3, x_4 \le 5 && \text{bounds}\\ & x^{(0)} = (1, 5, 5, 1) && \text{initial guess} \end{aligned}$$

Solve it twice:

• (i) providing the analytical gradient and Hessian of the objective (derive them by hand in the markdown cell);
• (ii) providing a first-order finite-difference implementation of the gradient and Hessian that you write yourself, i.e. $\partial_i f(x) \approx \bigl(f(x + h e_i) - f(x)\bigr)/h$ and the Hessian by differencing that gradient again.

Compare the results of (i) and (ii): the optimum found, the number of iterations, and the running time. What do you observe?

(The constraints' derivatives are simple; give them analytically in both runs so that the comparison isolates the objective's derivatives.)

TODO (gradient and Hessian of the objective, by hand)

$$\nabla f(x) = \begin{pmatrix} ? \\ ? \\ ? \\ ? \end{pmatrix}, \qquad \nabla^2 f(x) = \begin{pmatrix} ? & ? & ? & ? \\ ? & ? & ? & ? \\ ? & ? & ? & ? \\ ? & ? & ? & ? \end{pmatrix}$$

def f(x):
    return x[0] * x[3] * (x[0] + x[1] + x[2]) + x[2]


def grad_f(x):
    # TODO: analytical gradient (4,)
    ...


def hess_f(x):
    # TODO: analytical Hessian (4, 4)
    ...


def fd_grad(x, h=1e-6):
    # TODO: forward differences (f(x + h e_i) - f(x)) / h
    ...


def fd_hess(x, h=1e-4):
    # TODO: forward-difference the gradient: column i = (fd_grad(x + h e_i) - fd_grad(x)) / h
    ...


constraints = [
    {"type": "ineq", "fun": lambda x: np.prod(x) - 25,
     "jac": lambda x: np.array([x[1] * x[2] * x[3], x[0] * x[2] * x[3], x[0] * x[1] * x[3], x[0] * x[1] * x[2]])},
    {"type": "eq", "fun": lambda x: np.sum(x**2) - 40, "jac": lambda x: 2 * x},
]
bounds = [(1, 5)] * 4
x0 = np.array([1.0, 5.0, 5.0, 1.0])

# TODO: run minimize(...) twice (analytical, finite-difference), timing each;
# print x, f(x), number of iterations, running time.
...

TODO (your answer here)

Certify it

C1 — derivatives. At $x^{(0)}$, your finite-difference gradient must agree with the analytical one to about $10^{-5}$ and the finite-difference Hessian to about $10^{-3}$ (forward differences: error $\sim h$). A sign error or a swapped index shows up as a disagreement of order $1$.

C2 — the optimum. Both runs must satisfy the constraints and land on $f(x^*) = 17.014017$ to within $10^{-4}$, with $x^* = (1.0000, 4.7430, 3.8212, 1.3794)$ to within $10^{-3}$.

# TODO C1: assert max|fd_grad(x0) - grad_f(x0)| < 1e-5 and max|fd_hess(x0) - hess_f(x0)| < 1e-3
...

# TODO C2: for both results assert |f - 17.014017| < 1e-4, |x - x_star| < 1e-3,
#          prod(x) >= 25 - 1e-6, |sum(x^2) - 40| < 1e-6
...

Bonus — automatic differentiation

Deriving gradients and Hessians by hand is error-prone (that is what C1 is for). Modern tools do it exactly, by tracing the program: with JAX, jax.grad(f) and jax.hessian(f) return functions that evaluate the derivatives to machine precision — no $h$, no truncation error.

Task. Compute the gradient and Hessian of $f$ at $x^{(0)}$ with JAX and check them against your analytical ones (they should agree to $\sim 10^{-6}$ in single precision, $10^{-12}$ if you enable float64).

# TODO: import jax; jax.grad(f)(x0) and jax.hessian(f)(x0); compare with grad_f, hess_f
...