# Probability estimate using Monte Carlo methods

Jump to: Source code Text output Plots Library index

    function prob = montecarlo(A,b,Sigma,notrials);

    % Estimates the probability that a random vector x in R2
    % with mean zero and covariance Sigma satisfies Ax <= b,
    % based on 100 * notrials trials. Sigma must be PSD.

    m = size(A,1);

    R = chol(Sigma);   % Y = R^{-T}X has covariance I
    X = R'*randn(2,notrials);
    prob = length(find(sum(A*X - b(:,ones(1,notrials)) < 0) == m))/notrials;

    for i=1:99
    X = R'*randn(2,notrials);
    prob = 0.5*(prob + ...
      length(find(sum(A*X - b(:,ones(1,notrials)) < 0) == m))/notrials);
    end;