vectorizing challenge

[Heywood]

Hi all,

I'm having a hard time trying to vectorize the following code. Suppose I have a vector of interest V, and a reference vector R which is monotonically increasing (it's actually a cumulative distribution function). Both vectors contain only real numbers, except that I've modified R by tacking on -Inf at the beginning and +Inf at the end (the reason for which is hopefully obvious from the code below).

For each element of V, I want to find the first element of R that is greater than or equal to it.

To do this with a FOR loop is trivial (and yes, I realize the two lines in the loop can be combined into one; writing it this way just for clarity):

result = zeros(1,length(V));
for kk = 1:length(V);
ind = find( V(kk) <= R, 1, 'first' );
result(kk) = R(ind);
end;

But when V is very large, this is slow. Can anyone think of a clean way to 'vectorize' the FIND command?

Thanks very much in advance!

HJ

[Doug Schwarz]

In article <guvhl1$qif$1...@fred.mathworks.com>,
"Heywood " <heywo...@yahoo.com> wrote:

How about this?

% Create some sample data.
dist = rand(1,10);
R = [0,cumsum(dist)];
R = R/R(end);

% Find result.
n = length(R);
V = rand(1,100);
k = interp1(R,0:n,V);
result = R(ceil(k) + 1);

You are still finding the appropriate interval, but it's done
efficiently inside interp1.

--
Doug Schwarz
dmschwarz&ieee,org
Make obvious changes to get real email address.

[Matt]

"Heywood " <heywo...@yahoo.com> wrote in message <guvhl1$qif$1...@fred.mathworks.com>...

> Hi all,
>
> I'm having a hard time trying to vectorize the following code. Suppose I have a vector of interest V, and a reference vector R which is monotonically increasing (it's actually a cumulative distribution function). Both vectors contain only real numbers, except that I've modified R by tacking on -Inf at the beginning and +Inf at the end (the reason for which is hopefully obvious from the code below).
>
> For each element of V, I want to find the first element of R that is greater than or equal to it.

------

This one is probably not as efficient as Doug's, but is nice and compact

[dummy,inds] = max( bsxfun(@ge,R(:),V(:).') );
Result = R(inds);

[Bruno Luong]

"Heywood " <heywo...@yahoo.com> wrote in message <guvhl1$qif$1...@fred.mathworks.com>...

> result = zeros(1,length(V));

> for kk = 1:length(V);
> ind = find( V(kk) <= R, 1, 'first' );
> result(kk) = R(ind);
> end;

If speed is your main concern, don't use INTERP1 (slow), use HISTC or FIND_IDX on FEX

http://www.mathworks.com/matlabcentral/fileexchange/23049

Bruno

[Jos]

Doug Schwarz <s...@sig.for.address.edu> wrote in message <see-FBBAFE.2...@news.frontiernet.net>...

or

[dummy, k2] = histc(V,R)
result2 = R(k2+1) ;

Jos

[Heywood]

> [dummy, k2] = histc(V,R)
> result2 = R(k2+1) ;
>
> Jos

Hi Doug, Matt, Bruno, and Jos,

Thanks for all the suggestions! I learned a lot from all your ideas (I'd never seen bsxfun before, and only vaguely remembered interp1). Ultimately Jos's suggestion was the simplest and did the trick -- resulting in a ~20X speedup in execution time!

Thanks again everyone.

Cheers,

Heywood