Sample Code for Monte Carlo estimation of pi

Click For Summary
SUMMARY

The discussion focuses on using Monte Carlo sampling to estimate the value of pi through rejection sampling in MATLAB. The code provided demonstrates how to visualize the accepted and rejected samples within a unit circle inscribed in a square. The key point clarified is that the tilde (~) operator in MATLAB is used to create a logical negation, allowing the code to filter samples based on the rejection criteria defined by the vector 'reject'. The final estimate of pi is calculated as four times the mean of the accepted samples.

PREREQUISITES
  • Understanding of MATLAB programming and syntax
  • Familiarity with Monte Carlo methods for statistical estimation
  • Knowledge of logical indexing in MATLAB
  • Basic concepts of geometric probability and area estimation
NEXT STEPS
  • Explore MATLAB's logical indexing techniques in depth
  • Learn about advanced Monte Carlo methods and their applications
  • Investigate the implications of random number generation in MATLAB
  • Study the mathematical foundations of rejection sampling
USEFUL FOR

Data scientists, statisticians, and MATLAB users interested in numerical methods and probabilistic modeling, particularly those focusing on Monte Carlo simulations.

FallenApple
Messages
564
Reaction score
61
Below is the partial code to plot the rejection/accept regions using monte carlo sampling. The context is to estimate pi and plotting is just one part of it.

I'm not sure how the ~ symbol can work.

Code:
scatter(samples(1,~reject),samples(2,~reject),'b.')
scatter(samples(1,reject),samples(2,reject),'rx')

I get the idea that ~ means "not". But reject is a different vector. How can the samples matrix detect logicals on vectors that aren't even subsets of the matrix?

The entire code and website I got it from is below.
Code:
% DISPLAY A CIRCLE INSCRIBED IN A SQUARE

figure;
a = 0:.01:2*pi;
x = cos(a); y = sin(a);
hold on
plot(x,y,'k','Linewidth',2)

t = text(0.5, 0.05,'r');
l = line([0 1],[0 0],'Linewidth',2);
axis equal
box on
xlim([-1 1])
ylim([-1 1])
title('Unit Circle Inscribed in a Square')

pause;
rand('seed',12345)
randn('seed',12345)
delete(l); delete(t);

% DRAW SAMPLES FROM PROPOSAL DISTRIBUTION
samples = 2*rand(2,100000) - 1;

% REJECTION
reject = sum(samples.^2) > 1;

% DISPLAY REJECTION CRITERION
scatter(samples(1,~reject),samples(2,~reject),'b.')
scatter(samples(1,reject),samples(2,reject),'rx')
hold off
xlim([-1 1])
ylim([-1 1])

piHat = mean(sum(samples.*samples)<1)*4;

title(['Estimate of \pi = ',num2str(piHat)]);

https://theclevermachine.wordpress.com/2012/09/10/rejection-sampling/
 
Physics news on Phys.org
Orodruin said:
I suggest reading this: https://se.mathworks.com/help/matlab/math/matrix-indexing.html

If it is a tl;dr - scroll down directly to using logical arrays to access matrix elements.

Ok I see. How silly of me.

reject = sum(samples.^2) > 1;
is just a vector of logicals. True or False. So it can be used for conditioning.
 

Similar threads

  • · Replies 8 ·
Replies
8
Views
3K
  • · Replies 2 ·
Replies
2
Views
2K
  • · Replies 4 ·
Replies
4
Views
3K
  • · Replies 10 ·
Replies
10
Views
4K
  • · Replies 2 ·
Replies
2
Views
2K
  • · Replies 1 ·
Replies
1
Views
2K
  • · Replies 12 ·
Replies
12
Views
3K
Replies
1
Views
1K
  • · Replies 1 ·
Replies
1
Views
1K
  • · Replies 10 ·
Replies
10
Views
3K