diff --git a/exercises/solved_notebooks/P4_probmod/probmod_1-intro.jl b/exercises/solved_notebooks/P4_probmod/probmod_1-intro.jl
index 5f0b0b0..411d1f4 100644
--- a/exercises/solved_notebooks/P4_probmod/probmod_1-intro.jl
+++ b/exercises/solved_notebooks/P4_probmod/probmod_1-intro.jl
@@ -78,9 +78,9 @@ x(t) &= \mathrm{cos}(α) \, v_0 \, t \, ,
 ```
 
 where the parameters are:
--  $v_0$ is the initial velocity.
--  $\alpha$ is the launching angle.
--  $g = 9.8$ is the gravitational acceleration.
+-  $v_0$: the initial velocity.
+-  $\alpha$: the launching angle.
+-  $g = 9.8$: the gravitational acceleration.
 """
 
 # ╔═╡ ee9cd05c-7fa5-4412-b3b4-32b0a1774c87
@@ -91,11 +91,11 @@ t_f = \frac{x_f}{\mathrm{cos}(\alpha \, v_0)} \, ,
 ```
 evaluating it's y-coordinate there
 ```math
-y_d = \mathrm{sin}(α) \, v_0 \, t_f - \frac{g}{2} \, t_f^2 \, ,
+y_f = \mathrm{sin}(α) \, v_0 \, t_f - \frac{g}{2} \, t_f^2 \, ,
 ```
 and checking if it's in the range of our face
 ```math
-y_d ∈ [1.5, 1.7] \, .
+y_f ∈ [1.5, 1.7] \, .
 ```
 """
 
@@ -143,7 +143,13 @@ plot(
 )
 
 # ╔═╡ 314ba0ec-015f-488d-ac63-780b1dbf5bc2
-md"The latter is an **input variable**: when cycling, we can choose how far away we stay from the cyclist in front of us. The x-position of our face $x_f$ is therefore an input to the model."
+md"The third is an **input variable**: when cycling, we can choose how far away we stay from the cyclist in front of us. The x-position of our face $x_f$ is therefore an input to the model."
+
+# ╔═╡ 2878b32a-9e54-40b0-b151-031b95ad161b
+md"""
+!!! note
+	As you can see, to visualize a distribution you can simply call the `plot` function on the distribution (and set extra keyword arguments as desired).
+"""
 
 # ╔═╡ c070f3c9-dd80-49f9-b279-c77f46316116
 md"## Defining the model"
@@ -182,11 +188,11 @@ Our model can be defined as follows:
 	# get the time it takes the droplet to travel the distance x_f
 	t_f = x_f / (cos(α) * v0)
 	# get the droplet's altitude at this time
-	y_s = sin(α)*v0*t_f - g/2 * t_f^2
+	y_f = sin(α)*v0*t_f - g/2 * t_f^2
 	# to respect the ground, set altitude to 0 if it would be negative
-	y_s = max(0, y_s)
+	y_f = max(0, y_f)
 	
-	return y_s 
+	return y_f
 end
 
 # ╔═╡ 22144610-0c37-4501-b14c-83895d8b0eae
@@ -235,22 +241,22 @@ md"""
 md"### Getting the model output"
 
 # ╔═╡ 46f17c07-9aec-4fee-8590-4eeb5a13b93e
-md"To answer our question, we are interested in the `y_s` variable returned by the model. We can extract the output variable from our chain using the `generated_quantities` function:"
+md"To answer our question, we are interested in the `y_f` variable returned by the model. We can extract the output variable from our chain using the `generated_quantities` function:"
 
 # ╔═╡ 071fc405-7a33-4dab-a69a-1e68a0fb35ad
-y_s_samples = generated_quantities(splash_model, chain)
+y_f_samples = generated_quantities(splash_model, chain)
 
 # ╔═╡ 20e99145-e80b-4259-9bca-f7e14c2d355e
 md"It's often useful to visualize the distribution of our variable of interest. We can use the `histogram` function for this:"
 
 # ╔═╡ 45d44cf5-2289-4947-8d4e-1dcddb5447de
-histogram(y_s_samples, bins = 15)
+histogram(y_f_samples, bins = 15)
 
 # ╔═╡ 8d179748-40aa-4e88-b277-9234fe938251
 md"Checking whether the droplet altitude matches our face's can be done with a simple inequality operation:"
 
 # ╔═╡ 633a2b4c-379e-4e18-bbcd-b450930196a5
-face_hit_samples = [1.5 <= y_s <= 1.7 for y_s in y_s_samples]
+face_hit_samples = [1.5 <= y_f <= 1.7 for y_f in y_f_samples]
 
 # ╔═╡ ad8a6ea2-30ba-4104-824f-12a7ea8e9718
 md"Finally, we can estimate the probability of our face being hit by taking the mean amount of times we got hit in our sample."
@@ -264,7 +270,7 @@ md"""And that's the answer to our question! If we're 1m behind the cyclist in fr
 # ╔═╡ bf5642f4-09d3-491d-8eb9-47fcdb67c977
 md"""
 !!! note
-	Do you like your code sleek? You can also calculate this probability from your samples with just one line: `mean(1.5 .<= y_s_samples .<= 1.7)`
+	Do you like your code sleek? You can also calculate this probability from your samples with just one line: `mean(1.5 .<= y_f_samples .<= 1.7)`
 """
 
 # ╔═╡ 3067af6e-df69-48a0-a024-a6e219dcfc7c
@@ -335,6 +341,49 @@ md"""
 	We set a lot of keyword arguments of our plot to make the figure pretty here, but don't panic: this is still no programming course, so we will always clearly provide which one(s) you need if you need to plot something. If you ever do want to look for one yourself, you can find them on the function's help page (`?plot` or the `🔍 Live docs` in the bottom right corner).
 """
 
+# ╔═╡ 1de58e6a-4ee2-4d67-b97d-e63bbb9ed057
+md"## Sampling alternative"
+
+# ╔═╡ 83a8b9b6-72da-4aff-8364-37846b527a2f
+md"It's also possible to skip the chain construction entirely and simply call your model as a function to get a sample of the model output. Note that it's not possible to get the corresponding values of the stochastic variables using this approach, so only use this method if you **only** care about the model output."
+
+# ╔═╡ 771b9cb1-d809-46de-9380-2db73ccd014d
+# get a sample of the output 500 times
+y_f_samples_alternative = [splash_model() for _ in 1:500]
+
+# ╔═╡ 0f1cccd8-937d-4fc2-8657-c9967db483b8
+histogram(y_f_samples_alternative, bins = 15) # yup, that's the same distribution
+
+# ╔═╡ 7806eca0-a2ad-4910-8d6e-edc50073b37b
+md"## For-loops for many variables"
+
+# ╔═╡ 64ebc3ce-83ca-4cee-ab93-cddf51138543
+md"""
+Sometimes you need to define a large number of stochastic variables following the same distribution. Let's say, for example, you want to model the total mass of water being splashed on your face assuming that 10 droplets will hit you, each having a mass that is exponentially distributed with a mean of 30 mg. 
+
+This can be done using a **for-loop** as follows:
+"""
+
+# ╔═╡ bdf3c18c-b40b-4fb4-81e7-fe9d821cc166
+@model function splash_mass()
+ 	droplet_masses = zeros(10) # instantiate vector of 10 zeros
+	for i in 1:length(droplet_masses)
+		droplet_masses[i] ~ Exponential(30) # each element of the vector follows an Exponential distribution
+	end
+	total_mass = sum(droplet_masses)
+	
+	return total_mass
+end
+
+# ╔═╡ 2cbe2750-38c1-4e7f-8d03-9bed7c16861a
+mass_model = splash_mass();
+
+# ╔═╡ ad8b48c5-f892-49b3-8f19-64d93424485c
+mass_chain = sample(mass_model, Prior(), 1000);
+
+# ╔═╡ 8507f084-eff5-4351-800d-8eba7ca673bb
+total_mass_samples = generated_quantities(mass_model, mass_chain)
+
 # ╔═╡ befb8f05-6ab7-4594-9cf1-71df47c1df03
 md"## Essentials"
 
@@ -358,33 +407,32 @@ md"""
 # ╠═╡ show_logs = false
 let
 	# model definition (comment out for speed)
-	g = 9.8
 	@model function splash(x_f)
 		v0 ~ Normal(5, 1) 
 		α ~ TriangularDist(0, pi/2)
 	
 		t_f = x_f / (cos(α) * v0)
-		y_s = sin(α)*v0*t_f - g/2 * t_f^2
-		y_s = max(0, y_s)
-		return y_s 
+		y_f = sin(α)*v0*t_f - g/2 * t_f^2
+		y_f = max(0, y_f)
+		return y_f
 	end
 
 	# get samples
 	splash_model = splash(x_f)
 	chain = sample(splash_model, Prior(), 500);
 
-	y_s_samples = generated_quantities(splash_model, chain)
+	y_f_samples = generated_quantities(splash_model, chain)
 	v0_samples = chain["v0"]
 	α_samples = chain[:α]
 
 	# calculate interesting things
-	prob_face_hit = mean([1.5 <= y_s <= 1.7 for y_s in y_s_samples])
+	prob_face_hit = mean([1.5 <= y_f <= 1.7 for y_f in y_f_samples])
+
+	y(x, α, v0) = sin(α)/cos(α)*x - g/2*(x / (cos(α) * v0))^2
 	trajectories = [
-		x -> sin(α_samples[i]) / cos(α_samples[i]) * x -
-		g / 2 * ( x / (cos(α_samples[i])*v0_samples[i]) )^2
+		x -> y(x, α_samples[i], v0_samples[i])
 		for i in 1:length(v0_samples)
 	]
