Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Matlab Sample Code for Monte Carlo estimation of pi

  1. Jul 19, 2017 #1
    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 (Text):

    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 (Text):



    % 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/
     
  2. jcsd
  3. Jul 19, 2017 #2

    Orodruin

    User Avatar
    Staff Emeritus
    Science Advisor
    Homework Helper
    Gold Member

  4. Jul 19, 2017 #3
    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Sample Code for Monte Carlo estimation of pi
  1. Latex Code (Replies: 4)

  2. Matlab Code (Replies: 2)

Loading...