1import numpy as np
2from matplotlib import pyplot as plt
3import seaborn as sns
4import pandas as pd
5
6mu, std = 40.0, 5.0 # Mean and standard deviation of oxygen partial pressure (p)
7
8
9def map_samples(samples):
10 """
11 **TODO**: Transform samples of a normally distributed oxygen partial pressure (p)
12 using the Hill equation to approximate blood oxygen saturation.
13
14 .. math::
15 S(p) = \\frac{p^n}{p^n + K^n}
16
17 Note: You need to estimate the mean and variance of the samples by using
18 `np.mean <https://numpy.org/doc/2.2/reference/generated/numpy.mean.html>`_ and
19 `np.var <https://numpy.org/devdocs/reference/generated/numpy.var.html>`_.
20
21 Additionally compute an analytical approximation of the mean and variance
22 of the transformed values via first-order Taylor expansion.
23
24 :param samples:
25 A 1D array of sampled oxygen partial pressures (in mmHg), assumed to be
26 normally distributed with unknown mean and variance.
27
28 Returns:
29 --------
30 :return: 3-Tuple of mapped_samples, the estimated mean and the estimated variance of the saturation (via Taylor expansion)
31 """
32 N, K = 2.7, 26
33
34 # TODO: Apply the Hill function to each sample to get saturation values
35 mapped_samples = np.power(samples, N) / (np.power(samples, N) + np.power(K, N))
36
37 # TODO: Use np.mean and np.var to compute mean and variance of the original input samples
38 mu = np.mean(samples)
39 V = np.var(samples)
40
41 # TODO: Compute the derivative (Jacobian) of the Hill function at the mean
42 J = (
43 N
44 * np.power(K, N)
45 * np.power(mu, N - 1)
46 / ((np.power(mu, N) + np.power(K, N)) ** 2)
47 )
48
49 # TODO: Compute the transformed mean using the Hill function
50 mapped_mu = np.power(mu, N) / (np.power(mu, N) + np.power(K, N))
51
52 # TODO: Approximate the transformed variance via linear error propagation
53 mapped_var = J * V * J
54
55 # TODO: Return the mapped samples, the mapped mean and the mapped variance
56 return mapped_samples, mapped_mu, mapped_var
57
58
59# ---------------------------------------------------
60# There is no need to change anything below this line
61# ---------------------------------------------------
62
63# Main Program
64if __name__ == "__main__":
65 N, K = 2.7, 26
66
67 # Create three plots
68 fig, axs = plt.subplots(ncols=3, figsize=(8, 4))
69
70 # Draw random samples for the partial pressure
71 p = np.random.normal(mu, std, 10000)
72
73 # Plot histogram and density curve
74 df = pd.DataFrame(p, columns=["p"])
75 sns.histplot(df, bins=64, ax=axs[0], stat="density")
76 x = np.linspace(mu - 5.0 * std, mu + 5.0 * std, 100)
77 y = np.exp(-((x - mu) ** 2) / (2 * std**2)) / (std * np.sqrt(2.0 * np.pi))
78 axs[0].plot(x, y, color="#1f77b4")
79 axs[0].set_title("Sauerstoffpartialdruck - Verteilung")
80
81 # Draw the mapping function itself
82 x = np.linspace(0.0, 5.0 * K, 250)
83 y, _, _ = map_samples(x)
84 axs[1].plot(x, y, color="#1f77b4")
85 axs[1].set_ylim(0.0, 1.0)
86 axs[1].set_title("Hill-Gleichung")
87
88 # Draw tangent and range which is used from the samples
89 y0 = np.power(mu, N) / (np.power(mu, N) + np.power(K, N))
90 grad = (
91 N
92 * np.power(K, N)
93 * np.power(mu, N - 1)
94 / ((np.power(mu, N) + np.power(K, N)) ** 2)
95 )
96 y = (x - mu) * grad + y0
97 axs[1].plot(x, y, color="orange")
98 axs[1].plot(np.array([mu - 3 * std, mu - 3 * std]), np.array([0.0, 1.0]), "k--")
99 axs[1].plot(np.array([mu + 3 * std, mu + 3 * std]), np.array([0.0, 1.0]), "k--")
100
101 # Map the samples and retrieve estimated mean and variance
102 s, mu, var = map_samples(p)
103
104 # Plot histogram of mapped samples and estimated density
105 df = pd.DataFrame(s, columns=["s"])
106 sns.histplot(df, bins=64, ax=axs[2], stat="density")
107 x = np.linspace(mu - 5.0 * np.sqrt(var), mu + 5.0 * np.sqrt(var), 100)
108 y = np.exp(-((x - mu) ** 2) / (2 * var)) / (np.sqrt(2.0 * np.pi * var))
109 axs[2].plot(x, y, color="#1f77b4")
110 axs[2].set_title("Sauerstoffsättigung - Verteilung")
111
112 axs[2].plot(np.array([mu, mu]), np.array([0.0, np.max(y)]), "k--")
113
114 mu = np.mean(s)
115 axs[2].plot(np.array([mu, mu]), np.array([0.0, np.max(y)]), "o--")
116
117 plt.show()