Thread Subject:

speed up for loop inside for loop short question

Hi there,

Given three column vectors x,y,z I want to calculate the probabilities:

P(X<Y<Z), P(X=Y<Z), P(X<Y=Z), P(X=Y=Z).

The following code gets the job done but I would like it to be faster. Is there a way of getting rid of at least one of the for loops?

(I described the problem above instead of just giving the code because it might be helpful for someone to take a completely different approach, or check me (I think this is correct). Anyway if you are bored to check, any ideas regarding the speed would be more than welcome!).

Thanx in advance for any answers!

op=1;
for i=1:length(x)
    for j=1:length(y)
        xx=x(i).*ones(length(z),1);
        yy=y(j).*ones(length(z),1);
        k1(op)=length(find(xx<yy & yy<z));
        k2(op)=length(find(xx<yy & yy==z));
        k3(op)=length(find(xx==yy & yy<z));
        k4(op)=length(find(xx==yy & yy==z));
        op=op+1;
    end
end
op=op-1;
nx=length(x);ny=length(y);nz=length(z);
pith1=sum(k1)./(nx*ny*nz) %P(X<Y<Z)
pith2=sum(k2)./(nx*ny*nz) %P(X<Y=Z)
pith3=sum(k3)./(nx*ny*nz) %P(X=Y<Z)
pith4=sum(k4)./(nx*ny*nz) %P(X=Y=Z)

I think I found my way out!

nx=length(x);
ny=length(y);
nz=length(z);


X=reshape(repmat(x,[1 ny*nz])',1,numel(repmat(x,[1 ny*nz])))';
Y=reshape(repmat(y,[nx nz])',1,numel(repmat(y,[nx nz])))';
Z=repmat(z,[nx*ny,1]);

n=length(X);
pith1=length(find(X<Y & Y<Z ))./n;
pith2=length(find(X==Y & Y<Z ))./n;
pith3=length(find(X<Y & Y==Z))./n;
pith4=length(find(X==Y & Y==Z))./n;


Is the above correct?

"leo nidas" <bleonidas25@yahoo.gr> wrote in message <j11amh$o12$1@newscl02ah.mathworks.com>...
> Given three column vectors x,y,z I want to calculate the probabilities:
> P(X<Y<Z), P(X=Y<Z), P(X<Y=Z), P(X=Y=Z).
> ........
- - - - - - - - - - -
  I claim that with large nx, ny, and nz, if you initially sort your three vectors, the total number of required floating point operations can be dramatically reduced in spite of the larger number of lines of code needed. Try the following for large size vectors and see if it isn't faster.

 nx = length(x); ny = length(y); nz = length(z);
 X = [sort(x);inf]; Y = [sort(y);inf]; Z = [sort(z);inf];
 p1 = 0; p2 = 0; p3 = 0; p4 = 0;
 ix = 1; xd = 0; iz = 1; zg = nz;
 f = true;
 for iy = 1:ny
  Y0 = Y(iy);
  if f
   xe = 0;
   while X(ix) <= Y0
    xd = xd + 1;
    if X(ix) == Y0, xe = xe + 1; end
    ix = ix + 1;
   end
   ze = 0;
   while Z(iz) <= Y0
    zg = zg - 1;
    if Z(iz) == Y0, ze = ze + 1; end
    iz = iz + 1;
   end
  end
  p1 = p1 + (xd-xe)*zg; p2 = p2 + xe*zg;
  p3 = p3 + (xd-xe)*ze; p4 = p4 + xe*ze;
  f = Y0 < Y(iy+1);
 end
 n = nx*ny*nz;
 p1 = p1/n; p2 = p2/n; p3 = p3/n; p4 = p4/n;

  Note: My p1, p2, p3, and p4 are your pith1, pith2, pith3, and pith4, resp.

Roger Stafford

Dear leo nidas,

I think your and Roger's solution are fine already. But your loop method can be improved also: There is no need to expand xx and yy!

 n = length(x) * length(y);
 k1 = zeros(1, n); % Preallocate !!!
 k2 = zeros(1, n);
 k3 = zeros(1, n);
 k4 = zeros(1, n);
 op=1;
 for i=1:length(x)
     for j=1:length(y)
         xx=x(i);
         yy=y(j);
         k1(op)=sum(xx<yy & yy<z);
         k2(op)=sum(xx<yy & yy==z);
         k3(op)=sum(xx==yy & yy<z);
         k4(op)=sum(xx==yy & yy==z);
         op=op+1;
     end
 end

Next improvement: if xx<yy if TRUE, the number of (xx==yy & yy==z) is zero:
 ...
 for i=1:length(x)
     for j=1:length(y)
         xx=x(i);
         yy=y(j);
         if xx < yy
           k1(op)=sum(yy<z);
           k2(op)=sum(yy==z);
        elseif xx==yy
           k3(op)=sum(yy<z);
           k4(op)=sum(yy==z);
         end
         op=op+1;
     end
 end

Kind regards, Jan

Here is another algorithm, I'll restrict of calculating P(X<=Y<=Z) = p1+p2+p3+p4

%
x=rand(4000,1);
y=rand(5000,1);
z=rand(6000,1);

% Engine
n = [length(x) length(y) length(z)];
[a is] = sort([x(:); y(:); z(:)]);
bins = 1+[0 cumsum(n)];
[trash, loc] = histc(is, bins);
p = [1 0 0 0];

for j = loc.'
    p(j+1) = p(j+1)+p(j);
end
p = p(4)/prod(n)

%%
I'll let the calculation of p1, p2, p3 ans p4 separately as exercise for whom who is still interested in the problem :-)

% Bruno