wolffd@0: function Y = sample_cond_multinomial(X, M)
wolffd@0: % SAMPLE_MULTINOMIAL Sample Y(i) ~ M(X(i), :)
wolffd@0: % function Y = sample_multinomial(X, M)
wolffd@0: %
wolffd@0: % X(i) = i'th sample
wolffd@0: % M(i,j) = P(Y=j | X=i) = noisy channel model
wolffd@0: %
wolffd@0: % e.g., if X is a binary image,
wolffd@0: % Y = sample_multinomial(softeye(2, 0.9), X)
wolffd@0: % will create a noisy version of X, where bits are flipped with probability 0.1
wolffd@0: 
wolffd@0: if any(X(:)==0)
wolffd@0:   error('data must only contain positive integers')
wolffd@0: end
wolffd@0: 
wolffd@0: Y = zeros(size(X));
wolffd@0: for i=min(X(:)):max(X(:))
wolffd@0:   ndx = find(X==i);
wolffd@0:   Y(ndx) = sample_discrete(M(i,:), length(ndx), 1);
wolffd@0: end
wolffd@0: 
wolffd@0: