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

Question

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

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

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

Спасибо.

Answer 1

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

%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 является квадратным. Я предлагаю избегать симметричных величин на стадии разработки.