-
 	plot(trajectories, xlims = (0, 5), ylims = (0, 2), 
 		 legend = false, color = :blue, alpha = 0.3,
 		 title = "Probability of getting hit ≈ $(round(100*prob_face_hit, digits = 3))%"
@@ -411,10 +459,11 @@ end
 # ╟─e8821ac2-baa3-4674-868c-3a5b8593cc95
 # ╟─4811e1af-9ba4-49b3-a3a1-5d8ff8d5af92
 # ╟─01156ef7-b194-4372-8bc8-de9fad5e12f1
-# ╟─f30290bd-9d31-4939-b5d6-f89f066775d2
+# ╠═f30290bd-9d31-4939-b5d6-f89f066775d2
 # ╟─dab23eba-3eb1-434b-aef5-e08730e35114
-# ╟─299e402c-634a-4b18-b867-489bafaf4564
+# ╠═299e402c-634a-4b18-b867-489bafaf4564
 # ╟─314ba0ec-015f-488d-ac63-780b1dbf5bc2
+# ╟─2878b32a-9e54-40b0-b151-031b95ad161b
 # ╟─c070f3c9-dd80-49f9-b279-c77f46316116
 # ╟─84a08354-7bf3-4783-9285-a7e8850e798b
 # ╟─8f94567c-b87b-4eb4-8441-c6434da393a8
@@ -442,7 +491,7 @@ end
 # ╟─ad8a6ea2-30ba-4104-824f-12a7ea8e9718
 # ╠═e1f9d1e9-8a0c-4b28-ab67-a89e3287b50f
 # ╟─1c9e2de1-21a0-4d4d-918d-9e2f0fc4d912
-# ╟─bf5642f4-09d3-491d-8eb9-47fcdb67c977
+# ╠═bf5642f4-09d3-491d-8eb9-47fcdb67c977
 # ╟─3067af6e-df69-48a0-a024-a6e219dcfc7c
 # ╟─cbc04c5f-aa93-4340-a61e-1583d20d86b1
 # ╟─d36784b8-9f3b-4434-88fe-03b75b94754e
@@ -457,6 +506,16 @@ end
 # ╟─00572e1a-d1e6-4c84-b1b2-b53404087dd8
 # ╠═f550ed7a-1bba-40b5-9338-cdeb0520aae9
 # ╟─3b465fc0-3a9d-4e67-8d4a-a14b3cd61f4e
+# ╟─1de58e6a-4ee2-4d67-b97d-e63bbb9ed057
+# ╟─83a8b9b6-72da-4aff-8364-37846b527a2f
+# ╠═771b9cb1-d809-46de-9380-2db73ccd014d
+# ╠═0f1cccd8-937d-4fc2-8657-c9967db483b8
+# ╟─7806eca0-a2ad-4910-8d6e-edc50073b37b
+# ╟─64ebc3ce-83ca-4cee-ab93-cddf51138543
+# ╠═bdf3c18c-b40b-4fb4-81e7-fe9d821cc166
+# ╠═2cbe2750-38c1-4e7f-8d03-9bed7c16861a
+# ╠═ad8b48c5-f892-49b3-8f19-64d93424485c
+# ╠═8507f084-eff5-4351-800d-8eba7ca673bb
 # ╟─befb8f05-6ab7-4594-9cf1-71df47c1df03
 # ╟─52af5de0-b516-4947-84f7-dde63e4c513e
 # ╟─2ceede7e-8f0c-420c-bde2-f16b726cd204
diff --git a/exercises/solved_notebooks/P4_probmod/probmod_2-basics.jl b/exercises/solved_notebooks/P4_probmod/probmod_2-basics.jl
index 5551d01..fc49a42 100644
--- a/exercises/solved_notebooks/P4_probmod/probmod_2-basics.jl
+++ b/exercises/solved_notebooks/P4_probmod/probmod_2-basics.jl
@@ -1,5 +1,5 @@
 ### A Pluto.jl notebook ###
-# v0.20.4
+# v0.20.21
 
 using Markdown
 using InteractiveUtils
@@ -49,13 +49,13 @@ dpmodel = doublepoisson()
 dpchain = sample(dpmodel, Prior(), 10_000)
 
 # ╔═╡ f6d90b53-0f16-4e43-b77c-8eb05b8c6943
-spY = dpchain[:Y];
+Y_samples = dpchain[:Y];
 
 # ╔═╡ 90835b5d-9705-43d0-b449-09cde57d5394
 plot(Poisson(10))
 
 # ╔═╡ 431d9be9-edfc-48fb-9df7-6a2b961fe4b4
-histogram(spY)
+histogram(Y_samples)
 
 # ╔═╡ 13989ba4-bcf8-4fdd-8aee-ab58c8905bc9
 md"### 2: Probabilities"
@@ -65,16 +65,16 @@ md"""
 !!! tip
 	When comparing a vector of values to a single number, don't forget to use `.` to execute operations element-wise in Julia!
 	
-	✅ `spY .< 1` compares every element of `spY` to `1`
+	✅ `Y_samples .< 1` compares every element of `Y_samples` to `1`
 	
-	❌ `spY < 1` compares an entire vector with a single number → errors :(
+	❌ `Y_samples < 1` compares an entire vector with a single number → errors :(
 """
 
 # ╔═╡ 09a87037-9f4f-4dd9-bb6c-fe717db39ea6
-probXY1 = mean(3 .< spY .<= 10)
+probXY1 = mean(3 .< Y_samples .<= 10)
 
 # ╔═╡ 2e2f203a-34a6-4feb-8aa1-8a1c97a09165
-probXY2 = mean(spY.^2 .> 100)
+probXY2 = mean(Y_samples.^2 .> 100)
 
 # ╔═╡ c243ca59-191d-4905-825e-6d7825a3c8a4
 md"### 3: Variances"
@@ -84,91 +84,32 @@ md"""
 !!! hint
 	To create a sample of $X$ that is conditional on some value(s) of $Y$, you can start from a sample of $X$ and select only those elements for which the corresponding sample of $Y$ has the conditioned value(s).
 
-	In other words, you'll need to index `spX` based on `spY` (and vice versa for $\text{var}(Y ∣ X)$). 
+	In other words, you'll need to index `X_samples` based on `Y_samples` (and vice versa for $\text{var}(Y ∣ X)$). 
 """
 
 # ╔═╡ 04cc5b42-7ab2-4055-a959-ba894c598f69
-spX = dpchain[:X];
+X_samples = dpchain[:X];
 
 # ╔═╡ 597d97fb-cfbe-4cff-bfbb-57df9a2f42aa
-varXcondY = var(spX[spY .== 15])
+varXcondY = var(X_samples[Y_samples .== 15])
 
 # ╔═╡ 6ca31280-a997-4e84-98bb-e7784d0c555c
-varYcondX = var(spY[spX .== 15])
+varYcondX = var(Y_samples[X_samples .== 15])
 
 # ╔═╡ 9dcb122b-dd12-4e80-b391-f780e637d6fe
 var(Poisson(15)) 
 	# analytical answer of var(Y|X=15)
 	# Given X, Y follows a Poisson distribution with mean X
 
-# ╔═╡ ce7d57ed-4f31-4dcf-af3b-b37a2e2a9393
-md"## 2: Combinations"
-
-# ╔═╡ 087994ce-a26d-40c4-87eb-ef9f0ce7f1fb
-md"""
-Let `U ~ Uniform(0, 4)`, `V ∼ Normal(U, 1)` and `W ~ TriangularDist(0, 4, U)`. 
-1. Use sampling (n = 10_000) to make a histogram of `|V − W|`.
-2. Estimate `P(V > W)` and `P(V * W >= 10)`
-3. Are `V` and `W` independent?
-"""
-
-# ╔═╡ 9e5cc347-c74f-46a3-9534-c5ad812844bf
-md"### 1: Histogram"
-
-# ╔═╡ 47a43282-3892-4a9a-94b7-c359fa74e12b
-@model function combinations()
-	U ~ Uniform(0, 4)
-	V ~ Normal(U, 1.0)
-	W ~ TriangularDist(0, 4, U)
-end
-
-# ╔═╡ 890eaeb5-c63f-470c-85eb-d77cec68b740
-comb_model = combinations()
-
-# ╔═╡ 0d0297a3-60b6-47e5-b42b-bde99c09a16e
-comb_chain = sample(comb_model, Prior(), 10_000);
-
-# ╔═╡ c3cec28a-f472-46f5-9082-ddfdeb01ddf0
-spV = comb_chain[:V];
-
-# ╔═╡ 1911f031-d668-47df-995a-ca339ad4ae47
-spW = comb_chain[:W];
-
-# ╔═╡ 5c13b506-c05b-4258-95e7-93b648defa17
-spVW = abs.(spV - spW)
-
-# ╔═╡ d8ae7b0a-afee-41c0-af0c-f8b1d30086f9
-histogram(spVW)
-
-# ╔═╡ faa105b6-0700-4f4d-92fe-2bb72a4d6e44
-md"### 2: Probabilities"
-
-# ╔═╡ b4f156df-0b9c-495a-b95b-000c7f166243
-probVW1 = mean(spV .> spW)
-
-# ╔═╡ f7eb6faf-8884-4fea-9974-330d7fd555aa
-probVW2 = mean(spV .* spW .>= 10)
-
-# ╔═╡ 4decd959-aeb9-47d4-a381-14bbf4dbc5ab
-md"### 3: Independence"
-
-# ╔═╡ a465722f-49f9-4d37-9a2b-00d1120f5c06
+# ╔═╡ c976dd77-6a72-4acc-9b07-342cc01fc7ee
 md"""
-!!! hint
-	One way to prove dependence is showing that $E[V] \neq E[V \mid W \leq w]$ for at least one value $w$. 
-
-	If the expected value of $V$ can change based on some information about $W$, they can't be independent!
+!!! warning
+	Starting here, the entire structure of the answer will no longer be given. You are therefore expected to **add more code cells** yourself, using the `+` at the left in between two cells.
 """
 
-# ╔═╡ b2999f6e-2dbd-44ab-a303-97234f665d33
-mean(spV)
-
-# ╔═╡ 5ea80b00-941a-4226-87cb-66da7ee7d76d
-mean(spV[spW .<= 1]) # E[V] != E[V | W <= 1] => they are dependent
-
 # ╔═╡ 9740ea64-cd4f-46b1-a741-02e392280601
 md"""
-## 3: Dice
+## 2: Dice
 """
 
 # ╔═╡ 187854bb-9e30-454d-9e03-cccf77aebb6b
@@ -192,7 +133,7 @@ md"### 1: Watercube"
 
 # ╔═╡ 13cdd6ee-04ec-4085-8083-3d52e88aa35c
 md"""
-!!! hint
+!!! tip
 	Consider the humble `DiscreteUniform` distribution. Not sure how it works? Open the **🔍 Live Docs** at the bottom right of the screen for more information
 """
 
@@ -212,13 +153,13 @@ end
 watermodel = watercube()
 
 # ╔═╡ 97021710-4b26-4a1c-a794-265c9559ce4d
-sp_w = [watermodel() for i in 1:2000]
+watercube_samples = [watermodel() for i in 1:2000]
 
 # ╔═╡ c0cd75bd-2a46-4f62-a59a-9bb8d47d45f2
-p_watercube_kills = mean(sp_w .>= 50)
+p_watercube_kills = mean(watercube_samples .>= 50)
 
 # ╔═╡ bc5f5d70-e269-45f9-867b-a00876ce8c40
-histogram(sp_w, bins = 30)
+histogram(watercube_samples, bins = 30)
 
 # ╔═╡ 44a2cabc-aca9-4bb7-af9d-cd735c424d7b
 md"### 2: Dirtprism"
@@ -237,32 +178,32 @@ end
 dirtmodel = dirtprism()
 
 # ╔═╡ 700cc645-1123-4c8a-821e-1ae4fe8b0674
-sp_d = [dirtmodel() for i in 1:2000]
+dirtcube_samples = [dirtmodel() for i in 1:2000]
 
 # ╔═╡ 98d04790-e86f-438a-9389-cb9697934bbf
-p_dirtprism_kills = mean(sp_d .>= 50)
+p_dirtprism_kills = mean(dirtcube_samples .>= 50)
 
 # ╔═╡ d5413573-5202-4fa2-94f6-e92a613d63aa
-histogram(sp_d, bins = 30)
+histogram(dirtcube_samples, bins = 30)
 
 # ╔═╡ 49790a8f-9f53-4ba7-9543-d6a879b520e0
 md"### 3: Comparison"
 
 # ╔═╡ 6e098cb6-eddd-4924-a184-c2578dc28473
-p_watercube_is_better = mean(sp_w .> sp_d)
+p_watercube_is_better = mean(watercube_samples .> dirtcube_samples)
 
 # ╔═╡ 34f3014f-f4d4-43d1-b46f-bdca73aee33f
-md"## 4: Super eggs"
+md"## 3: Super eggs"
 
 # ╔═╡ 372436c4-262f-49b8-b1cf-626b043542bf
 md"""
-When a chicken lays an egg, there's a small chance it contains two egg yolks. This chance, as well as the number of eggs a chicken lays per year, go down as the chicken gets older. 
+When a chicken lays an egg, there's a small chance it contains two egg yolks. This chance, as well as the number of eggs a chicken lays per year, goes down as the chicken gets older. 
 """
 
 # ╔═╡ 20111742-008a-44c3-8c27-62791cce3e1e
 md"""
 You can make the following assumptions
-- The age $A$ of a random chicken (in years) is discrete and Uniformly distributed between 0 and 12.
+- The age $A$ of a random chicken (in years) is discrete and Uniformly distributed between 0 and 10.
 - The number of eggs $N$ an $A$-year old chicken lays in a year is Poisson distributed with mean $300 - 20 \, A$.
 - The probability $P$ of an $A$-year old chicken's egg having a double yolk is distributed as a `Beta(1, 800 + 100*A)`.
 """
@@ -271,7 +212,7 @@ You can make the following assumptions
 md"""
 !!! questions
     1. If someone hands you a random chicken, what is the probability it will lay 2 or more double eggs in a year? 
-    1. Compare the distributions of double eggs for 1-year old and 3-year old chickens.
+    1. Compare the distributions of double eggs for 1-year old and 5-year old chickens.
 """
 
 # ╔═╡ 6ebe676b-aee6-459b-bcda-de0fe14418a3
@@ -279,7 +220,7 @@ md"### 1: Probability"
 
 # ╔═╡ 34aad3d6-9deb-4298-b613-bf3a616f2b31
 md"""
-!!! hint
+!!! tip
 	In this exercise, the output variable (the number of double-yolked eggs) is **also a random variable**! In other words, it also follows some distribution.
 
 	When considering what distribution, consider that each of the $N$ eggs represents a "trial" with a $P$ chance of success for a double yolk.
@@ -287,7 +228,7 @@ md"""
 
 # ╔═╡ 2b3d930f-53d9-4869-9e4a-86a1a681b9d8
 @model function eggs()
-	A ~ DiscreteUniform(0, 12)
+	A ~ DiscreteUniform(0, 10)
 	N ~ Poisson(300 - 20*A)
 	P ~ Beta(1, 800 + 100*A)
 
@@ -302,10 +243,10 @@ egg_model = eggs();
 chain_egg = sample(egg_model, Prior(), 2000)
 
 # ╔═╡ 64d04ed9-6548-4251-8a4e-11b31e8f3143
-sp_egg = generated_quantities(egg_model, chain_egg);
+egg_samples = generated_quantities(egg_model, chain_egg);
 
 # ╔═╡ 00fd5517-5232-4d88-8ab0-3a0eb925eb3e
-p_multiple_double_eggs = mean(sp_egg .>= 2) 
+p_multiple_double_eggs = mean(egg_samples .>= 2) 
 	# a ~2% chance for two or more double eggs in one year
 
 # ╔═╡ 80e3544e-86ec-469c-be44-77f94fabfc28
@@ -315,13 +256,13 @@ md"### 2: Histograms"
 chicken_ages = chain_egg[:A];
 
 # ╔═╡ fe87f9ee-a07c-4b4c-b89d-c3bb4ee7ee7f
-histogram(sp_egg[chicken_ages .== 1], normalize = :probability)
+histogram(egg_samples[chicken_ages .== 1], normalize = :probability)
 
 # ╔═╡ 9d9b4091-cf44-46e6-b81b-74663c9c8dfe
-histogram(sp_egg[chicken_ages .== 3], normalize = :probability)
+histogram(egg_samples[chicken_ages .== 5], normalize = :probability)
 
 # ╔═╡ ff06c070-50a2-43d0-9729-1c47e728ff52
-md"## 5: Birthdays"
+md"## 4: Birthdays"
 
 # ╔═╡ 6ac2238a-16fd-4a8d-b779-8627d87367ed
 md"""
@@ -331,7 +272,7 @@ Sometimes, people are born on the same day of the year.
 # ╔═╡ 01648616-bf50-4f66-82fc-eaae3de22a38
 md"""
 !!! question
-	What is the probability that, in a class of 150 students, 3 or more share a birthday?
+	What is the probability that, in a class of 150 students, 3 or more share a birthday? Assume the probability for a person to be born is equal on every day of the year.
 """
 
 # ╔═╡ 29438eef-a725-4873-9372-194f2f899f12
@@ -366,10 +307,10 @@ end
 bday_model = birthdays(150)
 
 # ╔═╡ b0443045-74a1-4b48-8f6c-a1b5e15eace4
-sp_maxoccs = [bday_model() for i in 1:2000]
+bday_samples = [bday_model() for i in 1:2000]
 
 # ╔═╡ cac00880-15fa-483a-a09f-9b6d1219b0cf
-mean(sp_maxoccs .>= 3)
+mean(bday_samples .>= 3)
 
 # ╔═╡ 13002efe-15f1-4096-be4f-671432a8991e
 md"""
@@ -414,23 +355,7 @@ E_triplebday = p_triplebday * 365 # The expected value of the amount of triple b
 # ╠═597d97fb-cfbe-4cff-bfbb-57df9a2f42aa
 # ╠═6ca31280-a997-4e84-98bb-e7784d0c555c
 # ╠═9dcb122b-dd12-4e80-b391-f780e637d6fe
-# ╟─ce7d57ed-4f31-4dcf-af3b-b37a2e2a9393
-# ╟─087994ce-a26d-40c4-87eb-ef9f0ce7f1fb
-# ╟─9e5cc347-c74f-46a3-9534-c5ad812844bf
-# ╠═47a43282-3892-4a9a-94b7-c359fa74e12b
-# ╠═890eaeb5-c63f-470c-85eb-d77cec68b740
-# ╠═0d0297a3-60b6-47e5-b42b-bde99c09a16e
-# ╠═c3cec28a-f472-46f5-9082-ddfdeb01ddf0
-# ╠═1911f031-d668-47df-995a-ca339ad4ae47
-# ╠═5c13b506-c05b-4258-95e7-93b648defa17
-# ╠═d8ae7b0a-afee-41c0-af0c-f8b1d30086f9
-# ╟─faa105b6-0700-4f4d-92fe-2bb72a4d6e44
-# ╠═b4f156df-0b9c-495a-b95b-000c7f166243
-# ╠═f7eb6faf-8884-4fea-9974-330d7fd555aa
-# ╟─4decd959-aeb9-47d4-a381-14bbf4dbc5ab
-# ╟─a465722f-49f9-4d37-9a2b-00d1120f5c06
-# ╠═b2999f6e-2dbd-44ab-a303-97234f665d33
-# ╠═5ea80b00-941a-4226-87cb-66da7ee7d76d
+# ╟─c976dd77-6a72-4acc-9b07-342cc01fc7ee
 # ╟─9740ea64-cd4f-46b1-a741-02e392280601
 # ╟─187854bb-9e30-454d-9e03-cccf77aebb6b
 # ╟─747e3c0a-357a-448a-b479-d0fcbe44a6c0
diff --git a/exercises/solved_notebooks/P4_probmod/probmod_3-advanced.jl b/exercises/solved_notebooks/P4_probmod/probmod_3-advanced.jl
index 874e41d..7ac3e93 100644
--- a/exercises/solved_notebooks/P4_probmod/probmod_3-advanced.jl
+++ b/exercises/solved_notebooks/P4_probmod/probmod_3-advanced.jl
@@ -1,5 +1,5 @@
 ### A Pluto.jl notebook ###
-# v0.20.4
+# v0.20.21
 
 using Markdown
 using InteractiveUtils
@@ -7,7 +7,7 @@ using InteractiveUtils
 # This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
 macro bind(def, element)
     #! format: off
-    quote
+    return quote
         local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
         local el = $(esc(element))
         global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
@@ -52,7 +52,7 @@ Bacteria follow **logistic growth**, and you can use the following assumptions:
 md"""
 !!! questions
 	1. Plot the prior distribution of P0.
-	2. What is the probability your bacteria are in a splittable state 8 hours after inoculation?
+	2. What is the probability that your bacteria are in a splittable state 8 hours after inoculation?
 	3. Plot 100 of the sampled logistic growth curves from 0 to 12 hours.
 """
 
@@ -84,41 +84,31 @@ rand(dropletdist, 10000) |> histogram
 # ╔═╡ 107533fc-b300-4b8d-bea2-a3aa6a37938d
 md"### 2: Probability"
 
-# ╔═╡ 38c5b9cd-0baf-4a62-912e-4993614ddbf3
-md"""
-!!! tip
-	You can `return` the logistic function estimated within the model and retrieve it using `generated_quantities` to make plotting easier later on.
-
-	Remember: anonymous functions can be defined using `myfun = x -> ...`
-"""
-
 # ╔═╡ bde57599-1dca-41a4-94aa-498da72c2012
 logistic(t, P0, r, K) =  K / (1 + (K - P0)/P0 * exp(-r*t))
 
 # ╔═╡ 976cfb96-3196-4a61-bd5f-e4f5d24ba1e9
-@model function petrigrowth()
-	P0 ~ dropletdist
+@model function petrigrowth(t)
+	P0 ~ MixtureModel([Poisson(10), Poisson(30)], [0.75, 0.25])
     r ~ LogNormal(0.0, 0.3)
 	K ~ Normal(1e5, 1e4)
 
-	logfun = t -> logistic(t, P0, r, K)
-    return logfun
+	num_bacteria = logistic(t, P0, r, K)
+	splittable = 1e4 < num_bacteria < 1e5
+    return splittable
 end
 
 # ╔═╡ 4345b3dd-0731-4dd5-a319-627b4f91306e
-petri_model = petrigrowth();
+petri_model = petrigrowth(8);
 
 # ╔═╡ 9d747eef-a883-49b3-acb7-f0d077a2b902
 chain_petri = sample(petri_model, Prior(), 2000);
 
 # ╔═╡ ce8b53c3-0af1-4e9a-ad80-6fd7a2bf020b
-logfuns = generated_quantities(petri_model, chain_petri);
-
-# ╔═╡ e1d0c5d5-2b7a-4af4-a7ac-e1ca46e665c8
-sp_petri = [logfun(8.0) for logfun in logfuns];
+splittable_samples = generated_quantities(petri_model, chain_petri);
 
 # ╔═╡ e3092915-a083-425e-8fa1-b7bf370abc8a
-prob_splittable = mean((sp_petri .>= 1e4) .&& (sp_petri .<= 1e5))
+prob_splittable = mean(splittable_samples)
 
 # ╔═╡ e4f42f16-5cce-4fc2-aa01-8971f37c710e
 md"### 3: Plot"
@@ -126,11 +116,17 @@ md"### 3: Plot"
 # ╔═╡ 6253f255-e92d-4b5a-9d89-d4fe384fc438
 md"""
 !!! tip
-	Plotting a function `myfun` is as simple as entering `plot(myfun)`. The same syntax applies if `myfun` is a vector of functions. However, don't forget it was asked to plot only **100** functions.
+	Remember: anonymous functions can be defined using `myfun = x -> ...`, and can be visualized using `plot(myfun)`. The same syntax applies if `myfun` is a vector of functions. However, don't forget it was asked to plot only **100** functions.
 """
 
+# ╔═╡ 551e2c0c-faaf-476b-9021-7a10a75e92cf
+logfuns = [
+	t -> logistic(t, chain_petri[:P0][i], chain_petri[:r][i], chain_petri[:K][i])
+	for i in 1:length(chain_petri)
+];
+
 # ╔═╡ 14a8bc88-e10e-42d1-a4c0-7f3ebc5512a9
-plot(logfuns[1:100], color = :violet, alpha = 0.3, label = false, xlims = (0, 12))
+plot(logfuns[1:100], color = :violet, alpha = 0.5, label = false, xlims = (0, 12))
 	# extra arguments are for making the plot prettier but are not required
 
 # ╔═╡ c8941726-9e81-47fc-9b7e-cb3b5c0c61ca
@@ -138,7 +134,7 @@ md"# 2: Attraction"
 
 # ╔═╡ 431023df-3724-4325-b0ac-96dbf5e4fd20
 md"""
-Following a course on electromagnetism will teach one that computing the net force between 2 arbitrary shapes can be a terrifying task. Tragedy has it then, that this is a very general problem with application from making fusion reactors to space travel. We can ease the pain by turning it into a sampling problem.
+Following a course on electromagnetism will teach one that computing the net force between 2 arbitrary shapes can be a terrifying task. Tragedy has it then, that this is a very general problem with applications from making fusion reactors to space travel. We can ease the pain by turning it into a sampling problem.
 
 We'll start in a humble manner and simulate **the gravitational force between 2 cubes**. Both cubes are size 1. The first cube is in [0, 1] x [0, 1] x [0, 1], and the second cube in [1.1, 2.1] x [0, 1] x [0, 1], as shown in the figure below.
 """
@@ -149,7 +145,7 @@ begin
 	ye = [0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1]
 	ze = [0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1]
 
-	xe2 = [0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0] .+ 1.1
+	xe2 = xe .+ 1.1
 
 	plot(xlims = (-0.5, 2.5), ylims = (-1, 2), zlims = (0, 3))
 	plot!(xe, ye, ze; color = :blue, linewidth = 0.5, label = "cube 1")
@@ -158,18 +154,19 @@ end
 
 # ╔═╡ bb8e4cd6-6106-4d3f-9e35-0dfd8b2c45f5
 md"""
-The gravitational force can be estimated by **randomly sampling a point from both cubes** and using the formula for gravitational force between those points, ignoring all constants:
+You can sample the gravitational force between both cubes by **randomly sampling a point from both cubes** and using the formula for gravitational force between those points, ignoring all constants:
 
 ```math
-F = \frac{1}{r^2}
+F = \frac{1}{r^2} \, .
 ```
 """
 
 # ╔═╡ 058d2a08-7d40-46c5-8c5b-2b5ce9d2eb18
 md"""
 !!! questions
-	1. What is the estimated net force between the two cubes? Is this the same as if you had treated the cubes as point masses?
-	1. How many samples do you need to estimate this force reliably? Define a reliable estimator as one having a standard deviation of 0.1. Visualise the distribution of the estimator.
+	1. What is the estimated total force between the two cubes? Is this the same as if you had treated the cubes as point masses?
+	1. To estimate the net force, you take the average of $n$ samples. Of course, the result will vary every time you take a new sample: if you take the average of only 10 samples, your estimated total force will vary wildly! We can quantify how much the estimated total force varies by taking a **sample of sample averages** and then calculating the variance. Assuming you want your sample average to have a variance of no more than $0.01$, how many samples do you need? Visualise the distribution of the sample average.
+	    - EXTRA: How does the [central limit theorem](https://en.wikipedia.org/wiki/Central_limit_theorem) apply to this question? Can you use it to answer the question with less trial and error?
 """
 
 # ╔═╡ 3cd12379-5e7f-4f5d-a7bb-8b3a9d2e8497
@@ -199,7 +196,7 @@ cubemodel = cubeforce();
 force_sp = [cubemodel() for _ in 1:2000];
 
 # ╔═╡ 1ce8f808-b917-4831-89b4-4c318f05724d
-force_average = mean(force_sp)
+estimated_force = mean(force_sp)
 
 # ╔═╡ 420120c3-1c8c-46e4-a7c0-0f2405dd159a
 pointmass_force = 1 / 1.1^2
@@ -211,20 +208,58 @@ histogram(force_sp)
 md"### 2: Variance of Estimator"
 
 # ╔═╡ e38e8f51-90db-4136-84e2-f06cd03d502a
-@bind required_samples Slider(10:10:200, default = 140, show_value = true)
-	# manually change until standard deviation is about 0.1
+@bind n Slider(10:10:200, default = 20, show_value = true)
+	# manually change until variance is about 0.01
 
 # ╔═╡ 788f60fc-b1a1-4635-a36a-f954738bcffc
-estimator = mean([cubemodel() for _ in 1:required_samples])
+mean([cubemodel() for _ in 1:n]) # one sample of the estimated force
 
 # ╔═╡ a7e975a9-beb5-4169-9d22-3e970be2838b
-estimator_sp = [mean([cubemodel() for _ in 1:required_samples]) for _ in 1:2000];
+force_samples = [mean([cubemodel() for _ in 1:n]) for _ in 1:2000]; # "n" samples of the estimated force
 
 # ╔═╡ 7ff5cb25-47bd-432f-a404-359abaf44e3a
-estimator_sp_σ = std(estimator_sp)
+force_samples_var = var(force_samples)
 
 # ╔═╡ 99c7f6d6-eada-491b-b966-fdb1195fc111
-histogram(estimator_sp)
+histogram(force_samples)
+
+# ╔═╡ 5c8ba3b4-fee0-49bd-a5e4-4bef400807e0
+md"#### Extra: The central limit theorem"
+
+# ╔═╡ d695bdfb-1557-4923-9638-1f57689085ca
+md"""
+According to the central limit theorem, for samples $X_i$ following some distribution with mean $μ$ and variance $σ^2$, the distribution of the sample average $\bar{X}_n$ will converge to a normal distribution with mean $μ$ and variance $σ^2/n$.
+
+If we know the variance $σ_1^2$ for some sample size $n_1$, we can then find the required amount of samples $n_r$ for $σ_r^2 = 0.01$ using the following simple derivation:
+
+Given that
+```math
+\begin{align}
+σ^2_1 \approx σ^2 / n_1
+\\ σ^2_r \approx σ^2 / n_r
+\end{align}
+```
+it follows
+```math
+\begin{align}
+σ^2 &= σ^2
+\\ \Rightarrow n_1 σ^2_1 &\approx n_r σ^2_r
+\\ \Rightarrow n_r &\approx \frac{σ_1^2}{σ_r^2} n_1
+\end{align}
+```
+"""
+
+# ╔═╡ 9799aead-e2f5-418a-9079-34f985286fcb
+n_1 = 20
+
+# ╔═╡ f5769c0e-b452-49b6-9553-324eeeaad229
+var_1 = var([mean([cubemodel() for _ in 1:n_1]) for _ in 1:2000])
+
+# ╔═╡ 1a916042-ab26-4e52-9462-6880ca18072c
+var_r = 0.01
+
+# ╔═╡ 07ec86cc-84f1-497e-8cb8-5eaf3a004221
+n_r = var_1 / var_r * n_1 # we need at least ~ 130 samples
 
 # ╔═╡ Cell order:
 # ╟─52a38b60-178b-4a1d-ac32-e73fafd339f9
@@ -241,16 +276,15 @@ histogram(estimator_sp)
 # ╠═6517d682-b524-42e3-8a13-78ed5c1ce0dc
 # ╠═60eabb68-7d65-4f78-97d9-b1e20c2dbd0a
 # ╟─107533fc-b300-4b8d-bea2-a3aa6a37938d
-# ╟─38c5b9cd-0baf-4a62-912e-4993614ddbf3
 # ╠═bde57599-1dca-41a4-94aa-498da72c2012
 # ╠═976cfb96-3196-4a61-bd5f-e4f5d24ba1e9
 # ╠═4345b3dd-0731-4dd5-a319-627b4f91306e
 # ╠═9d747eef-a883-49b3-acb7-f0d077a2b902
 # ╠═ce8b53c3-0af1-4e9a-ad80-6fd7a2bf020b
-# ╠═e1d0c5d5-2b7a-4af4-a7ac-e1ca46e665c8
 # ╠═e3092915-a083-425e-8fa1-b7bf370abc8a
 # ╟─e4f42f16-5cce-4fc2-aa01-8971f37c710e
 # ╟─6253f255-e92d-4b5a-9d89-d4fe384fc438
+# ╠═551e2c0c-faaf-476b-9021-7a10a75e92cf
 # ╠═14a8bc88-e10e-42d1-a4c0-7f3ebc5512a9
 # ╟─c8941726-9e81-47fc-9b7e-cb3b5c0c61ca
 # ╟─431023df-3724-4325-b0ac-96dbf5e4fd20
@@ -268,5 +302,11 @@ histogram(estimator_sp)
 # ╠═e38e8f51-90db-4136-84e2-f06cd03d502a
 # ╠═788f60fc-b1a1-4635-a36a-f954738bcffc
 # ╠═a7e975a9-beb5-4169-9d22-3e970be2838b
-# ╠═7ff5cb25-47bd-432f-a404-359abaf44e3a
 # ╠═99c7f6d6-eada-491b-b966-fdb1195fc111
+# ╠═7ff5cb25-47bd-432f-a404-359abaf44e3a
+# ╟─5c8ba3b4-fee0-49bd-a5e4-4bef400807e0
+# ╟─d695bdfb-1557-4923-9638-1f57689085ca
+# ╠═9799aead-e2f5-418a-9079-34f985286fcb
+# ╠═f5769c0e-b452-49b6-9553-324eeeaad229
+# ╠═1a916042-ab26-4e52-9462-6880ca18072c
+# ╠═07ec86cc-84f1-497e-8cb8-5eaf3a004221
diff --git a/exercises/solved_notebooks/P4_probmod/probmod_4-review.jl b/exercises/solved_notebooks/P4_probmod/probmod_4-review.jl
index e23730f..849a80a 100644
--- a/exercises/solved_notebooks/P4_probmod/probmod_4-review.jl
+++ b/exercises/solved_notebooks/P4_probmod/probmod_4-review.jl
@@ -1,5 +1,5 @@
 ### A Pluto.jl notebook ###
-# v0.20.4
+# v0.20.21
 
 using Markdown
 using InteractiveUtils
@@ -15,7 +15,7 @@ md"# Review exercise: Buffon's needles"
 
 # ╔═╡ a6435b94-f1af-4609-abfc-93d88730d023
 md"""
-A wise man once said: ["there is no greater joy than estimating π"](https://en.wikipedia.org/wiki/Approximations_of_%CF%80). Next to throwing darts at the unit square, another method to accomplish this is using [Buffon's needle problem](https://en.wikipedia.org/wiki/Buffon%27s_needle_problem).
+A wise man once said: ["there is no greater joy than estimating π"](https://en.wikipedia.org/wiki/Approximations_of_%CF%80). One method to accomplish this is using [Buffon's needle problem](https://en.wikipedia.org/wiki/Buffon%27s_needle_problem).
 
 The experiment is as follows: consider a floor with parallel lines all a distance of 1 away from eachother. Now drop a needle of length 1 (and width ~0) on the floor with a **random position and angle**. What is the probability $P_{cross}$ that the needle will cross one of the lines?
 """
@@ -24,8 +24,8 @@ The experiment is as follows: consider a floor with parallel lines all a distanc
 md"The following image illustrates the problem (imagine $l$ = $t$ = 1) for two needles, where `a` crosses a line and `b` does not."
 
 # ╔═╡ c2ba2183-957f-4b94-a22e-0c2f3dd957ad
-md"""
-![Buffon's needles](https://upload.wikimedia.org/wikipedia/commons/thumb/5/58/Buffon_needle.svg/1920px-Buffon_needle.svg.png)
+html"""
+<img src="https://upload.wikimedia.org/wikipedia/commons/thumb/5/58/Buffon_needle.svg/1920px-Buffon_needle.svg.png" style="background-color:white;" alt:"Buffon's needles">
 """
 
 # ╔═╡ d71cf8dc-99e5-48e0-9abe-2242a6ccc30b
@@ -66,13 +66,13 @@ end
 needle_model = needles()
 
 # ╔═╡ 49a206fb-6390-45fa-900f-8031b77a2725
-sp_n = [needle_model() for i in 1:2000]
+needle_samples = [needle_model() for i in 1:2000]
 
 # ╔═╡ b6009fc4-2f67-4b56-acc9-ccc1a013f5e9
-mean(sp_n)
+mean(needle_samples)
 
 # ╔═╡ 63e0b367-6155-472d-8502-b405db10e979
-pi_est = 2 / mean(sp_n)
+pi_est = 2 / mean(needle_samples)
 
 # ╔═╡ Cell order:
 # ╠═882b942c-a5ae-445d-a946-6eb8cddc423e
diff --git a/exercises/solved_notebooks/P5_mcmc/MCMC_1-intro.jl b/exercises/solved_notebooks/P5_mcmc/MCMC_1-intro.jl
index 7c7f57c..c73f3df 100644
--- a/exercises/solved_notebooks/P5_mcmc/MCMC_1-intro.jl
+++ b/exercises/solved_notebooks/P5_mcmc/MCMC_1-intro.jl
@@ -281,10 +281,10 @@ Taking the sampled values of the mutation rate from the chain and plotting a his
 """
 
 # ╔═╡ b657217a-6ccc-4a41-b852-df4e39a7a10a
-sp_alpha = mutation_chain[:α];
+alpha_samples = mutation_chain[:α];
 
 # ╔═╡ 9cf08616-d599-42a3-82e8-e8c98853c1d8
-histogram(sp_alpha) # note the difference with the prior distribution!
+histogram(alpha_samples) # note the difference with the prior distribution!
 
 # ╔═╡ e02c42dc-627c-4bc0-8ebb-fdb2b0f15b64
 md"Plotting some sampled mutation rates from this distribution onto our data shows that they fit well:"
@@ -295,17 +295,17 @@ begin
 		ylabel = "Number of mutations", label = false, xlims = (0, 500), 
 		ylims = (0, 400), title = "A posteriori relationship between t and mean N"
 	);
