Convert a bivariate draw in a univariate draw in Matlab

Question

I have in mind the following experiment to run in Matlab and I am asking for an help to implement step (3). Any suggestion would be very appreciated.

(1) Consider the random variables X and Y both uniformly distributed on [0,1]

(2) Draw N realisation from the joint distribution of X and Y assuming that X and Y are independent (meaning that X and Y are uniformly jointly distributed on [0,1]x[0,1]). Each draw will be in [0,1]x[0,1].

(3) Transform each draw in [0,1]x[0,1] in a draw in [0,1] using the Hilbert space filling curve: under the Hilbert curve mapping, the draw in [0,1]x[0,1] should be the image of one (or more because of surjectivity) point(s) in [0,1]. I want pick one of these points. Is there any pre-built package in Matlab doing this?

I found this answer which I don't think does what I want as it explains how to obtain the Hilbert value of the draw (curve length from the start of curve to the picked point)

On wikipedia I found this code in C language (from (x,y) to d) which, again, does not fulfil my question.

Answer 1

I will focus only on your last point

(3) Transform each draw in [0,1]x[0,1] in a draw in [0,1] using the Hilbert space filling curve: under the Hilbert curve mapping, the draw in [0,1]x[0,1] should be the image of one (or more because of surjectivity) point(s) in [0,1]. I want pick one of these points. Is there any pre-built package in Matlab doing this?

As far as I know, there aren't pre-built packages in Matlab doing this, but the good news is that the code on wikipedia can be called from MATLAB, and it is as simple as putting together the conversion routine with a gateway function in a xy2d.c file:

#include "mex.h"

// source: https://en.wikipedia.org/wiki/Hilbert_curve
// rotate/flip a quadrant appropriately
void rot(int n, int *x, int *y, int rx, int ry) {
    if (ry == 0) {
        if (rx == 1) {
            *x = n-1 - *x;
            *y = n-1 - *y;
        }

        //Swap x and y
        int t  = *x;
        *x = *y;
        *y = t;
    }
}

// convert (x,y) to d
int xy2d (int n, int x, int y) {
    int rx, ry, s, d=0;
    for (s=n/2; s>0; s/=2) {
        rx = (x & s) > 0;
        ry = (y & s) > 0;
        d += s * s * ((3 * rx) ^ ry);
        rot(s, &x, &y, rx, ry);
    }
    return d;
}


/* The gateway function */
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[])
{
    int n;              /* input scalar */
    int x;              /* input scalar */
    int y;              /* input scalar */
    int *d;             /* output scalar */

    /* check for proper number of arguments */
    if(nrhs!=3) {
        mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nrhs","Three inputs required.");
    }
    if(nlhs!=1) {
        mexErrMsgIdAndTxt("MyToolbox:arrayProduct:nlhs","One output required.");
    }

    /* get the value of the scalar inputs  */
    n = mxGetScalar(prhs[0]);
    x = mxGetScalar(prhs[1]);
    y = mxGetScalar(prhs[2]);

    /* create the output */
    plhs[0] = mxCreateDoubleScalar(xy2d(n,x,y));

    /* get a pointer to the output scalar */
    d = mxGetPr(plhs[0]);
}

and compile it with mex('xy2d.c').

The above implementation

[...] assumes a square divided into n by n cells, for n a power of 2, with integer coordinates, with (0,0) in the lower left corner, (n-1,n-1) in the upper right corner.

In practice, a discretization step is required before applying the mapping. As in every discretization problem, it is crucial to choose the precision wisely. The snippet below puts everything together.

close all; clear; clc;

% number of random samples
NSAMPL = 100;

% unit square divided into n-by-n cells
% has to be a power of 2
n = 2^2;

% quantum
d = 1/n;

N = 0:d:1;

% generate random samples
x = rand(1,NSAMPL);
y = rand(1,NSAMPL);

% discretization
bX = floor(x/d);
bY = floor(y/d);

% 2d to 1d mapping
dd = zeros(1,NSAMPL);
for iid = 1:length(dd)
    dd(iid) = xy2d(n, bX(iid), bY(iid));
end


figure;
hold on;
axis equal;

plot(x, y, '.');
plot(repmat([0;1], 1, length(N)), repmat(N, 2, 1), '-r');
plot(repmat(N, 2, 1), repmat([0;1], 1, length(N)), '-r');


figure;
plot(1:NSAMPL, dd);
xlabel('# of sample')