Thread Subject:

Compute probability mass function on a predefined grid without loops

Hi, I have been trying for a while and with no success to solve the following simple problem without using "for" loops.
I have a predefined equally-spaced 1D grid of values , say
a=linspace(a_min,a_max,n);
I also have a vector of values s of the same size n with
max(s)<a_max and min(s)>a_min;
What I want to obtain is a probability mass function f(a) such that density associated to each s(i) is distributed to the grid values of a bracketing s(i) according to the relative distance of s(i) from such values.
In other words, suppose s(i) is bracketed by a(j) and a(j+1), then
f(a(j+1))=(s(i)-a(j))/(a(j+1)-a(j))
f(a(j))=1-f(a(j+1)).
I clearly was able to compute f(a) using for loops, but since both na and ns are pretty large (between 5000 and 10000) and since this computation is inside an fsolve I am looking for a smart way to vectorize it.
Thank you,
Massimo

"Massimo" <giovanma@bc.edu> wrote in message <j2h9i1$563$1@newscl01ah.mathworks.com>...
>
> What I want to obtain is a probability mass function f(a) such that density associated to each s(i) is distributed to the grid values of a bracketing s(i) according to the relative distance of s(i) from such values.
> In other words, suppose s(i) is bracketed by a(j) and a(j+1), then
> f(a(j+1))=(s(i)-a(j))/(a(j+1)-a(j))
> f(a(j))=1-f(a(j+1)).
===============

What happens if more than one
s(i) is bracketed by a(j) and a(j+1)?

"Matt J" wrote in message <j2i4ek$lb1$1@newscl01ah.mathworks.com>...

> What happens if more than one
> s(i) is bracketed by a(j) and a(j+1)?

Ok suppose both s(i) and s(i+1) lie between a(j) and a(j+1).
Then s(i) will be distributed to a(j) according to a weight w(i) [which I described above] that depends on the distance of s(i) from a(j) , while 1-w(i) will be distributed to a(j+1).
Similarly s(i+1) will be distributed to a(j) and a(j+1) according to a weight w(i+1) which measures the distance of s(i+1) from a(j). And so on.
At the end the density associated to a(j) will be just the sum of w(i) and w(i+1), similarly for the density at a(j+1).

I hope this helps.

Massimo

Here's one way, which uses this FEX tool

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



a=a(:); s=s(:); %ensure a and s are column vectors

idx=find_idx(s,a);

jlower=floor(idx);
jupper=jlower+1;

wtupper=(s-a(jlower))./(a(jupper)-a(jlower));
wtlower=1-wtupper;


f=accumarray([jlower;jupper],[wtlower;wtupper]),

"Matt J" wrote in message <j2k769$ker$1@newscl01ah.mathworks.com>...
> Here's one way, which uses this FEX tool
>
> http://www.mathworks.com/matlabcentral/fileexchange/23049-findidx
>
>
>
> a=a(:); s=s(:); %ensure a and s are column vectors
>
> idx=find_idx(s,a);
>
> jlower=floor(idx);
> jupper=jlower+1;
>
> wtupper=(s-a(jlower))./(a(jupper)-a(jlower));
> wtlower=1-wtupper;
>
>
> f=accumarray([jlower;jupper],[wtlower;wtupper]),

Matt, thanks. We are almost there...one question:
if I use find_idx and then apply floor(idx), couldn't I just use bin from
 [n bin] = histc(s,a)?
Massimo

No need to compute wtlower and wtupper from s, it can be computed directly from idx.

idx=find_idx(s,a);
jlower=floor(idx);
jupper=jlower+1;
wtlower = jupper-idx;
wtupper= 1-wtlower;

% Bruno

"Massimo" <giovanma@bc.edu> wrote in message <j2mc5l$kmq$1@newscl01ah.mathworks.com>...
>

>
> Matt, thanks. We are almost there...one question:
> if I use find_idx and then apply floor(idx), couldn't I just use bin from
> [n bin] = histc(s,a)?

I'm afraid don't see how that would apply, or why it would be desirable. The operations following floor(idx) are all much more trivial than HISTC in terms of computational expense

"Matt J" wrote in message <j2mi4e$bhm$1@newscl01ah.mathworks.com>...
>

> I'm afraid don't see how that would apply, or why it would be desirable. The operations following floor(idx) are all much more trivial than HISTC in terms of computational expense....
===============

....particularly with Bruno's proposed modification.