From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Please help me to generate a markov chain.
Date: Fri, 25 Feb 2011 19:48:05 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 31
Message-ID: <ik911l$86q$>
References: <hokljp$c4f$> <holf1i$92h$> <hp1e6n$6hj$> <hp65fe$d88$> <hpbf5d$f77$> <ik8llt$bnu$>
Reply-To: <HIDDEN>
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: 1298663285 8410 (25 Feb 2011 19:48:05 GMT)
NNTP-Posting-Date: Fri, 25 Feb 2011 19:48:05 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: comp.soft-sys.matlab:712443

"Stefan Depeweg" wrote in message <ik8llt$bnu$>...
> > "Roger Stafford" <> wrote in message <hp65fe$d88$>...
> > >  .......
> > >  [t,s_n] = histc(rand,[-inf;cumsum(p(1:n-1));inf]);
> > >  .......
> In:  histc(rand,[-inf;cumsum(p(1:n-1));inf]);
> Isn't it  p(1:n)?   p has the dimension  nx1. If you just  call it with p(1:n-1) you cut off the last possible state for your calculation.
> Stefan
- - - - - - - - - - -
  In reply to your question, Stefan, where (almost a year ago) I said

 [t,s_n] = histc(rand,[-inf;cumsum(p(1:n-1));inf]); ,

in accordance with the definition of histc, the n-th, next-to-last of the n+1 "edges" intervals there would be:

 sum(p(1:n-1)) <= rand < inf.

Since the p-values are assumed to all be non-negative and to have a sum of one, this interval would actually be

 sum(p(1:n-1)) <= rand <= 1 = sum(p(1:n))

which has width of p(n), and is just what was intended.  None of the n possible values is left out.  (The last "edges" value of course would be an impossible "rand == inf" and would never occur.)

  The reason for not just using cumsum(p) was to guard against possible round-off errors in calculating 'cumsum'.  Any erroneous final value slightly below one would provide a small chance for an invalid bin value for s_n of n+1 which I wanted to avoid.

  (As I mentioned before, using histc on a single rand value may not be the most efficient way of determining which interval rand falls in.  A simple binary search procedure might be better.)

Roger Stafford