Line | |
---|
1 | function h = sample_hist(p, n) |
---|
2 | %SAMPLE_HIST Sample from a multinomial distribution. |
---|
3 | % SAMPLE_HIST(P,N) returns a random matrix of counts, same size as P, |
---|
4 | % whose column sums are all N. Column j is sampled from a multinomial with |
---|
5 | % the probabilities p(:,j). |
---|
6 | % |
---|
7 | % Example: |
---|
8 | % sample_hist([0.2 0.4; 0.8 0.6],100) |
---|
9 | |
---|
10 | % The advantage of this alg is that the running time grows slowly with n. |
---|
11 | % It is the same alg used by BUGS. |
---|
12 | % If n is very small then it is faster to just take n samples from p. |
---|
13 | |
---|
14 | if nargin < 2 |
---|
15 | n = 1; |
---|
16 | end |
---|
17 | |
---|
18 | h = zeros(size(p)); |
---|
19 | z = repmat(1,1,cols(p)); |
---|
20 | n = repmat(n,1,cols(p)); |
---|
21 | js = 1:cols(p); |
---|
22 | % loop bins |
---|
23 | for i = 1:(rows(p)-1) |
---|
24 | % the count in bin i is a binomial distribution |
---|
25 | for j = js |
---|
26 | h(i,j) = randbinom(p(i,j)/z(j), n(j)); |
---|
27 | end |
---|
28 | n = n - h(i,:); |
---|
29 | z(js) = z(js) - p(i,js); |
---|
30 | js = find(z > 0); |
---|
31 | end |
---|
32 | h(end,:) = n; |
---|
Note: See
TracBrowser
for help on using the repository browser.