-	plot!([x -> αᵢ*x for αᵢ in sp_alpha[1:10:end]], color = :purple, opacity = 0.1, label = false)
+	plot!([x -> αᵢ*x for αᵢ in alpha_samples[1:10:end]], color = :purple, opacity = 0.1, label = false)
 end
 
 # ╔═╡ 87059440-2919-4f6d-9d32-0df3ce75e2a2
 md"And finally we can answer our first question:"
 
 # ╔═╡ 644cff58-68b9-4d4b-8896-617fcacc39c5
-mean(sp_alpha)
+mean(alpha_samples)
 
 # ╔═╡ 5f2f5a78-ca90-4baf-8bf0-ff1fae88d785
-sqrt(var(sp_alpha))
+sqrt(var(alpha_samples))
 
 # ╔═╡ bfbf811f-7b2f-473b-b132-0c7e965c8b0d
 md"""
diff --git a/exercises/solved_notebooks/P5_mcmc/MCMC_2-basics.jl b/exercises/solved_notebooks/P5_mcmc/MCMC_2-basics.jl
index a1bf6fa..8eee3ba 100644
--- a/exercises/solved_notebooks/P5_mcmc/MCMC_2-basics.jl
+++ b/exercises/solved_notebooks/P5_mcmc/MCMC_2-basics.jl
@@ -71,10 +71,10 @@ end
 molemodel = mole()
 
 # ╔═╡ b5bd5114-0a0c-4aa3-af00-4f3d05edc5e3
-spY = [molemodel() for i in 1:2000];
+Y_samples = [molemodel() for i in 1:2000];
 
 # ╔═╡ 4e8b9000-cb61-4f01-9ba9-17276ad0335e
-E_Y = mean(spY)
+E_Y = mean(Y_samples)
 
 # ╔═╡ 8777b133-7d7c-4a85-b89c-2f00093e9984
 md"### 3: Conditional expected value of X"
@@ -89,7 +89,7 @@ molechain = sample(cond_mole, NUTS(), 2000);
 plot(molechain)
 
 # ╔═╡ b79dafee-4f6f-4301-9267-1b9798c125ea
-spXcondY = molechain[:X];
+X_samplescondY = molechain[:X];
 
 # ╔═╡ fd540765-95e5-4071-81f7-e689b06cad0c
 E_XcondY = mean(molechain[:X])
diff --git a/exercises/solved_notebooks/P5_mcmc/MCMC_3-advanced.jl b/exercises/solved_notebooks/P5_mcmc/MCMC_3-advanced.jl
index 60facfa..82577cd 100644
--- a/exercises/solved_notebooks/P5_mcmc/MCMC_3-advanced.jl
+++ b/exercises/solved_notebooks/P5_mcmc/MCMC_3-advanced.jl
@@ -179,10 +179,10 @@ plot(petrichain)
 logfuns = generated_quantities(petrimodel, petrichain);
 
 # ╔═╡ c642574c-c01b-4398-aac9-43e514d7fa25
-sp_petri = [logfun(8.0) for logfun in logfuns]
+petri_samples = [logfun(8.0) for logfun in logfuns]
 
 # ╔═╡ d15fa091-839c-4239-96e2-eaf7335ce620
-prob_splittable = mean((sp_petri .>= 1e4) .&& (sp_petri .<= 1e5))
+prob_splittable = mean((petri_samples .>= 1e4) .&& (petri_samples .<= 1e5))
 
 # ╔═╡ 0ea95b67-d4da-4c5a-ad2e-05024ad074a3
 md"### 2"
@@ -227,8 +227,8 @@ let # so we dont need to rename all variables
 	petrimodel = petrigrowth🌟(5) | (P_obs = 21_000,)
 	petrichain = sample(petrimodel, NUTS(), 2_000)
 	logfuns = generated_quantities(petrimodel, petrichain);
-	sp_petri = [logfun(8.0) for logfun in logfuns]
-	prob_splittable = mean((sp_petri .>= 1e4) .&& (sp_petri .<= 1e5))
+	petri_samples = [logfun(8.0) for logfun in logfuns]
+	prob_splittable = mean((petri_samples .>= 1e4) .&& (petri_samples .<= 1e5))
 	println("The new `prob_splittable` is ", prob_splittable)
 	plot(petrichain)
 end
diff --git a/exercises/student_notebooks/P4_probmod/probmod_1-intro.jl b/exercises/student_notebooks/P4_probmod/probmod_1-intro.jl
index 5f0b0b0..411d1f4 100644
--- a/exercises/student_notebooks/P4_probmod/probmod_1-intro.jl
+++ b/exercises/student_notebooks/P4_probmod/probmod_1-intro.jl
@@ -78,9 +78,9 @@ x(t) &= \mathrm{cos}(α) \, v_0 \, t \, ,
 ```
 
 where the parameters are:
--  $v_0$ is the initial velocity.
--  $\alpha$ is the launching angle.
--  $g = 9.8$ is the gravitational acceleration.
+-  $v_0$: the initial velocity.
+-  $\alpha$: the launching angle.
+-  $g = 9.8$: the gravitational acceleration.
 """
 
 # ╔═╡ ee9cd05c-7fa5-4412-b3b4-32b0a1774c87
