CODEWITHSUNDEEP
matlabrandomprobability-density
mononono
mononono

Reputation: 161

How to generate a 2D random vector in MATLAB?

I have a non-negative function f defined on a unit square S = [0,1] x [0,1] such that enter image description here

My question is, how can I use MATLAB to generate a 2D random vector from S according to the probability density function f?

Upvotes: 4

Views: 536

Answers (2)

Daniel
Daniel

Reputation: 36720

ICDF

Besides rejection sampling there is a different approach to address this issue on a more mathematical level, but you need to sit down and do some math first to end up with a better solution. For 1 dimensional distributions you typically sample using the ICDF (inverted cumulative density function) function simply using ICDF(rand(n,1)) to get random samples.

If you manage to do the math, you could instead for your PDF function define two functions ICDF1 (ICDF for the first dimension) and ICDF2 (ICDF for the second dimension) in matlab.

The first ICDF1 would map unifrom random distributed samples to sample values for the first dimension of your random distribution.

The second ICDF2 would map the output if ICDF1 and uniform distributed samples to your intended solution.

Here is some matlab code assuming you already defined ICDF1 and ICDF2

samples=ICDF1(rand(n,1));
samples(:,2)=ICDF2(samples,rand(n,1));

The great advantage of this solution is, that it does not reject any samples, being potentially much faster.

Upvotes: 1

Daniel
Daniel

Reputation: 36720

Rejection Sampling

The suggestion Luis Mendo made is very good because it applies to nearly all distribution functions. Based on this answer I wrote code for m.

An important point when using rejection sampling this way is that you must know the maximum of your pdf within the range. If you over-estimate the maximum your code will only run slower. If you under-estimate it it will create wrong numbers!

The idea is that you sample many uniform distributed points and accept depending on the probability density for the points.

pdf=@(x).5.*x(:,1)+3./2.*x(:,2);
maximum=2; %Right maximum for THIS EXAMPLE. 
%If you are unable to determine the maximum of your 
%function within the [0,1]x[0,1] range, please give an example.
result=[];
n=10;
while (size(result,1)<n)
    %1. sample random point:
    val=rand(1,2);
    %2. Accept with probability pdf(val)/maximum
    if rand<pdf(val)/maximum
        %append to solution
        result(end+1,:)=val;
    end
end

I know that this solution is not a fast implementation, but I wanted to start with an implementation as simple as possible to make sure that the concept of rejection sampling becomes clear.

Upvotes: 1

Related Questions