我想到了在 Matlab 中运行的以下实验，我正在寻求帮助来实现步骤 (3)。任何建议将不胜感激。
(1) 考虑随机变量X和Y均均匀分布在[0,1]
(2)N从联合分布中提取实现，X假设Y和X是Y独立的（意味着X和Y均匀地联合分布在 上[0,1]x[0,1]）。每次抽奖将在[0,1]x[0,1].
(3)使用 Hilbert 空间填充曲线对 draw in[0,1]x[0,1]中的每个 draw in 进行变换[0,1]：在 Hilbert 曲线映射下，draw in[0,1]x[0,1]应该是 中的一个（或多个由于超射性）点的图像[0,1]。我想选择这些点之一。Matlab 中是否有任何预先构建的软件包可以做到这一点？
我找到了这个答案，我认为这不是我想要的，因为它解释了如何获得平局的希尔伯特值（从曲线起点到选取点的曲线长度）
在维基百科上，我在 C 语言中找到了这段代码（从(x,y)to d），这同样不能满足我的问题。
编辑 这个答案没有解决问题的更新版本，它明确询问构建希尔伯特曲线。相反，这个答案解决了一个关于双射映射构造的相关问题，以及与均匀分布的关系。
你的问题没有很好的定义。如果您只需要生成的分布是均匀的，那么没有什么能阻止您简单地选择f:(X,Y)->X. 无论X和Y是否相关，结果都是一致的。从您的帖子中，我只能假设您实际上想要的是使结果转换是双射的，或者在机器精度限制的情况下尽可能接近它。
值得注意的是，除非您需要最能保留局部性的算法（显然不需要结果分布是双射的，更不用说均匀了），否则无需费心构建您在问题中提到的希尔伯特曲线。它们与解决方案的关系与任何其他空间填充曲线一样多，并且计算量非常大。
因此，假设您正在寻找双射映射，您的问题相当于询问[unit] 正方形中的点集是否与 [unit] 线段中的点集具有相同的基数，如果是，如何构造该双射，即 1 对 1 对应。直觉说正方形应该有更高的基数，康托尔花了3 年时间试图证明这一点，最终证明完全相反——这些集合实际上是等数的。他对自己的发现感到非常惊讶，于是他写道：
我看到了，但我不相信！
满足**此标准的最常提到的双射如下。用它们的十进制形式表示x和y，即x=0。 x1 x2 x3 x4 x5...，并且y=0。 y1 y2 y3 y4 y5 ...，设z=0。x1 y1 x2 y2 x3 y3 x4 y4 x5 y5...，即交替两个数字的小数。双射背后的想法是微不足道的，尽管严格的证明需要相当多的先验知识。f:(X,Y)->Z
** 需要注意的是，如果我们采用 egx = 1/3 = 0.33333...和y = 1/5 = 0.199999... = 0.200000...，我们可以看到有两个序列对应于它们：z = 0.313939393939...和z = 0.323030303030...。为了克服这个障碍，我们必须证明将可数集添加到不可数集不会改变后者的基数。
实际上，我们必须处理机器精度而不是纯数学，这严格来说意味着这两个集合实际上都是有限的，因此不是等值的（假设您以与原始数字相同的精度存储结果）。x这意味着我们只是被迫做一些假设并丢失一些信息，例如，在这种情况下，和的有效数字的后半部分y。也就是说，除非我们使用与原始变量相比允许以双精度存储结果的不同数据类型。
最后，在 Matlab 中的示例实现：
x = rand();
y = rand();
chars = [num2str(x, '%.17f'); num2str(y, '%.17f')];
z = str2double(['0.' reshape(chars(:,3:end), 1, [])]);
>> cellstr(['x=' num2str(x, '%.17f'); 'y=' num2str(y, '%.17f'); 'z=' num2str(z, '%.17f')])
ans =
'x=0.65549803980353738'
'y=0.10975505072305158'
'z=0.61505947958500362'
编辑 这回答了在给定 x,y ~ U[0,1] 的情况下对转换 f(x,y) -> t ~ U[0,1] 的原始请求，另外还用于 x 和 y 相关。更新后的问题专门询问希尔伯特曲线 H(x,y) -> t ~ U[0,1] 并且仅适用于 x,y ~ U[0,1] 所以这个答案不再相关。
考虑 [0,1] r1, r2, r3, ...中的随机均匀序列。您正在将此序列分配给数字对 (x1,y1), (x2,y2), .... 您在问什么for 是对 (x,y) 对的变换，它在 [0,1] 中产生一个均匀随机数。
考虑随机子序列 r1, r3, ... 对应于 x1, x2, ...。如果您相信您的数字生成器在 [0,1] 中是随机且不相关的，那么子序列 x1, x2, ... 应该在 [0,1] 中也是随机且不相关的。因此，对问题第一部分的相当简单的答案是投影到 x 或 y 轴上。也就是说，只需选择 x。
接下来考虑 x 和 y 之间的相关性。由于您没有指定相关性的性质，我们假设对轴进行简单缩放，例如 x' => [0, 0.5], y' => [0, 3.0]，然后进行旋转。缩放不会引入任何相关性，因为 x' 和 y' 仍然是独立的。您可以使用矩阵乘法轻松生成它：
M1*p = [x_scale, 0; 0, y_scale] * [x; y]
对于矩阵 M1 和点 p。您可以通过采用这种拉伸形式并将其旋转 theta 来引入相关性：
M2*M1*p = [cos(theta), sin(theta); -sin(theta), cos(theta)]*M1*p
将它们与 theta = pi/4 或 45 度放在一起，您可以看到较大的 y 值与较大的 x 值相关：
cos_t = sin_t = cos(pi/4); % at 45 degrees, sin(t) = cos(t) = 1/sqrt(2)
M2 = [cos_t, sin_t; -sin_t, cos_t];
M1 = [0.5, 0.0; 0.0, 3.0];
p = random(2,1000);
p_prime = M2*M1*p;
plot(p_prime(1)', p_prime(2)', '.');
axis('equal');
结果图*以 45 度角显示了一个均匀分布的数字带：
使用剪切可以进行进一步的变换，如果你很聪明，可以使用平移（OpenGL 使用 4x4 变换矩阵，因此平移可以表示为线性变换矩阵，在变换步骤之前添加一个额外的维度，并在完成之前删除） .
Mk*...*M1 p = p_prime给定一个已知的仿射相关结构，您可以通过求解p从随机点 (x',y') 转换回点 (x,y)，其中 x 和 y 在 [0,1] 中是独立的p = inv(Mk*...*M1) * p_prime，在哪里p=[x;y]。同样，只需选择 x，它将在 [0,1] 中是一致的。如果变换矩阵是奇异的，则这不起作用，例如，如果您将投影矩阵 Mj 引入混合（尽管如果投影是第一步，您仍然可以恢复）。
* 您可能会注意到绘图来自 python 而不是 matlab。我现在没有 matlab 或 octave 坐在我面前，所以我希望我的语法细节是正确的。
您可以从 f(x,y)=z 计算希尔伯特曲线。基本上它是一个汉密尔顿路径遍历。您可以在 Nick 的空间索引希尔伯特曲线四叉树博客中找到很好的描述。或者看看单调 n 元格雷码。我在 php 中基于 Nick 的博客编写了一个实现：http: //monstercurves.codeplex.com。
我将只关注你的最后一点
(3)使用 Hilbert 空间填充曲线对 draw in
[0,1]x[0,1]中的每个 draw in 进行变换[0,1]：在 Hilbert 曲线映射下，draw in[0,1]x[0,1]应该是 中的一个（或多个由于超射性）点的图像[0,1]。我想选择这些点之一。Matlab 中是否有任何预先构建的软件包可以做到这一点？
据我所知，Matlab 中没有预先构建的包这样做，但好消息是可以从 MATLAB 调用维基百科上的代码，并且就像将转换例程与网关函数放在一起一样简单在一个xy2d.c文件中：
#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]);
}
并用mex('xy2d.c').
上述实现
[...] 假设一个正方形被分成n个单元格，n为 2 的幂，具有整数坐标，左下角为 (0,0)，上角为 ( n -1, n -1)右上角。
实际上，在应用映射之前需要一个离散化步骤。与每个离散化问题一样，明智地选择精度至关重要。下面的片段将所有内容放在一起。
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')