1import numpy as np
2from matplotlib import pyplot as plt
3import seaborn as sns
4from misc import draw_samples, draw_cov_ellipses
5
6STD_ALPHA = 5 # Standard deviation of angular values (degrees)
7
8
9def map_samples(samples, radius, alpha, cov):
10 """
11 **TODO**
12 Assume you have a normal distributed random variable :math:`\\boldsymbol{X}`
13 with mean
14
15 .. math::
16 \mu = (r, \\alpha)
17
18 and given covariance (`cov`-Parameter). :math:`\\boldsymbol{X}` represents
19 a random point in polar coordinates. We want to transform that point into cartesian coordinates by
20 applying the following transformation
21
22 .. math::
23 \\boldsymbol{Y} = \\begin{pmatrix}
24 r\cdot\cos(\\alpha)\\\\
25 r\cdot\sin(\\alpha)
26 \\end{pmatrix}
27
28 The samples parameter holds 512 samples of this random variable.
29
30 Apply the non-linear mapping to the samples and calculate the **exact** new mean and covariance of
31 :math:`\\boldsymbol{Y}` after linearization using the Jacobian Matrix ::math:`J`, namely
32
33 .. math::
34 E[\\boldsymbol{Y}] = f(\\boldsymbol{\mu})
35
36 .. math::
37 Cov[\\boldsymbol{Y}] = J\cdot Cov[\\boldsymbol{X}] \cdot J^T
38
39 Return the mapped samples as well as the
40 derived mean and covariance of the mapped random variable.
41 **Do not** estimate the mean and covariance from the mapped samples.
42
43 :param samples: (np.array 2x512) 512 Samples from :math:`\\boldsymbol{X}`
44 :param radius: The mean radius in polar coordinates
45 :param alpha: The mean angle in polar coordinates
46 :return: 3-tuple (mapped_samples, mapped_mu, mapped_cov)
47 """
48 # TODO: Map all samples from polar to cartesian coordinates
49 r, a = samples[0, :], samples[1, :]
50
51 mapped_samples = np.zeros_like(samples)
52 mapped_samples[0, :] = r * np.cos(a)
53 mapped_samples[1, :] = r * np.sin(a)
54
55 # TODO: Calculate Jacobian
56 J = np.array(
57 [
58 [np.cos(alpha), -radius * np.sin(alpha)],
59 [np.sin(alpha), radius * np.cos(alpha)],
60 ]
61 )
62
63 # TODO: Calculate mean and covariance of Y after linearization
64 mapped_mu = np.array([radius * np.cos(alpha), radius * np.sin(alpha)])
65 mapped_cov = J @ cov @ J.T
66
67 # TODO: Return your mapped samples, the mapped mean and the mapped covariance
68 return mapped_samples, mapped_mu, mapped_cov
69
70
71# ---------------------------------------------------
72# There is no need to change anything below this line
73# ---------------------------------------------------
74
75# Main Program
76if __name__ == "__main__":
77 # Set Seaborn display style
78 sns.set_style("whitegrid")
79
80 # Create figure for plots
81 fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(8, 8))
82
83 # Iterate different radians for plotting (0 to 2pi)
84 for radians in np.linspace(-np.deg2rad(45), np.deg2rad(45), 90):
85 # Generate 512 samples of a multivariate normal random variable (Shape 2 x 512)
86 cov = np.array([[0.5, 0.0], [0.0, np.deg2rad(STD_ALPHA) ** 2]])
87
88 samples = np.random.multivariate_normal(
89 mean=np.array([50, radians]), cov=cov, size=512
90 ).T
91
92 # Clear all axes
93 for a in ax:
94 a.clear()
95
96 # Generate mapped samples
97 mapped_samples, mapped_mu, mapped_cov = map_samples(samples, 50.0, radians, cov)
98
99 # Estimate mu and covariance from original samples as well as mapped samples
100 mu = np.mean(samples, axis=1)
101 cov = np.cov(samples)
102
103 estimated_mu = np.mean(mapped_samples, axis=1)
104 estimated_cov = np.cov(mapped_samples)
105
106 # Draw original sample point cloud as well as mapped samples point cloud
107 draw_samples(samples, ax[0])
108 draw_samples(mapped_samples, ax[1])
109
110 # Draw covariance ellipses
111 draw_cov_ellipses(
112 mu, cov, ax[0], edgecolor="lightblue", facecolor="none", linewidth=2
113 )
114 draw_cov_ellipses(
115 estimated_mu,
116 estimated_cov,
117 ax[1],
118 edgecolor="lightblue",
119 facecolor="none",
120 linewidth=2,
121 )
122 draw_cov_ellipses(
123 mapped_mu,
124 mapped_cov,
125 ax[1],
126 edgecolor="#1f77b4",
127 facecolor="none",
128 linewidth=2,
129 )
130
131 # Apply styling to plot
132 ax[0].set_title("Polar Coordinates")
133 ax[1].set_title("Cartesian Coordinates")
134
135 ax[0].set_xlabel("r")
136 ax[0].set_ylabel("alpha")
137 ax[0].set_xlim(20.0, 80.0)
138 ax[0].set_ylim(-np.deg2rad(60), np.deg2rad(60))
139
140 ax[1].set_xlim(20.0, 80.0)
141 ax[1].set_ylim(-45.0, 45.0)
142
143 # Wait shortly for animation to roll
144 plt.pause(0.1)
145 # plt.show()