モンテカルロシミュレーションで円周率を求める(D次元対応版)

pi_n.sas

%macro pi(dim=, rep=);

    data pi;
        do i = 1 to &rep.;
            %do j = 1 %to &dim.;
                x&j. = rand("uniform") ** 2;
            %end;
            dist = sqrt(sum(of x:));
            if dist <= 1 then cnt + 1;
            calc_pi = ((cnt / i) * (2 ** &dim.) * gamma(&dim. / 2 + 1)) ** (2 / &dim.);
            output;
        end;
        put calc_pi = ;
    run;
    
    %let pi = %sysfunc(constant(pi));

    title "&dim.次元";
    proc sgplot data = pi;
        series x = i y = calc_pi;
        refline &pi. / axis = y;
        yaxis min = 0;
    run;
    title;

%mend pi;

%pi(dim=2,  rep=10000);
%pi(dim=3,  rep=10000);
%pi(dim=4,  rep=10000);
%pi(dim=5,  rep=10000);
%pi(dim=6,  rep=10000);
%pi(dim=7,  rep=10000);
%pi(dim=8,  rep=10000);
%pi(dim=9,  rep=10000);
%pi(dim=10, rep=10000);
%pi(dim=11, rep=10000);
%pi(dim=12, rep=10000);
%pi(dim=13, rep=10000);
%pi(dim=14, rep=10000);
%pi(dim=15, rep=10000);