We’ll take a look at the design and execution of a simple Monte-Carlo simulation using Python.
A usually boring train commute of mine turned exciting when one of my friends challenged me to solve a problem from the prestigious Putnam Competition — a test that perplexes even the smartest math undergraduates out there.
The problem was as follows (paraphrased):If you choose 4 random points on a sphere and join them, what is the probability that the resulting tetrahedron contains the centre of the sphere?Now it’s certainly possible to ascertain the exact solution by applying a few concepts of geometry and probability.
I expect the math geniuses among you to be ready with the answer before you even finish reading this sentence.
For those of you (like me) who can’t, a video by 3Blue1Brown on YouTube provides a pretty elegant explanation of the analytical solution.
The word “random” in the problem, however, sparked off an alternate train of thought in my mind that eventually led to the solution that I’m going to present.
I considered the possibility of obtaining the solution through random sampling.
The idea was simple — draw a large number of random samples of tetrahedra and calculate the probability considering all those samples.
The AlgorithmThe broad algorithm for the simulation is as follows:Define the surface of a sphere of known radius as the sample space.
Define the resultant probability as num/den where num=0 and den=1.
Draw one sample of a tetrahedron inside the sphere from the sample space.
Determine whether the centre lies inside the tetrahedron.
Calculate the resultant probability by adding either a 1 or a 0 to num based on the success or failure of the current trial respectively and adding a 1 to den in any case.
Repeat steps 3 to 5 for the specified number of random samples to be chosen.
Display the resultant probability.
The CodeWe shall go step-by-step through the code and examine each step in detail.
Step 1: Import Required ModulesFirst, we import the required modules in Python — namely Scipy and Matplotlib.
Scipy has a sub-module that allows generation of random numbers.
import scipy as sciimport matplotlib.
pyplot as pltStep 2: Define Sample SpaceThe next step involves creating a sample space composed of a set number of points.
For this example, let the sample space have 10,000 points on the surface of a sphere having radius 1 unit.
How do we ensure that the points on the sample space lie on the sphere?.Three simple steps:Pick out any random 3D vector using sci.
random.
rand(1,3).
This will yield a 3D vector centered around 0.
5 with each coordinate having a value between 0 to 1.
Subtract 0.
5 from each coordinate so that it becomes centered around 0.
Calculate the norm (or magnitude) of the vector and divide the vector by its norm ensuring that it becomes a unit vector.
Any unit vector will necessarily lie on a sphere of radius 1.
An if-condition that ensures vectors with norms equal to 0 or greater than 1 do not undergo this transformation can be added for safety.
This animation was created using matplotlib; the code for this is not given in the article.
It is merely used for the purpose of representation of under-the-hood computations.
points=10000 #Number of points x=sci.
zeros((points,3)) #To store x-coordinates of sample pointsy=sci.
zeros((points,3)) #To store y-coordinates of sample pointsfor i in range(points): vector=sci.
random.
rand(1,3)-0.
5 if(sci.
linalg.
norm(vector)!=0 and sci.
linalg.
norm(vector)<=1.
0): x[i,:]=vector y[i,:]=x[i,:]/sci.
linalg.
norm(x[i,:])To be absolutely sure that all the points in the sample space lie on the surface of the sphere, we can add another checkpoint.
We compute the difference between the actual norms of the vectors in the sample space and the desired norm (equal to the radius of the sphere; here 1) and compare it with a tolerance (arbitrarily defined; here 1×10⁻¹⁰).
The norms that don’t match the desired norm are stored in danger_array.
If the size of danger_array is 0 (that is, no undesired norms exist), then we display the “All Clear” message and move ahead.
y_norms=sci.
zeros(points) #Array to store norms of sample vectorsfor i in range(points): y_norms[i]=sci.
linalg.
norm(y[i,:]) tol=1e-10 #Tolerance norm_diff=abs(y_norms-1) #Diff.
between actual and desired norm danger_array=y_norms[norm_diff>tol]if(len(danger_array)==0): print("All Clear")else: print("Danger")Step 3: Create functions to determine whether the centre lies inside a tetrahedronOnce all checks have been passed, we move on to the heart of the problem— determining whether the centre of the sphere lies in a random tetrahedron.
This can be simplified by using two functions.
The first checks whether the centre of the sphere and the fourth vertex of the tetrahedron lie on the same side of the plane formed by the remaining three vertices of the tetrahedron.
The second calls the first function for all the four faces of the tetrahedron and concludes whether the centre lies inside the tetrahedron.
To formulate the CheckSide function, we rely on the following principles of linear algebra:The cross product of two vectors gives a vector normal to the plane formed by the two.
The dot product of two vectors gives the magnitude of the component of one of the vectors in the direction of the other.
If the component is oriented in the opposite direction as compared to the second vector, the dot product comes out to be negative.
We use the first principle to find the normal to the plane formed by any three vertices.
Next, we find the dot product of the normal and a vector joining the fourth vertex with any of the three vertices — the sign of this dot product is stored.
Finally, we find the dot product of the normal and a vector joining the centre of the sphere and any of the three vertices — the sign of this dot product is compared with the previous one and if they match, then it can be concluded that the fourth vertex and the centre lie on the same side of the plane formed by the remaining three vertices.
This schematic will make it easier for you to visualize the application of the two principlesdef CheckSide(vertices,point): t1,t2,t3,t4=vertices p=point side_1=t2-t1 side_2=t3-t1 normal=sci.
cross(side_1,side_2) ref_vector=t4-t1 ref_sign=sci.
dot(normal,ref_vector) point_vector=p-t1 point_sign=sci.
dot(normal,point_vector) if(sci.
sign(ref_sign)==sci.
sign(point_sign)): return 1 else: return 0If we apply the CheckSide function to all the four faces of a tetrahedron, we will be able to accurately conclude whether the centre lies inside it.
Thus, we define a new function CheckTetrahedron that calls the CheckSide function four times.
Between every function call, we roll the vertices array by one so that the CheckSide function works on a new face on every next call.
Finally, the values returned by all the four CheckSide calls are added up and if they equal 4, then the centre is concluded to lie inside the tetrahedron; if not, it lies outside.
def CheckTetrahedron(vertices,point): vert=sci.
copy(vertices) check_1=CheckSide(vert,point) vert=sci.
roll(vert,1,axis=0) check_2=CheckSide(vert,point) vert=sci.
roll(vert,1,axis=0) check_3=CheckSide(vert,point) vert=sci.
roll(vert,1,axis=0) check_4=CheckSide(vert,point) sum_check=check_1+check_2+check_3+check_4 if(sum_check==4.
): return 1 else: return 0Step 4: InitializationNow that we have all the requisite functions defined, it’s time to turn to the actual run.
Before we begin the run, we need to initialize a few parameters for the purposes of the run and for plotting.
The array check_point will store ‘0’s and ‘1’s to denote the success or failure of every iteration.
Instead of storing just a single value of resultant probability, we shall create an array to store its value for every iteration so that we can plot a curve that shows the variation in probability with the number of iterations.
centre=[0,0,0]number_of_samples=10000sample_span=sci.
arange(0,number_of_samples,1)check_point=sci.
zeros(number_of_samples) prob=sci.
zeros(number_of_samples)Step 5: The RunIt’s time to run the simulation!.In every iteration, we do the following three things:Pick out 4 different random points from the sample space.
Apply the CheckTetrahedron function to those 4 points and the centre to determine whether the centre lies in the tetrahedron formed by those 4 points.
Depending upon the success or failure of the trial, store either 1 or 0 respectively in check_point and update the value of resultant probability in the prob array.
This animation was created using matplotlib; the code for this is not given in the article.
It is merely used for the purpose of representation of under-the-hood computations.
for i in range(number_of_samples): indices=sci.
random.
randint(0,points,4) vertex_list=y[indices] check_point[i]=CheckTetrahedron(vertex_list,centroid) prob[i]=len(check_point[check_point==1.
])/(i+1)Step 6: PlottingFinally, its time to plot the results and find out the resultant probability after running the simulation for the given number of iterations.
We add a few bells and whistles to the plot to make it look attractive and also to ensure that it conveys all the required information lucidly.
#Plot blank figureplt.
figure(figsize=(15,10))#Plot resultant probability from simulationplt.
plot(sample_span,prob,color="navy",linestyle="-",label="Resultant probability from simulation")#Plot resultant probability from analytical solutionplt.
plot(sample_span,[0.
125]*len(sample_span),color="red",linestyle="-",label="Actual resultant probability from analytical solution (0.
125)")#Plot value of final resultant probability in textplt.
text(sample_span[int(number_of_samples/2)],0.
05,f"The final probability is {prob[-1]:.
4f}",fontsize=16)#Display axis labelsplt.
xlabel("Number of iterations",fontsize=14)plt.
ylabel("Probability",fontsize=14)#Display legendplt.
legend(loc="upper right",fontsize=14)#Display title of the plotplt.
title("Variation of resultant probability with increase in the number of iterations",fontsize=14)The final plot is given below.
From the nature of the plot, it is quite clear that the probability converges to the true value (that is, the analytical solution) as the number of iterations goes up.
Plot generated by matplotlib detailing the variation of probability with the number of iterationsConcluding ThoughtsIt is possible that your simulation may not converge to the true value in the same number of iterations that mine did due to the associated randomness in the method.
In that case, you can keep increasing the number of iterations till your solution converges to the true value.
The analytical solution is 1/8 or 0.
1250 and the Monte-Carlo simulation converged to a value of 0.
1248, which is quite close enough.
We have created an apparent order out of randomness through the simulation — the fact that repeated random sampling can lead to a single convergent solution is astounding.
The method may seem inelegant — almost like a brute-force use of computational power — but it is very effective and does give accurate answers.
The code has been documented and executed in a Jupyter notebook that can be viewed online.
I hope I have been able to explain the method and code with clarity.
You are welcome to respond in the comments, ask questions or provide any suggestions that you might have.
You can also follow me on Twitter here.
.