Search code examples
matlaboctavemontecarlo

How to plot a point in some color based on the result of a comparaison, not using a loop?

I'm using random points to determine the area below a curve (Monte-Carlo):

  • X: 1xn vector of x values for the function
  • Y: 1xn vector of y = f(x)
  • RP: mxn matrix of m random y for each x

I would like to split RP into RPA and RPB depending on it being above or below the curve. The idea is then to plot RPA and RPB against X, in different colors. This code doesn't work because RPA and RPB number of columns is not the same than X:

clf
f = @(x) sin(x/10) + cos(x/60); % Function
xMin = 1; xMax = 100; % x interval

X = [xMin:xMax];
Y = f(X);
plot(X,Y), hold on % Plot function

yMin = min(Y); yMax = max(Y); % Axes limits
set(gca, 'xlim', [xMin, xMax], 'ylim', [yMin, yMax])

m = 20; % Random points per x value
RP = rand(m, columns(X)) .* (yMax-yMin) .+ yMin;

% Split points (doesn't work)
RPA = RP(RP>Y);
RPB = RP(RP<=Y);

br = size(RPB) / size(RP) % Ratio of points below
a = (xMax - xMin) * (yMax - yMin) * br % Area below

% Plot points
plot(X, RPA, 'r.') % Above
plot(X, RPB, 'g.') % Below

Is there a possibility to build RPA and RPB so that they are the same size than RP, with the excluded points y being NaN or something similar, which can be counted and plotted?

Solution

  • You gave a good answer yourself. You can build RPA and RPB with strategic NaNs:

    % Split points (works!)
RPA = RP;
RPA(RP<=Y) = NaN;
RPB = RP;
RPB(RPB > Y) = NaN;

    And than calculating the ration as the not-NaN:

    br = sum(~isnan(RPB)) / sum(~isnan(RP)) % Ratio of points below

    I get this nice image: enter image description here