Rev | Line | |
---|
[37] | 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.