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