|
15 | 15 | # See the License for the specific language governing permissions and
|
16 | 16 | # limitations under the License.
|
17 | 17 |
|
| 18 | +""" |
| 19 | +This script is meant to be copy-pasted into an IPython REPL session. |
| 20 | +
|
| 21 | +Upon execution, enter `X` to print the contents of the generated array of random |
| 22 | +values. |
| 23 | +
|
| 24 | +To generate an array containing the biased covariance, set the `bias` kwarg to |
| 25 | +`True`. |
| 26 | +""" |
| 27 | + |
| 28 | +import numpy as np |
| 29 | + |
18 | 30 | # Define the list of known means:
|
19 |
| -mean = [ 1.0, 1.0 ]; |
| 31 | +mean = [1.0, 1.0] |
20 | 32 |
|
21 | 33 | # Create a 2x3 matrix of mean values:
|
22 |
| -broadcasted_mean = np.broadcast_to(mean, (3,2)).T; |
| 34 | +broadcasted_mean = np.broadcast_to(mean, (3, 2)).T |
23 | 35 |
|
24 | 36 | # Define the matrix of standard deviations:
|
25 | 37 | sigma = [
|
26 |
| - [ 1.0, 0.7 ], |
27 |
| - [ 0.7, 1.0 ] |
28 |
| -]; |
| 38 | + [1.0, 0.7], |
| 39 | + [0.7, 1.0] |
| 40 | +] |
29 | 41 |
|
30 | 42 | # Generate a random sample of normally distributed numbers:
|
31 |
| -X = np.random.default_rng().multivariate_normal(mean, sigma, 3).T; |
| 43 | +X = np.random.default_rng().multivariate_normal(mean, sigma, 3).T |
32 | 44 |
|
33 | 45 | # Center the generated values:
|
34 | 46 | for n in range(X.shape[0]):
|
35 | 47 | X[n] = X[n] - X[n].mean()
|
36 | 48 |
|
37 | 49 | # Transform the generated values via Cholesky decomposition (see https://stats.stackexchange.com/questions/120179/generating-data-with-a-given-sample-covariance-matrix and https://www.r-bloggers.com/2011/10/simulating-data-following-a-given-covariance-structure/)...
|
38 |
| -L_inv = np.linalg.inv(np.linalg.cholesky(np.cov(X, bias = False))); |
39 |
| -X = np.dot(L_inv, X); |
| 50 | +L_inv = np.linalg.inv(np.linalg.cholesky(np.cov(X, bias=False))) |
| 51 | +X = np.dot(L_inv, X) |
40 | 52 |
|
41 |
| -L = np.linalg.cholesky(sigma); |
42 |
| -X = np.dot(L, X); |
| 53 | +L = np.linalg.cholesky(sigma) |
| 54 | +X = np.dot(L, X) |
43 | 55 |
|
44 | 56 | # Un-center the transformed values:
|
45 |
| -X += broadcasted_mean; |
46 |
| - |
47 |
| -# Display the generated values: |
48 |
| -X |
| 57 | +X += broadcasted_mean |
49 | 58 |
|
50 | 59 | # Confirm that the generated values have the expected covariance matrix:
|
51 |
| -np.cov(X, bias = False) |
| 60 | +np.cov(X, bias=False) |
0 commit comments