在Matlab中将双变量绘图转换为单变量绘图
Convert a bivariate draw in a univariate draw in Matlab
我想到了以下在 Matlab 中 运行 的实验，我请求帮助来实现步骤 (3)。任何建议将不胜感激。
(1) 考虑随机变量X和Y均均匀分布在[0,1]
上
(2) 从X和Y的联合分布得出N实现，假设X和Y是独立的（意思是X 和 Y 在 [0,1]x[0,1] 上均匀联合分布）。每次抽奖将在 [0,1]x[0,1] 内进行。
(3)使用希尔伯特space填充曲线将[0,1]x[0,1]中的每个绘图变换成[0,1]中的绘图：在希尔伯特曲线映射下，[=20=中的绘图] 应该是 [0,1] 中一个（或多个，因为满射）点的图像。我想选择其中一个点。 Matlab 中是否有任何预构建的包可以做到这一点？
我找到了 this 的答案，我认为这不是我想要的，因为它解释了如何获得绘图的希尔伯特值（从曲线起点到选取点的曲线长度）
在维基百科上，我找到了 this C 语言代码（从 (x,y) 到 d），这又不能满足我的问题。
您可以根据 f(x,y)=z 计算希尔伯特曲线。基本上它是哈密顿路径遍历。您可以在 Nick 的空间索引希尔伯特曲线四叉树博客中找到很好的描述。或者看看单调 n 元格雷码。我在 php:http://monstercurves.codeplex.com.
中基于 Nick 的博客编写了一个实现
编辑 这个答案没有解决问题的更新版本，它明确询问有关构建希尔伯特曲线的问题。相反，这个答案解决了一个关于构建双射映射的相关问题，以及与均匀分布的关系。
您的问题定义不明确。如果您只需要结果分布是均匀的，没有什么可以阻止您简单地选择 f:(X,Y)->X。无论 X 和 Y 是否相关，结果都是统一的。从你的post我只能假设你想要的，事实上，是为了让结果转换成为bijective，或者尽可能接近它机器精度限制。
值得注意的是，除非您需要最能保留 locality 的算法（这显然 not 需要结果分布双射的，更不用说制服了)，没有必要费心构建你在问题中提到的 Hilbert curves 。它们与解决方案的关系与任何其他 space 填充曲线一样多，并且计算量非常大。
因此，假设您正在寻找双射映射，您的问题等同于询问 [unit] 正方形中的点集是否具有相同的 cardinality as the set of points in a [unit] line segment, and if it is, how to construct that bijection, i.e. 1-to-1 correspondence. The intuition says the square should have a higher cardinality, and Cantor spent 3 years trying to prove that, eventually proving quite the opposite - that these sets are in fact equinumerous。他对自己的发现感到非常惊讶，以至于他写道：
I see it, but I don't believe it!
满足**此标准的最常提到的双射如下。用十进制形式表示 x 和 y，即 x=0. x1 x2 x3 x4 x5... 和 y=0 . y1 y2 y3 y4 y5...，令 f:(X,Y)->Z 为 z=0. x1 y1 x2 y2 x3 y3 x4 y4 x5 y5... ，即交替两个数字的小数点。双射背后的想法是微不足道的，尽管严格的证明需要相当多的先验知识。
** 需要注意的是，如果我们采用例如x = 1/3 = 0.33333...和y = 1/5 = 0.199999... = 0.200000...，我们可以看到对应的有两个序列：z = 0.313939393939...和z = 0.323030303030...。为了克服这个障碍，我们必须证明 adding a countable set to an uncountable one does not change the cardinality of the latter.
实际上我们必须处理机器精度而不是纯数学，严格来说这意味着这两个集合实际上是有限的因此不是等数的（假设您以与原始数字相同的精度存储结果）。这意味着我们只是被迫做一些假设并丢失一些信息，例如，在这种情况下，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'
Edit 这回答了最初的转换请求 f(x,y) -> t ~ U[0,1] given x,y ~ 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), ... . 你要的是对 (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 变换矩阵，因此平移可以表示为线性变换矩阵，在变换步骤之前添加一个额外的维度并删除在他们完成之前）。
给定一个已知的仿射相关结构，您可以从随机点 (x',y') 转换回点 (x,y)，其中 x 和 y 在 [0,1] 中是独立的，方法是求解 Mk*...*M1 p = p_prime 对于 p，或者等效地，通过设置 p = inv(Mk*...*M1) * p_prime，其中 p=[x;y]。同样，只需选择 x，它将在 [0,1] 中统一。如果变换矩阵是奇异的，则这不起作用，例如，如果您将投影矩阵 Mj 引入混合中（尽管如果投影是第一步，您仍然可以恢复）。
* 您可能会注意到该图来自 python 而不是 matlab。我现在面前没有 matlab 或 octave，所以我希望我的语法细节是正确的。
我只会关注你的最后一点
(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?
据我所知，在 Matlab 中没有预构建的包可以做到这一点，但好消息是维基百科上的 code 可以从 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')编译它。
以上实现
[...] 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.
实际上，在应用映射之前需要 离散化 步骤。与每个离散化问题一样，明智地选择精度至关重要。下面的代码片段将所有内容放在一起。
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')
我想到了以下在 Matlab 中 运行 的实验，我请求帮助来实现步骤 (3)。任何建议将不胜感激。
(1) 考虑随机变量X和Y均均匀分布在[0,1]
(2) 从X和Y的联合分布得出N实现，假设X和Y是独立的（意思是X 和 Y 在 [0,1]x[0,1] 上均匀联合分布）。每次抽奖将在 [0,1]x[0,1] 内进行。
(3)使用希尔伯特space填充曲线将[0,1]x[0,1]中的每个绘图变换成[0,1]中的绘图：在希尔伯特曲线映射下，[=20=中的绘图] 应该是 [0,1] 中一个（或多个，因为满射）点的图像。我想选择其中一个点。 Matlab 中是否有任何预构建的包可以做到这一点？
我找到了 this 的答案，我认为这不是我想要的，因为它解释了如何获得绘图的希尔伯特值（从曲线起点到选取点的曲线长度）
在维基百科上，我找到了 this C 语言代码（从 (x,y) 到 d），这又不能满足我的问题。
您可以根据 f(x,y)=z 计算希尔伯特曲线。基本上它是哈密顿路径遍历。您可以在 Nick 的空间索引希尔伯特曲线四叉树博客中找到很好的描述。或者看看单调 n 元格雷码。我在 php:http://monstercurves.codeplex.com.
中基于 Nick 的博客编写了一个实现编辑 这个答案没有解决问题的更新版本，它明确询问有关构建希尔伯特曲线的问题。相反，这个答案解决了一个关于构建双射映射的相关问题，以及与均匀分布的关系。
您的问题定义不明确。如果您只需要结果分布是均匀的，没有什么可以阻止您简单地选择 f:(X,Y)->X。无论 X 和 Y 是否相关，结果都是统一的。从你的post我只能假设你想要的，事实上，是为了让结果转换成为bijective，或者尽可能接近它机器精度限制。
值得注意的是，除非您需要最能保留 locality 的算法（这显然 not 需要结果分布双射的，更不用说制服了)，没有必要费心构建你在问题中提到的 Hilbert curves 。它们与解决方案的关系与任何其他 space 填充曲线一样多，并且计算量非常大。
因此，假设您正在寻找双射映射，您的问题等同于询问 [unit] 正方形中的点集是否具有相同的 cardinality as the set of points in a [unit] line segment, and if it is, how to construct that bijection, i.e. 1-to-1 correspondence. The intuition says the square should have a higher cardinality, and Cantor spent 3 years trying to prove that, eventually proving quite the opposite - that these sets are in fact equinumerous。他对自己的发现感到非常惊讶，以至于他写道：
I see it, but I don't believe it!
满足**此标准的最常提到的双射如下。用十进制形式表示 x 和 y，即 x=0. x1 x2 x3 x4 x5... 和 y=0 . y1 y2 y3 y4 y5...，令 f:(X,Y)->Z 为 z=0. x1 y1 x2 y2 x3 y3 x4 y4 x5 y5... ，即交替两个数字的小数点。双射背后的想法是微不足道的，尽管严格的证明需要相当多的先验知识。
** 需要注意的是，如果我们采用例如x = 1/3 = 0.33333...和y = 1/5 = 0.199999... = 0.200000...，我们可以看到对应的有两个序列：z = 0.313939393939...和z = 0.323030303030...。为了克服这个障碍，我们必须证明 adding a countable set to an uncountable one does not change the cardinality of the latter.
实际上我们必须处理机器精度而不是纯数学，严格来说这意味着这两个集合实际上是有限的因此不是等数的（假设您以与原始数字相同的精度存储结果）。这意味着我们只是被迫做一些假设并丢失一些信息，例如，在这种情况下，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'
Edit 这回答了最初的转换请求 f(x,y) -> t ~ U[0,1] given x,y ~ 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), ... . 你要的是对 (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 变换矩阵，因此平移可以表示为线性变换矩阵，在变换步骤之前添加一个额外的维度并删除在他们完成之前）。
给定一个已知的仿射相关结构，您可以从随机点 (x',y') 转换回点 (x,y)，其中 x 和 y 在 [0,1] 中是独立的，方法是求解 Mk*...*M1 p = p_prime 对于 p，或者等效地，通过设置 p = inv(Mk*...*M1) * p_prime，其中 p=[x;y]。同样，只需选择 x，它将在 [0,1] 中统一。如果变换矩阵是奇异的，则这不起作用，例如，如果您将投影矩阵 Mj 引入混合中（尽管如果投影是第一步，您仍然可以恢复）。
* 您可能会注意到该图来自 python 而不是 matlab。我现在面前没有 matlab 或 octave，所以我希望我的语法细节是正确的。
我只会关注你的最后一点
(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?
据我所知，在 Matlab 中没有预构建的包可以做到这一点，但好消息是维基百科上的 code 可以从 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')编译它。
以上实现
[...] 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.
实际上，在应用映射之前需要 离散化 步骤。与每个离散化问题一样，明智地选择精度至关重要。下面的代码片段将所有内容放在一起。
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')