Andrew Cooke: C[omp]ute

Uniformly Random, Correlated Numbers in Matlab/Octave

From: andrew cooke <andrew@...>

Date: Sat, 23 Jul 2011 11:27:19 -0400

Someone just asked me how to implement the work at
http://acooke.org/random.pdf in Matlab.  Here's the code for Octave, which is
very similar:


n = 1000;  # the number of points
t = 0.1; # the 't' parameter from the paper

sq2 = sqrt(2);

# first, generate points in the u,v coordinates as in figure 4.
# u goes from 0 to sq2, v from -t*sq2 to t*sq2
uv = rand(n,2) * [sq2,0;0,2*t*sq2] - repmat([0,t*sq2],n,1);
# now rotate that by 45 degrees to xy
xy = uv * [1/sq2,1/sq2;-1/sq2,1/sq2];
# finally, "fold in" points outside the unit square
for i = 1:n
  #printf("%d x %f y %f\n", i, xy(i,1), xy(i,2));
  if (xy(i,1) < 0) xy(i,1) = -xy(i,1); endif
  if (xy(i,2) < 0) xy(i,2) = -xy(i,2); endif
  if (xy(i,1) > 1) xy(i,1) = 2-xy(i,1); endif
  if (xy(i,2) > 1) xy(i,2) = 2-xy(i,2); endif
endfor;

#xy
plot(xy(:,1),xy(:,2),'*')


If anyone with more experience than me know of a better way of folding the
points (ie without the explicit "for" loop), please let me know!

Andrew

Avoiding For Loops

From: Jonathan Dursi <ljdursi@...>

Date: Sat, 23 Jul 2011 12:26:57 -0400

I think this is pretty easily done with the vector-friendly min & max 
operations:

function xy = randminmax(n)
     t = 0.1; # the 't' parameter from the paper

     sq2 = sqrt(2);

     # first, generate points in the u,v coordinates as in figure 4.
     # u goes from 0 to sq2, v from -t*sq2 to t*sq2
     uv = rand(n,2) * [sq2,0;0,2*t*sq2] - repmat([0,t*sq2],n,1);
     # now rotate that by 45 degrees to xy
     xy = uv * [1/sq2,1/sq2;-1/sq2,1/sq2];

     # finally, "fold in" points outside the unit square
     xy = max(xy,-xy);
     xy = min(xy,2-xy);
endfunction


t1=now(); xy=randfor(100000); t2=now(); t2-t1
ans =  1.8261e-05
t1=now(); xy2=randminmax(100000); t2=now(); t2-t1
ans =  1.4284e-07

-- 
Jonathan Dursi | SciNet, Compute/Calcul Canada

Re: Avoiding For Loops

From: andrew cooke <andrew@...>

Date: Sat, 23 Jul 2011 13:50:57 -0400

Thanks!  I was looking for something like a map - never ocurred to me to look
for specific functions.  I'll forward that to the person who asked.

(And even when I read your suggestion, before I looked at the code, I couldn't
see how the upper bound would work.  I'm clearly not getting the matrix-based
approach at all... guess I better practice!)

Andrew

Comment on this post

C[omp]ute

Personal Projects

Last 100 entries

Uniformly Random, Correlated Numbers in Matlab/Octave

Avoiding For Loops

Re: Avoiding For Loops