@@ -91,11 +91,11 @@ t_f = \frac{x_f}{\mathrm{cos}(\alpha \, v_0)} \, ,
 ```
 evaluating it's y-coordinate there
 ```math
-y_d = \mathrm{sin}(α) \, v_0 \, t_f - \frac{g}{2} \, t_f^2 \, ,
+y_f = \mathrm{sin}(α) \, v_0 \, t_f - \frac{g}{2} \, t_f^2 \, ,
 ```
 and checking if it's in the range of our face
 ```math
-y_d ∈ [1.5, 1.7] \, .
+y_f ∈ [1.5, 1.7] \, .
 ```
 """
 
@@ -143,7 +143,13 @@ plot(
 )
 
 # ╔═╡ 314ba0ec-015f-488d-ac63-780b1dbf5bc2
-md"The latter is an **input variable**: when cycling, we can choose how far away we stay from the cyclist in front of us. The x-position of our face $x_f$ is therefore an input to the model."
+md"The third is an **input variable**: when cycling, we can choose how far away we stay from the cyclist in front of us. The x-position of our face $x_f$ is therefore an input to the model."
+
+# ╔═╡ 2878b32a-9e54-40b0-b151-031b95ad161b
+md"""
+!!! note
+	As you can see, to visualize a distribution you can simply call the `plot` function on the distribution (and set extra keyword arguments as desired).
+"""
 
 # ╔═╡ c070f3c9-dd80-49f9-b279-c77f46316116
 md"## Defining the model"
@@ -182,11 +188,11 @@ Our model can be defined as follows:
 	# get the time it takes the droplet to travel the distance x_f
 	t_f = x_f / (cos(α) * v0)
 	# get the droplet's altitude at this time
-	y_s = sin(α)*v0*t_f - g/2 * t_f^2
+	y_f = sin(α)*v0*t_f - g/2 * t_f^2
 	# to respect the ground, set altitude to 0 if it would be negative
-	y_s = max(0, y_s)
+	y_f = max(0, y_f)
 	
-	return y_s 
+	return y_f
 end
 
 # ╔═╡ 22144610-0c37-4501-b14c-83895d8b0eae
@@ -235,22 +241,22 @@ md"""
 md"### Getting the model output"
 
 # ╔═╡ 46f17c07-9aec-4fee-8590-4eeb5a13b93e
-md"To answer our question, we are interested in the `y_s` variable returned by the model. We can extract the output variable from our chain using the `generated_quantities` function:"
+md"To answer our question, we are interested in the `y_f` variable returned by the model. We can extract the output variable from our chain using the `generated_quantities` function:"
 
 # ╔═╡ 071fc405-7a33-4dab-a69a-1e68a0fb35ad
-y_s_samples = generated_quantities(splash_model, chain)
+y_f_samples = generated_quantities(splash_model, chain)
 
 # ╔═╡ 20e99145-e80b-4259-9bca-f7e14c2d355e
 md"It's often useful to visualize the distribution of our variable of interest. We can use the `histogram` function for this:"
 
 # ╔═╡ 45d44cf5-2289-4947-8d4e-1dcddb5447de
-histogram(y_s_samples, bins = 15)
+histogram(y_f_samples, bins = 15)
 
 # ╔═╡ 8d179748-40aa-4e88-b277-9234fe938251
 md"Checking whether the droplet altitude matches our face's can be done with a simple inequality operation:"
 
 # ╔═╡ 633a2b4c-379e-4e18-bbcd-b450930196a5
-face_hit_samples = [1.5 <= y_s <= 1.7 for y_s in y_s_samples]
+face_hit_samples = [1.5 <= y_f <= 1.7 for y_f in y_f_samples]
 
 # ╔═╡ ad8a6ea2-30ba-4104-824f-12a7ea8e9718
 md"Finally, we can estimate the probability of our face being hit by taking the mean amount of times we got hit in our sample."
@@ -264,7 +270,7 @@ md"""And that's the answer to our question! If we're 1m behind the cyclist in fr
 # ╔═╡ bf5642f4-09d3-491d-8eb9-47fcdb67c977
 md"""
 !!! note
-	Do you like your code sleek? You can also calculate this probability from your samples with just one line: `mean(1.5 .<= y_s_samples .<= 1.7)`
+	Do you like your code sleek? You can also calculate this probability from your samples with just one line: `mean(1.5 .<= y_f_samples .<= 1.7)`
 """
 
 # ╔═╡ 3067af6e-df69-48a0-a024-a6e219dcfc7c
@@ -335,6 +341,49 @@ md"""
 	We set a lot of keyword arguments of our plot to make the figure pretty here, but don't panic: this is still no programming course, so we will always clearly provide which one(s) you need if you need to plot something. If you ever do want to look for one yourself, you can find them on the function's help page (`?plot` or the `🔍 Live docs` in the bottom right corner).
 """
 
+# ╔═╡ 1de58e6a-4ee2-4d67-b97d-e63bbb9ed057
+md"## Sampling alternative"
+
+# ╔═╡ 83a8b9b6-72da-4aff-8364-37846b527a2f
+md"It's also possible to skip the chain construction entirely and simply call your model as a function to get a sample of the model output. Note that it's not possible to get the corresponding values of the stochastic variables using this approach, so only use this method if you **only** care about the model output."
+
+# ╔═╡ 771b9cb1-d809-46de-9380-2db73ccd014d
+# get a sample of the output 500 times
+y_f_samples_alternative = [splash_model() for _ in 1:500]
+
+# ╔═╡ 0f1cccd8-937d-4fc2-8657-c9967db483b8
+histogram(y_f_samples_alternative, bins = 15) # yup, that's the same distribution
+
+# ╔═╡ 7806eca0-a2ad-4910-8d6e-edc50073b37b
+md"## For-loops for many variables"
+
+# ╔═╡ 64ebc3ce-83ca-4cee-ab93-cddf51138543
+md"""
+Sometimes you need to define a large number of stochastic variables following the same distribution. Let's say, for example, you want to model the total mass of water being splashed on your face assuming that 10 droplets will hit you, each having a mass that is exponentially distributed with a mean of 30 mg. 
+
+This can be done using a **for-loop** as follows:
+"""
+
+# ╔═╡ bdf3c18c-b40b-4fb4-81e7-fe9d821cc166
+@model function splash_mass()
+ 	droplet_masses = zeros(10) # instantiate vector of 10 zeros
+	for i in 1:length(droplet_masses)
+		droplet_masses[i] ~ Exponential(30) # each element of the vector follows an Exponential distribution
+	end
+	total_mass = sum(droplet_masses)
+	
+	return total_mass
+end
+
+# ╔═╡ 2cbe2750-38c1-4e7f-8d03-9bed7c16861a
+mass_model = splash_mass();
+
+# ╔═╡ ad8b48c5-f892-49b3-8f19-64d93424485c
+mass_chain = sample(mass_model, Prior(), 1000);
+
+# ╔═╡ 8507f084-eff5-4351-800d-8eba7ca673bb
+total_mass_samples = generated_quantities(mass_model, mass_chain)
+
 # ╔═╡ befb8f05-6ab7-4594-9cf1-71df47c1df03
 md"## Essentials"
 
@@ -358,33 +407,32 @@ md"""
 # ╠═╡ show_logs = false
 let
 	# model definition (comment out for speed)
-	g = 9.8
 	@model function splash(x_f)
 		v0 ~ Normal(5, 1) 
 		α ~ TriangularDist(0, pi/2)
 	
 		t_f = x_f / (cos(α) * v0)
-		y_s = sin(α)*v0*t_f - g/2 * t_f^2
-		y_s = max(0, y_s)
-		return y_s 
+		y_f = sin(α)*v0*t_f - g/2 * t_f^2
+		y_f = max(0, y_f)
+		return y_f
 	end
 
 	# get samples
 	splash_model = splash(x_f)
 	chain = sample(splash_model, Prior(), 500);
 
-	y_s_samples = generated_quantities(splash_model, chain)
+	y_f_samples = generated_quantities(splash_model, chain)
 	v0_samples = chain["v0"]
 	α_samples = chain[:α]
 
 	# calculate interesting things
-	prob_face_hit = mean([1.5 <= y_s <= 1.7 for y_s in y_s_samples])
+	prob_face_hit = mean([1.5 <= y_f <= 1.7 for y_f in y_f_samples])
+
+	y(x, α, v0) = sin(α)/cos(α)*x - g/2*(x / (cos(α) * v0))^2
 	trajectories = [
-		x -> sin(α_samples[i]) / cos(α_samples[i]) * x -
-		g / 2 * ( x / (cos(α_samples[i])*v0_samples[i]) )^2
+		x -> y(x, α_samples[i], v0_samples[i])
 		for i in 1:length(v0_samples)
 	]
-
 	plot(trajectories, xlims = (0, 5), ylims = (0, 2), 
 		 legend = false, color = :blue, alpha = 0.3,
 		 title = "Probability of getting hit ≈ $(round(100*prob_face_hit, digits = 3))%"
@@ -411,10 +459,11 @@ end
 # ╟─e8821ac2-baa3-4674-868c-3a5b8593cc95
 # ╟─4811e1af-9ba4-49b3-a3a1-5d8ff8d5af92
 # ╟─01156ef7-b194-4372-8bc8-de9fad5e12f1
-# ╟─f30290bd-9d31-4939-b5d6-f89f066775d2
+# ╠═f30290bd-9d31-4939-b5d6-f89f066775d2
 # ╟─dab23eba-3eb1-434b-aef5-e08730e35114
-# ╟─299e402c-634a-4b18-b867-489bafaf4564
+# ╠═299e402c-634a-4b18-b867-489bafaf4564
 # ╟─314ba0ec-015f-488d-ac63-780b1dbf5bc2
+# ╟─2878b32a-9e54-40b0-b151-031b95ad161b
 # ╟─c070f3c9-dd80-49f9-b279-c77f46316116
 # ╟─84a08354-7bf3-4783-9285-a7e8850e798b
 # ╟─8f94567c-b87b-4eb4-8441-c6434da393a8
@@ -442,7 +491,7 @@ end
 # ╟─ad8a6ea2-30ba-4104-824f-12a7ea8e9718
 # ╠═e1f9d1e9-8a0c-4b28-ab67-a89e3287b50f
 # ╟─1c9e2de1-21a0-4d4d-918d-9e2f0fc4d912
-# ╟─bf5642f4-09d3-491d-8eb9-47fcdb67c977
+# ╠═bf5642f4-09d3-491d-8eb9-47fcdb67c977
 # ╟─3067af6e-df69-48a0-a024-a6e219dcfc7c
 # ╟─cbc04c5f-aa93-4340-a61e-1583d20d86b1
 # ╟─d36784b8-9f3b-4434-88fe-03b75b94754e
@@ -457,6 +506,16 @@ end
 # ╟─00572e1a-d1e6-4c84-b1b2-b53404087dd8
 # ╠═f550ed7a-1bba-40b5-9338-cdeb0520aae9
 # ╟─3b465fc0-3a9d-4e67-8d4a-a14b3cd61f4e
+# ╟─1de58e6a-4ee2-4d67-b97d-e63bbb9ed057
+# ╟─83a8b9b6-72da-4aff-8364-37846b527a2f
+# ╠═771b9cb1-d809-46de-9380-2db73ccd014d
+# ╠═0f1cccd8-937d-4fc2-8657-c9967db483b8
+# ╟─7806eca0-a2ad-4910-8d6e-edc50073b37b
+# ╟─64ebc3ce-83ca-4cee-ab93-cddf51138543
+# ╠═bdf3c18c-b40b-4fb4-81e7-fe9d821cc166
+# ╠═2cbe2750-38c1-4e7f-8d03-9bed7c16861a
+# ╠═ad8b48c5-f892-49b3-8f19-64d93424485c
+# ╠═8507f084-eff5-4351-800d-8eba7ca673bb
 # ╟─befb8f05-6ab7-4594-9cf1-71df47c1df03
 # ╟─52af5de0-b516-4947-84f7-dde63e4c513e
 # ╟─2ceede7e-8f0c-420c-bde2-f16b726cd204
diff --git a/exercises/student_notebooks/P4_probmod/probmod_2-basics.jl b/exercises/student_notebooks/P4_probmod/probmod_2-basics.jl
index 27f2d54..8a31b13 100644
--- a/exercises/student_notebooks/P4_probmod/probmod_2-basics.jl
+++ b/exercises/student_notebooks/P4_probmod/probmod_2-basics.jl
@@ -1,5 +1,5 @@
 ### A Pluto.jl notebook ###
-# v0.20.4
+# v0.20.21
 
 using Markdown
 using InteractiveUtils
@@ -45,11 +45,17 @@ end
 # ╔═╡ 42c18a70-efb3-436b-83bb-b586280d4a4e
 dpmodel = doublepoisson();
 
+# ╔═╡ a21b2f44-c8d7-4838-926f-8fa39c36884e
+dpchain = missing
+
 # ╔═╡ b22a18d5-f70b-42a5-a09e-3de515148a6d
-spY = missing
+Y_samples = missing
+
+# ╔═╡ a11629df-31d2-4ac7-bf8e-80b910cab2fb
+missing # plot of X
 
 # ╔═╡ 34bd60df-272f-49d7-9346-fb4d125fe89b
-histogram(spY)
+missing # histogram of Y
 
 # ╔═╡ 13989ba4-bcf8-4fdd-8aee-ab58c8905bc9
 md"### 2: Probabilities"
@@ -59,9 +65,9 @@ md"""
 !!! tip
 	When comparing a vector of values to a single number, don't forget to use `.` to execute operations element-wise in Julia!
 	
-	✅ `spY .< 1` compares every element of `spY` to `1`
+	✅ `Y_samples .< 1` compares every element of `Y_samples` to `1`
 	
-	❌ `spY < 1` compares an entire vector with a single number → errors :(
+	❌ `Y_samples < 1` compares an entire vector with a single number → errors :(
 """
 
 # ╔═╡ 39e4f3eb-b7a9-4ece-846f-eb02dbd77860
@@ -78,11 +84,11 @@ md"""
 !!! hint
 	To create a sample of $X$ that is conditional on some value(s) of $Y$, you can start from a sample of $X$ and select only those elements for which the corresponding sample of $Y$ has the conditioned value(s).
 
-	In other words, you'll need to index `spX` based on `spY` (and vice versa for $\text{var}(Y ∣ X)$). 
+	In other words, you'll need to index `X_samples` based on `Y_samples` (and vice versa for $\text{var}(Y ∣ X)$). 
 """
 
 # ╔═╡ 263048e5-206c-499e-836c-bfe489ed9b74
-spX = missing;
+X_samples = missing;
 
 # ╔═╡ aedd0fe8-da3e-4463-b0cf-7c4f9a22db52
 varXcondY = missing
@@ -93,53 +99,15 @@ varYcondX = missing
 # ╔═╡ ce9e2ce2-26e0-4f17-adf5-c922ba98239d
 missing # analytical answer of (missing)
 
-# ╔═╡ ce7d57ed-4f31-4dcf-af3b-b37a2e2a9393
-md"## 2: Combinations"
-
-# ╔═╡ 087994ce-a26d-40c4-87eb-ef9f0ce7f1fb
-md"""
-Let `U ~ Uniform(0, 4)`, `V ∼ Normal(U, 1)` and `W ~ TriangularDist(0, 4, U)`. 
-1. Use sampling (n = 10_000) to make a histogram of `|V − W|`.
-2. Estimate `P(V > W)` and `P(V * W >= 10)`
-3. Are `V` and `W` independent?
-"""
-
-# ╔═╡ 9e5cc347-c74f-46a3-9534-c5ad812844bf
-md"### 1: Histogram"
-
-# ╔═╡ 47a43282-3892-4a9a-94b7-c359fa74e12b
-@model function combinations()
-	U ~ missing
-	V ~ missing
-	W ~ missing
-end
-
-# ╔═╡ b0502f20-17af-4ecc-be80-a26b3e42d57f
-spVW = missing
-
-# ╔═╡ faa105b6-0700-4f4d-92fe-2bb72a4d6e44
-md"### 2: Probabilities"
-
-# ╔═╡ 97aec4a9-2954-4967-bc15-c0123bac2e75
-probVW1 = missing
-
-# ╔═╡ 30166840-d5f5-4a2a-acc2-cde76a87e95a
-probVW2 = missing
-
-# ╔═╡ 4decd959-aeb9-47d4-a381-14bbf4dbc5ab
-md"### 3: Independence"
-
-# ╔═╡ aa953baa-5105-49d7-82e6-94ca462624f7
+# ╔═╡ 19e16f77-ea34-45f6-83eb-8d256e5fd10d
 md"""
-!!! hint
-	One way to prove dependence is showing that $E[V] \neq E[V \mid W \leq w]$ for at least one value $w$. 
-
-	If the expected value of $V$ can change based on some information about $W$, they can't be independent!
+!!! warning
+	Starting here, only part of the answer's structure will be given. You are therefore expected to **add more code cells** yourself, using the `+` at the left in between two cells.
 """
 
 # ╔═╡ 9740ea64-cd4f-46b1-a741-02e392280601
 md"""
-## 3: Dice
+## 2: Dice
 """
 
 # ╔═╡ 187854bb-9e30-454d-9e03-cccf77aebb6b
@@ -163,7 +131,7 @@ md"### 1: Watercube"
 
 # ╔═╡ 4fefa78b-746d-4e25-baee-d791eb22d930
 md"""
-!!! hint
+!!! tip
 	Consider the humble `DiscreteUniform` distribution. Not sure how it works? Open the **🔍 Live Docs** at the bottom right of the screen for more information
 """
 
@@ -178,11 +146,11 @@ md"""
     return dicesum
 end
 
-# ╔═╡ 02bc37e2-5cf6-404a-b3fe-b2120671adb2
-watermodel = watercube()
+# ╔═╡ 5fe1080e-8c2c-4c6b-84b4-3857dc1f399e
+watercube_samples = missing # samples of dicesum for watercube
 
 # ╔═╡ b5886255-5c1d-4d84-b7ee-6c690fa526dc
-p_watercube_kills = missing
+p_watercube_kills = missing # probability that watercube kills the monster
 
 # ╔═╡ 84e06162-e3f1-4fd7-baf6-5095172413d2
 missing # histogram
@@ -194,15 +162,12 @@ md"### 2: Dirtprism"
 @model function dirtprism()
 	# check the "For-loops for many variables" section from the intro notebook!
 	
-	dicesum = missing # consider the `sum` function
+	dicesum = missing
 	return dicesum
 end
 
-# ╔═╡ 83afb3c1-9e6a-4d18-b0c7-05ed0173df40
-dirtmodel = dirtprism()
-
 # ╔═╡ d9152416-8a7b-480c-ba9f-7ab15404b7a6
-p_dirtprism_kills = missing
+p_dirtprism_kills = missing # probability that dirtprism kills the monster
 
 # ╔═╡ 90e058d7-b3fe-4c42-a652-3c42bf9d851a
 missing # histogram
@@ -214,17 +179,17 @@ md"### 3: Comparison"
 p_watercube_is_better = missing
 
 # ╔═╡ 34f3014f-f4d4-43d1-b46f-bdca73aee33f
-md"## 4: Super eggs"
+md"## 3: Super eggs"
 
 # ╔═╡ 372436c4-262f-49b8-b1cf-626b043542bf
 md"""
-When a chicken lays an egg, there's a small chance it contains two egg yolks. This chance, as well as the number of eggs a chicken lays per year, go down as the chicken gets older.
+When a chicken lays an egg, there's a small chance it contains two egg yolks. This chance, as well as the number of eggs a chicken lays per year, goes down as the chicken gets older.
 """
 
 # ╔═╡ 20111742-008a-44c3-8c27-62791cce3e1e
 md"""
 You can make the following assumptions
-- The age $A$ of a random chicken (in years) is discrete and Uniformly distributed between 0 and 12.
+- The age $A$ of a random chicken (in years) is discrete and Uniformly distributed between 0 and 10.
 - The number of eggs $N$ an $A$-year old chicken lays in a year is Poisson distributed with mean $300 - 20 \, A$.
 - The probability $P$ of an $A$-year old chicken's egg having a double yolk is distributed as a `Beta(1, 800 + 100*A)`.
 """
@@ -233,7 +198,7 @@ You can make the following assumptions
 md"""
 !!! questions
     1. If someone hands you a random chicken, what is the probability it will lay 2 or more double eggs in a year? 
-    1. Compare the distributions of double eggs for 1-year old and 3-year old chickens.
+    1. Compare the distributions of double eggs for 1-year old and 5-year old chickens.
 """
 
 # ╔═╡ 98fcfcef-bf63-4eae-a325-ed4cef6d4fdd
@@ -241,7 +206,7 @@ md"### 1: Probability"
 
 # ╔═╡ a0341046-c14a-494d-a9bd-a60c207c9e76
 md"""
-!!! hint
+!!! tip
 	In this exercise, the output variable (the number of double-yolked eggs) is **also a random variable**! In other words, it also follows some distribution.
 
 	When considering what distribution, consider that each of the $N$ eggs represents a "trial" with a $P$ chance of success for a double yolk.
@@ -265,7 +230,7 @@ missing # histogram 1
 missing # histogram 2
 
 # ╔═╡ ff06c070-50a2-43d0-9729-1c47e728ff52
-md"## 5: Birthdays"
+md"## 4: Birthdays"
 
 # ╔═╡ 6ac2238a-16fd-4a8d-b779-8627d87367ed
 md"""
@@ -275,7 +240,7 @@ Sometimes, people are born on the same day of the year.
 # ╔═╡ 01648616-bf50-4f66-82fc-eaae3de22a38
 md"""
 !!! question
-	What is the probability that, in a class of 150 students, 3 or more share a birthday?
+	What is the probability that, in a class of 150 students, 3 or more share a birthday? Assume the probability for a person to be born is equal on every day of the year.
 """
 
 # ╔═╡ da44d18c-8be3-446e-a5c2-905af545d2c6
@@ -305,7 +270,9 @@ end
 # ╟─def27d26-2205-4a66-94f3-eddbc17483bf
 # ╠═ff38df99-f843-414d-8e45-b46e06a65c22
 # ╠═42c18a70-efb3-436b-83bb-b586280d4a4e
+# ╠═a21b2f44-c8d7-4838-926f-8fa39c36884e
 # ╠═b22a18d5-f70b-42a5-a09e-3de515148a6d
+# ╠═a11629df-31d2-4ac7-bf8e-80b910cab2fb
 # ╠═34bd60df-272f-49d7-9346-fb4d125fe89b
 # ╟─13989ba4-bcf8-4fdd-8aee-ab58c8905bc9
 # ╟─aba42086-224f-44a0-b616-f8f651afdd18
@@ -317,28 +284,18 @@ end
 # ╠═aedd0fe8-da3e-4463-b0cf-7c4f9a22db52
 # ╠═7ef53a87-e5df-4724-b896-3d1d46214c68
 # ╠═ce9e2ce2-26e0-4f17-adf5-c922ba98239d
-# ╟─ce7d57ed-4f31-4dcf-af3b-b37a2e2a9393
-# ╟─087994ce-a26d-40c4-87eb-ef9f0ce7f1fb
-# ╟─9e5cc347-c74f-46a3-9534-c5ad812844bf
-# ╠═47a43282-3892-4a9a-94b7-c359fa74e12b
-# ╠═b0502f20-17af-4ecc-be80-a26b3e42d57f
-# ╟─faa105b6-0700-4f4d-92fe-2bb72a4d6e44
-# ╠═97aec4a9-2954-4967-bc15-c0123bac2e75
-# ╠═30166840-d5f5-4a2a-acc2-cde76a87e95a
-# ╟─4decd959-aeb9-47d4-a381-14bbf4dbc5ab
-# ╟─aa953baa-5105-49d7-82e6-94ca462624f7
+# ╟─19e16f77-ea34-45f6-83eb-8d256e5fd10d
 # ╟─9740ea64-cd4f-46b1-a741-02e392280601
 # ╟─187854bb-9e30-454d-9e03-cccf77aebb6b
 # ╟─747e3c0a-357a-448a-b479-d0fcbe44a6c0
 # ╟─eda95f45-083a-4e65-b57e-bd9890da1f9c
 # ╟─4fefa78b-746d-4e25-baee-d791eb22d930
 # ╠═36477423-5628-4ed3-b54d-9a050557f6b7
-# ╠═02bc37e2-5cf6-404a-b3fe-b2120671adb2
+# ╠═5fe1080e-8c2c-4c6b-84b4-3857dc1f399e
 # ╠═b5886255-5c1d-4d84-b7ee-6c690fa526dc
 # ╠═84e06162-e3f1-4fd7-baf6-5095172413d2
 # ╟─a1b933ac-5d1b-4800-a6e8-e942846b19d8
 # ╠═6b009ba3-83a8-4176-86d4-dd9f70ed29ec
-# ╠═83afb3c1-9e6a-4d18-b0c7-05ed0173df40
 # ╠═d9152416-8a7b-480c-ba9f-7ab15404b7a6
 # ╠═90e058d7-b3fe-4c42-a652-3c42bf9d851a
 # ╟─49790a8f-9f53-4ba7-9543-d6a879b520e0
diff --git a/exercises/student_notebooks/P4_probmod/probmod_3-advanced.jl b/exercises/student_notebooks/P4_probmod/probmod_3-advanced.jl
index 8dac90d..f61e298 100644
--- a/exercises/student_notebooks/P4_probmod/probmod_3-advanced.jl
+++ b/exercises/student_notebooks/P4_probmod/probmod_3-advanced.jl
@@ -1,5 +1,5 @@
 ### A Pluto.jl notebook ###
-# v0.20.4
+# v0.20.21
 
 using Markdown
 using InteractiveUtils
@@ -7,7 +7,7 @@ using InteractiveUtils
 # This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
 macro bind(def, element)
     #! format: off
-    quote
+    return quote
         local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
         local el = $(esc(element))
         global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
@@ -52,7 +52,7 @@ Bacteria follow **logistic growth**, and you can use the following assumptions:
 md"""
 !!! questions
 	1. Plot the prior distribution of P0.
-	2. What is the probability your bacteria are in a splittable state 8 hours after inoculation?
+	2. What is the probability that your bacteria are in a splittable state 8 hours after inoculation?
 	3. Plot 100 of the sampled logistic growth curves from 0 to 12 hours.
 """
 
@@ -72,44 +72,20 @@ md"""
 """
 
 # ╔═╡ eb6dc4e7-e779-4bfc-b865-3defa3894181
-dropletdist = missing;
+missing # plot
 
 # ╔═╡ 107533fc-b300-4b8d-bea2-a3aa6a37938d
 md"### 2: Probability"
 
-# ╔═╡ 38c5b9cd-0baf-4a62-912e-4993614ddbf3
-md"""
-!!! tip
-	You can `return` the logistic function estimated within the model and retrieve it using `generated_quantities` to make plotting easier later on.
-
-	Remember: anonymous functions can be defined using `myfun = x -> ...`
-"""
-
 # ╔═╡ bde57599-1dca-41a4-94aa-498da72c2012
 logistic(t, P0, r, K) =  K / (1 + (K - P0)/P0 * exp(-r*t))
 
 # ╔═╡ 976cfb96-3196-4a61-bd5f-e4f5d24ba1e9
-@model function petrigrowth()
-	P0 ~ dropletdist
-    r ~ missing
-	K ~ missing
-
-	logfun = missing
-    return logfun
+@model function petrigrowth(t)
+	missing
+    return splittable
 end
 
-# ╔═╡ 4345b3dd-0731-4dd5-a319-627b4f91306e
-petri_model = missing
-
-# ╔═╡ 9d747eef-a883-49b3-acb7-f0d077a2b902
-chain_petri = missing
-
-# ╔═╡ ce8b53c3-0af1-4e9a-ad80-6fd7a2bf020b
-logfuns = missing
-
-# ╔═╡ e1d0c5d5-2b7a-4af4-a7ac-e1ca46e665c8
-sp_petri = missing
-
 # ╔═╡ e3092915-a083-425e-8fa1-b7bf370abc8a
 prob_splittable = missing
 
@@ -119,9 +95,12 @@ md"### 3: Plot"
 # ╔═╡ 2b39e0bd-91cd-456e-9053-7b8fe5a395fb
 md"""
 !!! tip
-	Plotting a function `myfun` is as simple as entering `plot(myfun)`. The same syntax applies if `myfun` is a vector of functions. However, don't forget it was asked to plot only **100** growth curves.
+	Remember: anonymous functions can be defined using `myfun = x -> ...`, and can be visualized using `plot(myfun)`. The same syntax applies if `myfun` is a vector of functions. However, don't forget it was asked to plot only **100** functions.
 """
 
+# ╔═╡ ac37dc49-ff1d-41c6-a2e3-74be2aaabd4c
+logfuns = missing
+
 # ╔═╡ d6c193d7-9fe2-413b-801f-ebc33c772ee9
 missing # plot
 
@@ -130,7 +109,7 @@ md"# 2: Attraction"
 
 # ╔═╡ 431023df-3724-4325-b0ac-96dbf5e4fd20
 md"""
-Following a course on electromagnetism will teach one that computing the net force between 2 arbitrary shapes can be a terrifying task. Tragedy has it then, that this is a very general problem with application from making fusion reactors to space travel. We can ease the pain by turning it into a sampling problem.
+Following a course on electromagnetism will teach one that computing the net force between 2 arbitrary shapes can be a terrifying task. Tragedy has it then, that this is a very general problem with applications from making fusion reactors to space travel. We can ease the pain by turning it into a sampling problem.
 
 We'll start in a humble manner and simulate **the gravitational force between 2 cubes**. Both cubes are size 1. The first cube is in [0, 1] x [0, 1] x [0, 1], and the second cube in [1.1, 2.1] x [0, 1] x [0, 1], as shown in the figure below.
 """
@@ -141,27 +120,28 @@ begin
 	ye = [0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1]
 	ze = [0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1]
 
-	xe2 = [0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0] .+ 1.1
+	xe2 = xe .+ 1.1
 
 	plot(xlims = (-0.5, 2.5), ylims = (-1, 2), zlims = (0, 3))
 	plot!(xe, ye, ze; color = :blue, linewidth = 0.5, label = "cube 1")
 	plot!(xe2, ye, ze; color = :orange, lw = 0.5, label = "cube 2")
 end
 
-# ╔═╡ bb8e4cd6-6106-4d3f-9e35-0dfd8b2c45f5
+# ╔═╡ 4cb55330-8a8a-4622-af4a-ab6f0b643123
 md"""
-The gravitational force can be estimated by **randomly sampling a point from both cubes** and using the formula for gravitational force between those points, ignoring all constants:
+You can sample the gravitational force between both cubes by **randomly sampling a point from both cubes** and using the formula for gravitational force between those points, ignoring all constants:
 
 ```math
