Генерация случайных выборок из произвольной дискретной функции плотности вероятности в Matlab

Question

У меня есть произвольная функция плотности вероятности, дискретизированная как матрица в Matlab, это означает, что для каждой пары x, y вероятность сохраняется в матрице: A(x,y) = вероятность

Это матрица 100x100, и я хотел бы иметь возможность генерировать случайные выборки двух измерений (x, y) из этой матрицы, а также, если возможно, иметь возможность вычислять среднее и другие моменты PDF. Я хочу сделать это, потому что после повторной выборки я хочу подогнать образцы к приближенной модели смеси Гаусса.

Я искал повсюду, но я не нашел ничего такого особенного, как этот. Я надеюсь, что вы сможете мне помочь.

Спасибо.

Answer 1

Если у вас действительно есть дискретная, вероятно, функция плотности, определяемая как A (в отличие от непрерывной функции плотности вероятности, которая просто описывается A), вы можете "обмануть", превратив свою двумерную задачу в одномерную.

%define the possible values for the (x,y) pair
row_vals = [1:size(A,1)]'*ones(1,size(A,2));  %all x values
col_vals = ones(size(A,1),1)*[1:size(A,2)];  %all y values

%convert your 2D problem into a 1D problem
A = A(:);
row_vals = row_vals(:);
col_vals = col_vals(:);

%calculate your fake 1D CDF, assumes sum(A(:))==1
CDF = cumsum(A); %remember, first term out of of cumsum is not zero

%because of the operation we're doing below (interp1 followed by ceil)
%we need the CDF to start at zero
CDF = [0; CDF(:)];

%generate random values
N_vals = 1000;  %give me 1000 values
rand_vals = rand(N_vals,1);  %spans zero to one

%look into CDF to see which index the rand val corresponds to
out_val = interp1(CDF,[0:1/(length(CDF)-1):1],rand_vals); %spans zero to one
ind = ceil(out_val*length(A));

%using the inds, you can lookup each pair of values
xy_values = [row_vals(ind) col_vals(ind)];

Надеюсь, это поможет!

чип

Answer 2

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

Поскольку половина ответа лучше, чем отсутствие ответа, вот как вы можете вычислить среднее и старшее моменты численно с помощью matlab:

%generate some dummy input
xv=linspace(-50,50,101);
yv=linspace(-30,30,100);
[x y]=meshgrid(xv,yv);

%define a discretized two-hump Gaussian distribution
A=floor(15*exp(-((x-10).^2+y.^2)/100)+15*exp(-((x+25).^2+y.^2)/100));
A=A/sum(A(:)); %normalized to sum to 1

%plot it if you like
%figure;
%surf(x,y,A)

%actual half-answer starts here    

%get normalized pdf
weight=trapz(xv,trapz(yv,A));
A=A/weight; %A normalized to 1 according to trapz^2

%mean
mean_x=trapz(xv,trapz(yv,A.*x));
mean_y=trapz(xv,trapz(yv,A.*y));

Итак, дело в том, что вы можете выполнить двойной интеграл на прямоугольной сетке, используя два последовательных вызова trapz, Это позволяет вам вычислять интеграл от любой величины, которая имеет ту же форму, что и ваша сетка, но недостатком является то, что векторные компоненты должны вычисляться независимо. Если вы хотите вычислить только то, что может быть параметризовано x а также y (которые, естественно, имеют такой же размер, как и у вас), тогда вы можете обойтись без дополнительных размышлений.

Вы также можете определить функцию для интеграции:

function res=trapz2(xv,yv,A,arg)

if ~isscalar(arg) && any(size(arg)~=size(A))
    error('Size of A and var must be the same!')
end

res=trapz(xv,trapz(yv,A.*arg));

end

Таким образом, вы можете вычислить такие вещи, как

weight=trapz2(xv,yv,A,1);
mean_x=trapz2(xv,yv,A,x);

ПРИМЕЧАНИЕ: причина, по которой я использовал сетку 101x100 в этом примере, заключается в том, что двойной вызов trapz должны быть выполнены в правильном порядке. Если вы меняете xv а также yv в вызовах вы получаете неправильный ответ из-за несоответствия определению A, но это не будет очевидно, если A квадратный Я предлагаю избегать симметричных величин на стадии разработки.