# C[omp]ute

Welcome to my blog, which was once a mailing list of the same name and is still generated by mail. Please reply via the "comment" links.

Always interested in offers/projects/new ideas. Eclectic experience in fields like: numerical computing; Python web; Java enterprise; functional languages; GPGPU; SQL databases; etc. Based in Santiago, Chile; telecommute worldwide. CV; email.

© 2006-2013 Andrew Cooke (site) / post authors (content).

## 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