-F = \frac{1}{r^2}
+F = \frac{1}{r^2} \, .
 ```
 """
 
 # ╔═╡ ceeeee76-e31c-4429-8ed0-e1c503433dbf
 md"""
 !!! questions
-	1. What is the estimated net force between the two cubes? Is this the same as if you had treated the cubes as point masses?
-	1. How many samples do you need to estimate this force reliably? Define a reliable estimator as one having a standard deviation of 0.1. Visualise the distribution of the estimator.
+	1. What is the estimated total force between the two cubes? Is this the same as if you had treated the cubes as point masses?
+	1. To estimate the net force, you take the average of $n$ samples. Of course, the result will vary every time you take a new sample: if you take the average of only 10 samples, your estimated total force will vary wildly! We can quantify how much the estimated total force varies by taking a **sample of sample averages** and then calculating the variance. Assuming you want your sample average to have a variance of no more than $0.01$, how many samples do you need? Visualise the distribution of the sample average.
+	    - EXTRA: How does the [central limit theorem](https://en.wikipedia.org/wiki/Central_limit_theorem) apply to this question? Can you use it to answer the question with less trial and error?
 """
 
 # ╔═╡ 9be28327-d87f-4d50-bbcc-91e799f14dbf
@@ -169,27 +149,15 @@ md"### 1: Net Force"
 
 # ╔═╡ 06d0e92f-1f07-41fd-b6ee-e94eb539627d
 @model function cubeforce()
-    x1 ~ missing
-    y1 ~ missing
-    z1 ~ missing
-
-	x2 ~ missing
-    y2 ~ missing
-    z2 ~ missing
-
 	F = missing
-	
     return F
 end
 
 # ╔═╡ 43119060-bd35-4c4c-8831-d5c5df1d8dd5
 cubemodel = cubeforce();
 
-# ╔═╡ d26ad946-2bb4-4383-860f-d601903ce1be
-force_sp = missing
-
 # ╔═╡ 7b117453-0875-4b41-a228-866c6c0a8208
-force_average = missing
+estimated_force = missing
 
 # ╔═╡ b28cfbae-2fac-4b38-b234-53f71e381bcd
 pointmass_force = missing # (doesn't require Turing, only maths)
@@ -198,22 +166,18 @@ pointmass_force = missing # (doesn't require Turing, only maths)
 md"### 2: Variance of Estimator"
 
 # ╔═╡ e38e8f51-90db-4136-84e2-f06cd03d502a
-@bind required_samples Slider(10:10:200, show_value = true)
-
-# ╔═╡ b753fa82-28ab-46d1-9085-4c65301db046
-estimator = mean([cubemodel() for _ in 1:required_samples]) 
-	# a single estimation of the force given `required_samples` samples
+@bind n Slider(10:10:200, show_value = true)
 
 # ╔═╡ 55cf127f-5534-4926-955a-487ea9553b70
-estimator_sp = missing 
-	# a sample of estimations given `required_samples` samples
-
-# ╔═╡ f6a67c1e-3677-4e51-b6dd-5d11a5146ea7
-estimator_sp_σ = missing
-	# standard deviation of the force estimator
+force_samples = missing 
+	# multiple samples of your estimated force, using `n` samples
 
 # ╔═╡ 45f0813b-c4c9-4f13-8d66-1e58293c4422
-missing # histogram
+missing # histogram of estimated force samples
+
+# ╔═╡ f6a67c1e-3677-4e51-b6dd-5d11a5146ea7
+force_samples_var = missing
+	# variance of the estimated force
 
 # ╔═╡ Cell order:
 # ╟─52a38b60-178b-4a1d-ac32-e73fafd339f9
@@ -228,31 +192,25 @@ missing # histogram
 # ╟─bebe6f6b-9603-4062-8685-363ed7ff3dcb
 # ╠═eb6dc4e7-e779-4bfc-b865-3defa3894181
 # ╟─107533fc-b300-4b8d-bea2-a3aa6a37938d
-# ╟─38c5b9cd-0baf-4a62-912e-4993614ddbf3
 # ╠═bde57599-1dca-41a4-94aa-498da72c2012
 # ╠═976cfb96-3196-4a61-bd5f-e4f5d24ba1e9
-# ╠═4345b3dd-0731-4dd5-a319-627b4f91306e
-# ╠═9d747eef-a883-49b3-acb7-f0d077a2b902
-# ╠═ce8b53c3-0af1-4e9a-ad80-6fd7a2bf020b
-# ╠═e1d0c5d5-2b7a-4af4-a7ac-e1ca46e665c8
 # ╠═e3092915-a083-425e-8fa1-b7bf370abc8a
 # ╟─e4f42f16-5cce-4fc2-aa01-8971f37c710e
 # ╟─2b39e0bd-91cd-456e-9053-7b8fe5a395fb
+# ╠═ac37dc49-ff1d-41c6-a2e3-74be2aaabd4c
 # ╠═d6c193d7-9fe2-413b-801f-ebc33c772ee9
 # ╟─c8941726-9e81-47fc-9b7e-cb3b5c0c61ca
 # ╟─431023df-3724-4325-b0ac-96dbf5e4fd20
 # ╟─599ac984-ef1d-4c7a-8e87-9d4ddb1aa710
-# ╟─bb8e4cd6-6106-4d3f-9e35-0dfd8b2c45f5
+# ╟─4cb55330-8a8a-4622-af4a-ab6f0b643123
 # ╟─ceeeee76-e31c-4429-8ed0-e1c503433dbf
 # ╟─9be28327-d87f-4d50-bbcc-91e799f14dbf
 # ╠═06d0e92f-1f07-41fd-b6ee-e94eb539627d
 # ╠═43119060-bd35-4c4c-8831-d5c5df1d8dd5
-# ╠═d26ad946-2bb4-4383-860f-d601903ce1be
 # ╠═7b117453-0875-4b41-a228-866c6c0a8208
 # ╠═b28cfbae-2fac-4b38-b234-53f71e381bcd
 # ╟─304d6052-6d3d-487c-8c01-dab259276d6f
 # ╠═e38e8f51-90db-4136-84e2-f06cd03d502a
-# ╠═b753fa82-28ab-46d1-9085-4c65301db046
 # ╠═55cf127f-5534-4926-955a-487ea9553b70
-# ╠═f6a67c1e-3677-4e51-b6dd-5d11a5146ea7
 # ╠═45f0813b-c4c9-4f13-8d66-1e58293c4422
+# ╠═f6a67c1e-3677-4e51-b6dd-5d11a5146ea7
diff --git a/exercises/student_notebooks/P4_probmod/probmod_4-review.jl b/exercises/student_notebooks/P4_probmod/probmod_4-review.jl
index f2b9423..c12d084 100644
--- a/exercises/student_notebooks/P4_probmod/probmod_4-review.jl
+++ b/exercises/student_notebooks/P4_probmod/probmod_4-review.jl
@@ -1,5 +1,5 @@
 ### A Pluto.jl notebook ###
-# v0.20.4
+# v0.20.21
 
 using Markdown
 using InteractiveUtils
@@ -15,7 +15,7 @@ md"# Review exercise: Buffon's needles"
 
 # ╔═╡ a6435b94-f1af-4609-abfc-93d88730d023
 md"""
-A wise man once said: ["there is no greater joy than estimating π"](https://en.wikipedia.org/wiki/Approximations_of_%CF%80). Next to throwing darts at the unit square, another method to accomplish this is using [Buffon's needle problem](https://en.wikipedia.org/wiki/Buffon%27s_needle_problem).
+A wise man once said: ["there is no greater joy than estimating π"](https://en.wikipedia.org/wiki/Approximations_of_%CF%80). One method to accomplish this is using [Buffon's needle problem](https://en.wikipedia.org/wiki/Buffon%27s_needle_problem).
 
 The experiment is as follows: consider a floor with parallel lines all a distance of 1 away from eachother. Now drop a needle of length 1 (and width ~0) on the floor with a **random position and angle**. What is the probability $P_{cross}$ that the needle will cross one of the lines?
 """
@@ -24,8 +24,8 @@ The experiment is as follows: consider a floor with parallel lines all a distanc
 md"The following image illustrates the problem (imagine $l$ = $t$ = 1) for two needles, where `a` crosses a line and `b` does not."
 
 # ╔═╡ c2ba2183-957f-4b94-a22e-0c2f3dd957ad
-md"""
-![Buffon's needles](https://upload.wikimedia.org/wikipedia/commons/thumb/5/58/Buffon_needle.svg/1920px-Buffon_needle.svg.png)
+html"""
+<img src="https://upload.wikimedia.org/wikipedia/commons/thumb/5/58/Buffon_needle.svg/1920px-Buffon_needle.svg.png" style="background-color:white;" alt:"Buffon's needles">
 """
 
 # ╔═╡ d71cf8dc-99e5-48e0-9abe-2242a6ccc30b
diff --git a/exercises/student_notebooks/P5_mcmc/MCMC_1-intro.jl b/exercises/student_notebooks/P5_mcmc/MCMC_1-intro.jl
index 46f856a..c40109a 100644
--- a/exercises/student_notebooks/P5_mcmc/MCMC_1-intro.jl
+++ b/exercises/student_notebooks/P5_mcmc/MCMC_1-intro.jl
@@ -281,10 +281,10 @@ Taking the sampled values of the mutation rate from the chain and plotting a his
 """
 
 # ╔═╡ b657217a-6ccc-4a41-b852-df4e39a7a10a
-sp_alpha = mutation_chain[:α];
+alpha_samples = mutation_chain[:α];
 
 # ╔═╡ 9cf08616-d599-42a3-82e8-e8c98853c1d8
-histogram(sp_alpha) # note the difference with the prior distribution!
+histogram(alpha_samples) # note the difference with the prior distribution!
 
 # ╔═╡ e02c42dc-627c-4bc0-8ebb-fdb2b0f15b64
 md"Plotting some sampled mutation rates from this distribution onto our data shows that they fit well:"
@@ -295,17 +295,17 @@ begin
 		ylabel = "Number of mutations", label = false, xlims = (0, 500), 
 		ylims = (0, 400), title = "A posteriori relationship between t and mean N"
 	);
-	plot!([x -> αᵢ*x for αᵢ in sp_alpha[1:10:end]], color = :purple, opacity = 0.1, label = false)
+	plot!([x -> αᵢ*x for αᵢ in alpha_samples[1:10:end]], color = :purple, opacity = 0.1, label = false)
 end
 
 # ╔═╡ 87059440-2919-4f6d-9d32-0df3ce75e2a2
 md"And finally we can answer our first question:"
 
 # ╔═╡ 644cff58-68b9-4d4b-8896-617fcacc39c5
-mean(sp_alpha)
+mean(alpha_samples)
 
 # ╔═╡ 5f2f5a78-ca90-4baf-8bf0-ff1fae88d785
-sqrt(var(sp_alpha))
+sqrt(var(alpha_samples))
 
 # ╔═╡ bfbf811f-7b2f-473b-b132-0c7e965c8b0d
 md"""
diff --git a/exercises/student_notebooks/P5_mcmc/MCMC_2-basics.jl b/exercises/student_notebooks/P5_mcmc/MCMC_2-basics.jl
index 8c9f21e..e49044f 100644
--- a/exercises/student_notebooks/P5_mcmc/MCMC_2-basics.jl
+++ b/exercises/student_notebooks/P5_mcmc/MCMC_2-basics.jl
@@ -71,7 +71,7 @@ end
 molemodel = mole()
 
 # ╔═╡ 61b8ef63-a613-4dcb-8f7a-936e0d59b862
-spY = missing
+Y_samples = missing
 
 # ╔═╡ 4e8b9000-cb61-4f01-9ba9-17276ad0335e
 E_Y = missing
@@ -89,7 +89,7 @@ molechain = missing
 missing # plot molechain
 
 # ╔═╡ af4811ce-c6e7-458a-8f89-0562355b3eb0
-spXcondY = missing
+X_samplescondY = missing
 
 # ╔═╡ fd540765-95e5-4071-81f7-e689b06cad0c
 E_XcondY = missing