Преобразование двумерной ничьей в одномерном розыгрыше в Matlab

Question

Я имею в виду следующий эксперимент, который нужно запустить в Matlab, и я прошу помочь выполнить шаг (3). Любое предложение было бы очень оценено.

(1) Рассмотрим случайные величины X и Y и равномерно распределены по [0,1]

(2) Draw N реализацию от совместного распределения X и Y в предположении, что X и Y независимы (это означает, что X и Y равномерно распределены по [0,1]x[0,1]). Каждая ничья будет в [0,1]x[0,1].

(3) Преобразование каждый дро в [0,1]x[0,1] в розыгрыше в [0,1] с использованием кривого заполнения пространства Гильберта: при отображении кривого Гильберта, втягивания [0,1]x[0,1] должна быть изображением одного (или более из-за сюрьективность) точек (ы) в [0,1]. Я хочу выбрать один из этих пунктов. Есть ли какой-либо готовый пакет в Matlab?

Я нашел this ответ, который я не думаю, что делает то, что я хочу, как он объясняет, как получить значение Гильберта жеребьевке (длина кривой от начала кривой до подобранной точки)

На википедии Я нашел this код на языке C (от (x,y) до d), который, опять же, не отвечает моему вопросу.

Answer 1

я остановлюсь только на последней точке

(3) преобразование каждого розыгрыша в [0,1]x[0,1] вничью в [0,1] с использованием кривой заполнения пространства Гильберта: при отображении кривой Гильберта, втягивания [0,1]x[0,1] должен быть изображение одной (или более из-за сюръективности) точек (точек) в [0,1]. Я хочу выбрать один из этих пунктов. Есть ли в Matlab готовый пакет?

Насколько я знаю, там не заранее построенные пакеты в 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').

выше реализация

[...] предполагает квадрат, разделенный на п по п клеток, для н степенью 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') 

Answer 2

Вы можете вычислить гильбертовую кривую из f (x, y) = z. В основном это гамильтоновский обход пути. Вы можете найти хорошее описание в блоге Quadtree кривой пространственного индекса Ника. Или взгляните на монотонный n-мерный серый код. Я написал реализацию, основанную на блоге Ника в php: http://monstercurves.codeplex.com.

Answer 3

EDITВ этом ответе не рассматривается обновленная версия вопроса, в которой явно задается вопрос о построении кривой Гильберта. Вместо этого в этом ответе рассматривается связанный вопрос о построении биективного отображения и отношение к равномерному распределению.

Ваша проблема не совсем определена. Если вам нужно только, чтобы получившееся распределение было равномерным, ничто не мешает вам просто выбрать f:(X,Y)->X. Результат будет равномерным независимо от того, коррелированы ли X и Y. Из вашего поста я могу только предположить, что на самом деле вы хотите, чтобы получившееся преобразование было bijective, или как можно ближе к ним, с учетом заданных ограничений точности машины.

Стоит отметить, что, если не нужен алгоритм, который лучше всего в сохранении местности (что явно не требуется для результирующего распределения биективных, не говоря уже о униформе), нет необходимости беспокоиться построениями Hilbert curves, что вы упоминание в вашем вопросе. Они так же связаны с решением, как и любая другая кривая заполнения пространства, и невероятно интенсивно вычислительны.

Так что, если вы ищете биективное сопоставление, ваш вопрос эквивалентен запросу , имеет ли набор точек в квадрате [unit] тот же номер cardinality как набор точек в сегменте линии [unit] и если да, то как построить эту биекцию, т. е. соответствие 1-к-1. Интуиция говорит, что квадрат должен иметь более высокую мощность, а Cantor потратил 3 years, пытаясь доказать, что, в конце концов, доказывая совершенно противоположное - эти наборы фактически являются equinumerous. Он был настолько удивлен его открытием, что писал:

Я вижу это, но я этому не верю!

Наиболее часто упоминаемый биекция, выполняющая ** эти критерии, заключается в следующем. Представляют x и y в их десятичной форме, то есть x = 0.x1x2x3x4x5 ... и у = 0.у1 у2у3у4 Y5 ..., и пусть f:(X,Y)->Z быть г = 0.x1у1 x2у2 х3у3x4 у4x5у5 ..., т.е. чередуя десятичные двух чисел. Идея биекции тривиальна, хотя строгое доказательство требует довольно много предварительного знания.

** Предостережение таково, что, если мы возьмем, например, 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' 

Answer 4

РедактироватьЭто отвечает на первоначальный запрос для преобразования Р (х, у) -> т ~ U [0,1] Данные х, у ~ U [0,1], а также для х и у коррелированы. В обновленном вопросе конкретно задана кривая Гильберта 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. Вы можете ввести корреляцию, приняв эту растянутую форму и вращая ее тета:

M2*M1*p = [cos(theta), sin(theta); -sin(theta), cos(theta)]*M1*p 

Собираем все вместе с тета = пи/4, или 45 градусов, вы можете увидеть, что большие значения у коррелируют с более крупными значения х:

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 или октавы, поэтому я надеюсь, что я правильно понял синтаксические данные.