Path: news.mathworks.com!not-for-mail
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$1@fred.mathworks.com>
References: <hokljp$c4f$1@fred.mathworks.com> <holf1i$92h$1@fred.mathworks.com> <hp1e6n$6hj$1@fred.mathworks.com> <hp65fe$d88$1@fred.mathworks.com> <hpbf5d$f77$1@fred.mathworks.com> <ik8llt$bnu$1@fred.mathworks.com>
Reply-To: <HIDDEN>
NNTP-Posting-Host: webapp-02-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1298663285 8410 172.30.248.37 (25 Feb 2011 19:48:05 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Fri, 25 Feb 2011 19:48:05 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: news.mathworks.com comp.soft-sys.matlab:712443

"Stefan Depeweg" wrote in message <ik8llt$bnu$1@fred.mathworks.com>...
> > "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hp65fe$d88$1@fred.mathworks.com>...
> > >  .......
> > >  [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