Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am trying to learn how to script/program in mathematica and I am retyping some simple C and FORTRAN programs into Mathematica as a way to learn. I am writing a Monte Carlo code to calculate the value of pi. In the method attached a 10000x2 dimensional array is created using table with values ranging from 0 to 1. The magnitude of each row of the array is determined inside a For loop and if it is greater than 0 it is rejected and if it is less it is used in the running tally that is divided by the total number of points and then multiplied by 4 to give an estimate of pi. On top of calculating pi, I would like to create two more arrays to use in a list plot. In a second array I would like to store values of the first array that met the rejection criteria and in third array I would like to store values that were rejected, so I can plot the rejected points as red and the accepted points as blue. However, no matter what I try Mathematica will not seem to let me transfer points from the first array titled "vector". I am attaching the current code to this question that effectively calculates the estimate of pi. If anyone can show me how to transfer the rejected points to another array as determined inside the For loop I would be very grateful.

[counter, samples, a, vector, tally, pi, tally2] counter = 1; samples = 10000; tally = 1; a := Random[Real, {0, 1}]; Timing[vector = Table[a, {i, samples}, {j, 2}]; For[i = 1, i <= samples, i++, {If[Norm[vector[[i]]] < 1, tally++, ""]}]; pi = (tally/samples)*4 // N] `

share|improve this question
    
Please take a look at this mathematica.stackexchange.com/q/18393/193, and this mathematica.stackexchange.com/q/18/193. You're using the wrong paradigms ... –  belisarius Oct 25 '13 at 5:06

2 Answers 2

up vote 2 down vote accepted

This should get you started:

Generate your 2D points:

arr = RandomReal[1, {100000, 2}];

From there Select points that fall within a quarter of a circle:

cir = Select[arr, Norm[#] < 1 &];

Now compute Pi

N[4 Length[cir]/Length[arr]]

This gives:

3.14292

As you increase the number of points, your approximation of Pi will be better.

You can even plot the two regions as follows:

Graphics[{If[Norm[#] < 1, Red, Darker[Green]], Point[#]} & /@ arr]

Mathematica graphics

share|improve this answer
    
Thank you, I was not even aware that the Select function even existed much less how to use it. I also did not even think to use the length function. Your comment was a big help in my learning process, thank you. –  Jon Oct 25 '13 at 5:40
    
@Jon. Glad I could help. Note that you can up vote an answer if you like it and you can also select one of the answers that you feel was most helpful to you. Welcome to Mathematica stackexchange. –  RunnyKine Oct 25 '13 at 5:51
    
I will do that. Would you mind also helping solve a minor mystery for me? I know the general theory behind the # and & symbols but there are a few things perplexing me in your example. As I understand it the code will reference a # to a & symbol and the variable to the right of the & symbol. However, in your second line of code there is nothing to the right of the & symbol. Also in the last line there is a /@ to the right of it. How are these symbols working? –  Jon Oct 25 '13 at 6:03
    
Also it will not let me add a like to your comment, I have to have a 15 reputation in order for the web site to allow this. –  Jon Oct 25 '13 at 6:04
    
@Jon see this: mathematica.stackexchange.com/questions/18393/… –  RunnyKine Oct 25 '13 at 6:26

Just for a minor alternative:

arr = RandomReal[1, {100000, 2}];
g = GatherBy[arr, Norm@# < 1 &];
pi = 4. #1/(#1 + #2) & @@ (Length /@ g)
ListPlot[g, PlotStyle -> {Red, Green}, AspectRatio->1,Frame->True![\]][1]

wherepi is $\pi$ approximation

enter image description here

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.