Monte Carlo Simulations in Python

I wrote a few lines of code below in python. My intentions are to generate 1000 different values of 'φ' and 'Λ' and in each case it should choose a new value of 'a' from a random range (0,1).

Then convert the values "φ" and "Λ" from cylindrical coordinates to cartesian coordinates and show the plot in the graph.

How do I generate different values of "φ" and "Λ" in this case? And I tried to convert to Cartesian but it rather shows an empty plot. Please what do I need to do at this point in order to arrive at my results.

I sincerely appreciate your past contributions and I look forward to your help, thanks.

# import the necessary libraries
import math
import random
from math import pi

N = 1000

for num in range(N):
a = random.uniform(0,1)
# I am using these symbols for: λ - lambda null, Λ - mean free path, φ = angle
λ = -1
φ = 2*pi*a
Λ = λ*math.log(a)

x = Λ*math.cos(φ)
y = Λ*math.sin(φ)
plt.plot(x,y)
plt.show()

• Try printing x and y, and not just plotting them. You may find the problem is with matplotlib.pyplot, not the rest.
– J.G.
Sep 26, 2022 at 14:23
• little suggestion: do not use special characters in scripting or coding. Use "phi" instead of "$\phi$" Sep 26, 2022 at 14:26
• I've tried your code. x and y are each only one number, and that's not what you want. If you replace each math and random function with its numpy counterpart, or make lists and append to them in the for loop, it'll work.
– J.G.
Sep 26, 2022 at 14:28
• Thank you very much @J.G for your valuable inputs. I have printed x and y as you rightly suggested, I got a single values for x and y. My intention is to generate 1000 different values of 'phi' and 'lambda' and to plot them. Please how may I do that? Sep 26, 2022 at 14:31
• The problem with your for loop is that only the line a = random.uniform(0,1) is inside it. What comes after it must be indented to match the columns of a. Effectively you see the last value of a in your code so far. I think only plt.show() can stay on the left where it is. Sep 26, 2022 at 17:49

2 Answers

Are you sure your transformation is correct? You are using the single random variable $$a$$ to generate points $$\phi,\,\lambda$$, which means that you are tracing a path along a log-polar coordinate system (not 2D cylindrical, as you are using $$\log(a)$$ instead of $$\sqrt{a}$$). For instance, running JG's code (but using scatter() instead of plot()) I get,

Which seems wrong to me, as you are doing a transformation $$\mathbb{R}^1\to\mathbb{R}^2$$. I would think you need a transformation of $$\mathbb{R}^2\to\mathbb{R}^2$$ in order to plot a randomly/uniformly distributed series of points on the domain. This would be done by using two arrays of random values,

import matplotlib.pyplot as plt
import numpy as np

N = 1000
a = np.random.uniform(0, 1, size=N)
b = np.random.uniform(-1, 1, size=N)
lam = -np.log(a)
phi = np.arccos(b)
plt.scatter(lam * np.cos(phi), lam * np.sin(phi))
plt.plot()


which gives a plot more like,

Note that using $$\phi=\cos^{-1}(x)$$ for $$x\in\operatorname{Uniform}(-1,1)$$ is for removing biases from sampling $$\operatorname{Uniform}(0,\,2\pi)$$ (cf. this site, among others for sampling the hemisphere).

• Thank you very much @Kyle Kanos for your useful inputs. They are greatly appreciated. Some other concerns are, (a) the scattering points only concentrates between (0 and pi) that is from 0 to 180 degrees. How do I make the scattering symmetrical such that the points are all scattered in all directions covering between (0 to 2Pi)? (b) How do I also show the distances between one point to another in a visible chat? Thank you and everyone for your valuable contributions. I truly appreciate them all. Sep 28, 2022 at 0:57
• For (a), I would look at the site I linked at the bottom of the answer as it indicates how one needs to think about such transformations. For (b), I am not sure what you are asking. Do you mean how to pick two points and connect them via a straight line in the chart? Or simply how to calculate the distances and display them? Sep 28, 2022 at 12:16
• Thank you again for your continued help. I have checked this site blog.thomaspoulet.fr/uniform-sampling-on-unit-hemisphere and other references but I have not been able to arrive at the desired answer yet. As for, the (b) part, I understand how to do it now. Thank you very much for every help. Sep 28, 2022 at 22:29

Expanding on my comments while also taking advice from @basics, you want something like

from matplotlib import pyplot as plt
import numpy as np

a = np.sort(np.random.uniform(size=1000))
Lambda, phi = -np.log(a), 2*np.pi*a
plt.plot(Lambda*np.cos(phi), Lambda*np.sin(phi))
plt.show()


A numpy-free solution would be something like

from matplotlib import pylot as plt
import math, random

a = sorted(random.random() for _ in range(1000))
Lambda, phi = [-math.log(x) for x in a], [2*math.pi*x for x in a]
plt.plot(*zip(*((L*math.cos(p), L*math.sin(p)) for L, p in zip(Lambda, phi))))
plt.show()