user2005253
user2005253

Reputation:

Sampling in Matlab

So let me start off by saying that I do not have the statistics toolbox for Matlab so I am trying to find a way to work around this. In any case, what I am trying to do is to replicate the R sample function. For example, in R

> x = sample(1:5,20,replace=T,prob=c(.1,.1,.1,.1,.6))
> x
 [1] 5 5 5 4 5 2 5 5 1 5 5 5 5 5 5 3 5 1 5 5

so I am sampling the integers 1,2,3,4,5 with replacement. But furthermore, I am sampling each integer with a certain proportion, i.e., the integer 5 should be sampled about 60% of the time.

So my question that I would like to find a solution to is how to achieve this in Matlab?

Upvotes: 1

Views: 1476

Answers (3)

shoelzer
shoelzer

Reputation: 10708

Here's one way to make your own sample function:

function x = sample(v, n, p)

pc = cumsum(p) / sum(p);
r = rand(1,n);
x = zeros(1,n);
for i = length(pc):-1:1
    x(r<pc(i)) = v(i);
end

It's not exactly efficient, but it does what you want. Call it like so:

v = [1 2 3 4 5];
p = [.1 .1 .1 .1 .6];
n = 20;
x = sample(v,n,p);

Upvotes: 1

Jonas
Jonas

Reputation: 74940

Here's how you can perform weighted sampling with replacement (something Matlab's randsample doesn't support, btw);

function r = sample(pop,n,weights)
%# each weight creates a "bin" of defined size. If the value of a random number
%# falls into the bin, we pick the value

%# turn weights into a normed cumulative sum
csWeights = cumsum(weights(:))/sum(weights);
csWeights = [0;csWeights(1:end-1)];

%# for each value: pick a random number, check against weights
idx = sum(bsxfun(@ge,rand(1,n),csWeights),1);

r = pop(idx);

Upvotes: 4

Richie Cotton
Richie Cotton

Reputation: 121077

The unweighted case is easy using randi.

function r = sample(pop, n)

imax = length(pop);
index = randi(imax, n, 1);
r = pop(index);

end

In the weighted case, something like this should do the trick:

function r = sample(pop, n, prob)

cumprob = cumsum(prob);
r = zeros(1, n);
for i = 1:n
  index = find(rand < cumprob, 1, 'last');
  r(i) = pop(index);
end

end

Upvotes: 3

Related Questions