your function seems to complicate the problem a little bit.
The following line is enough to do the jod
[~,x] = histc(rand(1,n),[0;cumsum(p(:))/sum(p)]);
Jos,
I have checked with randp. Though they seem offering similar functionalities, however, the efficiency is drastically different, especially in very large scale monte carlo simulation, say you need to draw thousands or millions of samples from a distribution over thousands or even millions of states, which is not unusual in real engineering applications.
With randp, it would be incur obvious latency when you want to draw thousands of samples from thousands of states, or even run out of memory (thus resulting an empty matrix), as the algorithm implemented by randp is of complexity O(k n), where k is the number of states in the sample space. However, even million-state-level sampling can be accomplished by this function within milliseconds, as its complexity is only O(n logk).
Hope my explanation can clarify the differences between these two files.
Comment only