From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Compute a bivariate probability mass function on a predefined grid without loops
Date: Sun, 22 Apr 2012 17:12:09 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 32
Message-ID: <jn1e59$oi9$>
References: <jms4dg$7db$> <jmtam1$e0q$> <jmu49f$j4h$> <jmu8rb$6fh$> <jmuhhc$97p$> <jmujvn$i1b$> <jmum1d$pfn$> <jn0gib$4dv$>
Reply-To: <HIDDEN>
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: 1335114729 25161 (22 Apr 2012 17:12:09 GMT)
NNTP-Posting-Date: Sun, 22 Apr 2012 17:12:09 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1187260
Xref: comp.soft-sys.matlab:765522

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <jn0gib$4dv$>...
> .......
> The question I wonder is : does the direct calculation 
> p = (Sk-k(1))*((N-1)/(k(N)-k(1)))+1
> always lead to the correct index from an input vector generated by LINSPACE due to numerical accuracy.
> We have discussed in the past the discrepancy between LINSPACE and COLON-generated vectors. Similarly I guess in some vary extreme cases, p could return a wrong interval (?), whereas dichotomy search will always return the "right" interval, thus more robust. Your opinion?
> .......
"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <jn0iar$ar0$>...
> ......
> p = ( ( (Sk-k(1)) * (N-1) ) ./ (k(N)-k(1)) ) + 1;
> should undo accurately what LINSPACE did.
- - - - - - - - - - -
  Hi again, Bruno.  My reason for suggesting the formulas

 p = (Sk-k(1))*((N-1)/(k(N)-k(1)))+1;
 q = (Sm-m(1))*((N-1)/(m(N)-m(1)))+1;

was based, as I have said, on the supposition that 'interp1' would have to use an order(N*log(N)) algorithm to accomplish the same task unless it first tested for uniformity of distribution of k and m, and would therefore take more time for large N.

  As to its accuracy, that is quite another matter.  Certainly there is the possibility that Sk or Sm values that happen to be exceedingly close to k or m values will end up in the "wrong" rectangle.  What this really amounts to is a difference of opinion, so to speak, between these formulas and where 'linspace' in generating k and m would place the rectangle boundaries.  Both are subject to rounding errors.

  With the weighting functions that you and I have recommended in this thread, there is continuity as you cross these rectangular boundaries, so the consequences of being in the wrong rectangle are very small, of the same magnitude as rounding errors.  That is, as you cross a boundary, two of the weighting factors are simultaneously dropping down to zero while two others on an adjacent rectangle on the opposite side are starting up from zero.  In other words, I am suggesting that perhaps robustness in this respect isn't really necessary in this weighting problem.

  As to your addendum, yes, if you know exactly how 'linspace' does its computation, then using precisely the same formula with Sk and Sm with identical parentheses would (presumably) ensure that the correct rectangle is always selected.  It would follow from the IEEE 754 rounding specifications (assuming no weird actions by the compiler.)  That assumes of course that the user actually used 'linspace' to generate k and m rather than, say, 'colon'.

  I chose the other parentheses arrangement to ensure that the compiler only computes the factor (N-1)/(k(N)-k(1)) once, but a sufficiently smart compiler really ought to realize that anyway.

Roger Stafford