import numpy as np import matplotlib.pyplot as plt if __name__ == "__main__": n_dim = 10 # We will sample a standard Gaussian random vector X in n_dim-dimensional Euclidean space. n_samples = 100000 # We will sample it n_samples times. r_min = 0.0 # We will look at the distribution of the Euclidean 2-norm r of X between r_min and r_max, in bins of width delta_r r_max = 3.0 delta_r = 0.10 bins = np.arange(r_min, r_max + delta_r, delta_r) counts = np.zeros((len(bins) - 1,)) # Set up an array to count how many samples land in each bin for n in range(n_samples): x = np.random.normal(size=(n_dim,)) # Draw a standard Gaussian r = np.linalg.norm(x) # Calculate its Euclidean norm normalised_r = r / np.sqrt(n_dim) # We expect concentration around the sphere of radius n^{1/2}, so let's normalise # Now we check into which bin the sample draw falls, and increment the appropriate counter for i in range(len(bins))[:-1]: if normalised_r >= bins[i] and normalised_r < bins[i + 1]: counts[i] += 1.0 masses = counts / n_samples # Divide by the number of samples to get (empirical) probability masses # And now plot the results fig, ax = plt.subplots() histo = ax.bar(bins[1:], masses, width=delta_r) ax.set_xlabel("|| X ||_2") ax.set_ylabel("pdf") ax.set_title("%s-sample PDF for || X ||_2 with X ~ N(0, I_%s)" % (n_samples, n_dim)) plt.show() # EOF