Random simulation (and stimulation?) Kiyoshi Akima Senior Member Posts: 325 Threads: 18 Joined: Jul 2006 07-17-2006, 02:23 PM Most modern calculators have random number generators (RNGs) built in, and they are useful things. Of course, most calculator (and computer) RNGs aren't truly random, though there are hardware RNGs (usually with a Geiger counter or other external input). Even software pseudo-RNGs can be useful. And not just for games. Much of probability theory is built upon the study of random events. Consider a large flat surface ruled with parallel lines d apart. Onto this drop a needle of length LRUN # Points Integral -------------------------------- 1 10 .186416796277 2 100 .206225523207 3 1000 .206581885362 4 10000 .20584278312 5 100000 .205729644257 6 1000000 .205804443536 7 10000000 .20580445907 ``` so your 8-dimensional integral evaluates to I = 0.20580, correct to all five digits plus or minus 1 ulp. If you alwo want to compute a running estimate of the error, this slightly longer version will do, still at 3 lines: ``` 1 DESTROY ALL @ RANDOMIZE PI @ FOR K=1 TO 7 @ M=10^K @ S=0 @ E=0 2 FOR N=1 TO M @ T=1/(1+RND+RND+RND+RND+RND+RND+RND+RND) @ S=S+T @ E=E+T*T 3 NEXT N @ I=S/M @ E=SQR((E/M-I*I)/M) @ DISP K;M,I,E @ NEXT K >RUN # Points Integral Estimated error -------------------------------------------------------- 1 10 .186416796277 7.87240479015E-3 2 100 .206225523207 3.63497557268E-3 3 1000 .206581885362 1.1943736672E-3 4 10000 .20584278312 3.62667298278E-4 5 100000 .205729644257 1.15685295773E-4 6 1000000 .205804443536 3.66089195702E-5 7 10000000 .20580445907 1.15814079667E-5 ``` so we can now confidently state that I = 0.205804 +- 0.000012 which arguably is "a little more precise than 1/9 < I < 1", as you requested. As is typical of Monte Carlo methods, the error is proportional to 1/Sqrt(N), where N is the number of random multidimensional points used, so the error diminishes by a factor of Sqrt(10) = 3.16+ for each increase of an order of magnitude in the number of points used. This is evident in the above data, where each Estimated Error approximatey equals the previous one divided by 3.16+. Faster convergence (order 1/N and better) may be easily obtained by a simple modification to the above 3-liners but for the purposes of computing "a little more precise" value, this will do. Your first question is a simple case of B*****'s needle, which I'll leave for others to deal with, and as for the World Cup I couldn't care less ! :-) Best regards from V. Edited: 18 July 2006, 2:21 p.m. RDunn Junior Member Posts: 3 Threads: 0 Joined: Jan 1970 07-20-2006, 02:04 PM For ten tosses of Buffon's Needle, the most likely outcome as given by the binomial formula is 7, which is a little surprising as 6 gives the best approximation to the expected probability, which is 2/pi. Here's a first stab at Buffon's Needle for the HP 48GX: ```<< 1 0 { 0 0 0 0 0 0 0 0 0 0 0 } -> n h p o << i n START 1 10 START IF RAND RAND pi * SIN < THEN 'h' INCR DROP END NEXT 'o' h GET 1 + 'o' h ROT PUT h 'p' STO+ 1 'h' STO NEXT o p n - n 10 * / >> >> ``` The program produces a list of the number of occurences of each of the eleven possible outcomes and the proportion of hits. Here's the output of a run of 100 trials, which took about 45 seconds: ```{ 0 0 0 2 8 14 24 28 14 7 3 } .653 ``` which gives 3.0627871363 as an estimate for pi. Kiyoshi Akima Senior Member Posts: 325 Threads: 18 Joined: Jul 2006 07-24-2006, 12:27 PM As has been mentioned, this is simply Buffon's Needle. The probability of a needle of length L intersecting one of a set of lines spaced d is 2L/(pi d). A fuller discussion, including programs in RPN, BASIC, and RPL, is available here in PDF format. Kiyoshi Akima Senior Member Posts: 325 Threads: 18 Joined: Jul 2006 07-24-2006, 12:54 PM Yes, the Monte Carlo is the obvious approach. But it's not necessarily the best. How long did it take for 1E7 points? Emu71 took about an hour on my machine and I estimate one to two weeks on a real 71B. It's true that the error for the simple Monte Carlo is proportional to 1/sqrt(N), and there are ways to improve this to 1/N. However, there is a way to evaluate this integral with an error proportional to 1/N2. And this scheme does at least as well as a 1E8-point Monte Carlo with N<1000. This gives a calculator or a real-speed emulator a chance to get five or six digits in time measurable by a clock instead of a calendar. BTW, what you call the "error estimate" is actually a one-standard-deviation bound on the error, not a rigorous bound. And since there's no guarantee that the errors follow a Gaussian distribution, it can only be used as a rough indication of probable error. Rodger Rosenbaum Senior Member Posts: 305 Threads: 17 Joined: Jun 2007 07-28-2006, 01:42 AM It has been more than a week now, and I'm eagerly awaiting your solution, especially if you're going to expound on: "It's true that the error for the simple Monte Carlo is proportional to 1/sqrt(N), and there are ways to improve this to 1/N. However, there is a way to evaluate this integral with an error proportional to 1/N2. And this scheme does at least as well as a 1E8-point Monte Carlo with N<1000. This gives a calculator or a real-speed emulator a chance to get five or six digits in time measurable by a clock instead of a calendar." Kiyoshi Akima Senior Member Posts: 325 Threads: 18 Joined: Jul 2006 07-28-2006, 12:28 PM As stated originally, the multi-dimensional integral will be covered in part 3. Due to the less than overwhelming amount of interest out there, I'll be taking it slowly. Watch for part 2 next week. But I did promise an answer to the Extra Extra Credit question tying the problem into the recent World Cup. As RDunn pointed out, Comte du Buffon was French. In addition, the goalkeeper for the winning Italian side was named Buffon. Rodger Rosenbaum Senior Member Posts: 305 Threads: 17 Joined: Jun 2007 07-28-2006, 07:08 PM The "extra credit" integral can be integrated without any special tricks other than lots of algebra, to give a final result of: (11760*Log[27] - 44914*Log[2] - 2352*Log[729/256] - 290829*Log[3] + 3525221*Log[4] + 451150*Log[5] - 15676416*Log[6] + 3136*Log[27/4] + 23059204*Log[7] - 16777216*Log[8] + 4782969*Log[9] - 107163*Log[256/27] + 1003520*Log[3125/256] + 7840*Log[16])/5040 which evaluates to 0.2057972212989439 (if I haven't made a typo). I don't think this is the technique you had in mind, but it does provide a benchmark to compare the accuracy of the Montecarlo results. I am aware of the so-called quasi-Montecarlo technique which can give convergence of roughly 1/N, but I'll have to wait to see the trick for 1/N^2 « Next Oldest | Next Newest »

 Possibly Related Threads... Thread Author Replies Views Last Post HP-70 Simulation for iPhone Willy R. Kunz 2 720 11-17-2013, 07:12 AM Last Post: Namir HP Prime: RANDOM Alberto Candel 4 873 10-18-2013, 09:18 PM Last Post: Alberto Candel Question about Simulation Namir 7 1,145 09-12-2013, 08:02 PM Last Post: Namir Some random wish list items for the 43s Marcel Samek 18 2,256 07-11-2013, 05:35 PM Last Post: Paul Dale Virtual HP-IL 40 col. video interface simulation Christoph Giesselink 10 1,532 08-19-2012, 06:46 PM Last Post: Richard Wagoner USB-41 82143A simulation update: MAN-TRACE modes Diego Diaz 0 464 06-09-2012, 08:46 PM Last Post: Diego Diaz. HP-41 > USB printer simulation (82143A) Diego Diaz 25 2,967 02-17-2012, 06:38 PM Last Post: Diego Diaz Emulator updates & HP82240B printer simulation Christoph Giesselink 14 1,819 10-28-2011, 05:26 AM Last Post: Nick_S OT: Help with Testing Random Number Generators Namir 5 902 08-01-2011, 11:20 PM Last Post: Paul Dale Question about random number generators Namir 20 2,060 07-08-2011, 11:05 AM Last Post: David Hayden

Forum Jump: