2016-10-25 4 views
0

Я нашел Matlab файл (https://in.mathworks.com/matlabcentral/fileexchange/9700-random-vectors-with-fixed-sum), который генерирует матрицу вероятностей с равномерной вероятностью, и каждый столбец имеет сумму 1. файл состоит в следующемКак сделать вероятностную матрицу с каждой строкой элементами суммы 1 и каждая строка имеет равномерную вероятность в MATLAB

function [x,v] = randfixedsum(n,m,s,a,b) 


%[x,v] = randfixedsum(n,m,s,a,b) 

% 
% This generates an n by m array x, each of whose m columns 
% contains n random values lying in the interval [a,b], but 
% subject to the condition that their sum be equal to s. The 
% scalar value s must accordingly satisfy n*a <= s <= n*b. The 
% distribution of values is uniform in the sense that it has the 
% conditional probability distribution of a uniform distribution 
% over the whole n-cube, given that the sum of the x's is s. 
% 
% The scalar v, if requested, returns with the total 
% n-1 dimensional volume (content) of the subset satisfying 
% this condition. Consequently if v, considered as a function 
% of s and divided by sqrt(n), is integrated with respect to s 
% from s = a to s = b, the result would necessarily be the 
% n-dimensional volume of the whole cube, namely (b-a)^n. 
% 
% This algorithm does no "rejecting" on the sets of x's it 
% obtains. It is designed to generate only those that satisfy all 
% the above conditions and to do so with a uniform distribution. 
% It accomplishes this by decomposing the space of all possible x 
% sets (columns) into n-1 dimensional simplexes. (Line segments, 
% triangles, and tetrahedra, are one-, two-, and three-dimensional 
% examples of simplexes, respectively.) It makes use of three 
% different sets of 'rand' variables, one to locate values 
% uniformly within each type of simplex, another to randomly 
% select representatives of each different type of simplex in 
% proportion to their volume, and a third to perform random 
% permutations to provide an even distribution of simplex choices 
% among like types. For example, with n equal to 3 and s set at, 
% say, 40% of the way from a towards b, there will be 2 different 
% types of simplex, in this case triangles, each with its own 
% area, and 6 different versions of each from permutations, for 
% a total of 12 triangles, and these all fit together to form a 
% particular planar non-regular hexagon in 3 dimensions, with v 
% returned set equal to the hexagon's area. 
% 
% Roger Stafford - Jan. 19, 2006 

% Check the arguments. 
if (m~=round(m))|(n~=round(n))|(m<0)|(n<1) 
error('n must be a whole number and m a non-negative integer.') 
elseif (s<n*a)|(s>n*b)|(a>=b) 
error('Inequalities n*a <= s <= n*b and a < b must hold.') 
end 

% Rescale to a unit cube: 0 <= x(i) <= 1 
s = (s-n*a)/(b-a); 

% Construct the transition probability table, t. 
% t(i,j) will be utilized only in the region where j <= i + 1. 
k = max(min(floor(s),n-1),0); % Must have 0 <= k <= n-1 
s = max(min(s,k+1),k); % Must have k <= s <= k+1 
s1 = s - [k:-1:k-n+1]; % s1 & s2 will never be negative 
s2 = [k+n:-1:k+1] - s; 
w = zeros(n,n+1); w(1,2) = realmax; % Scale for full 'double' range 
t = zeros(n-1,n); 
tiny = 2^(-1074); % The smallest positive matlab 'double' no. 
for i = 2:n 
tmp1 = w(i-1,2:i+1).*s1(1:i)/i; 
tmp2 = w(i-1,1:i).*s2(n-i+1:n)/i; 
w(i,2:i+1) = tmp1 + tmp2; 
tmp3 = w(i,2:i+1) + tiny; % In case tmp1 & tmp2 are both 0, 
tmp4 = (s2(n-i+1:n) > s1(1:i)); % then t is 0 on left & 1 on right 
t(i-1,1:i) = (tmp2./tmp3).*tmp4 + (1-tmp1./tmp3).*(~tmp4); 
end 

% Derive the polytope volume v from the appropriate 
% element in the bottom row of w. 
v = n^(3/2)*(w(n,k+2)/realmax)*(b-a)^(n-1); 

% Now compute the matrix x. 
x = zeros(n,m); 
if m == 0, return, end % If m is zero, quit with x = [] 
rt = rand(n-1,m); % For random selection of simplex type 
rs = rand(n-1,m); % For random location within a simplex 
s = repmat(s,1,m); 
j = repmat(k+1,1,m); % For indexing in the t table 
sm = zeros(1,m); pr = ones(1,m); % Start with sum zero & product 1 
for i = n-1:-1:1 % Work backwards in the t table 
e = (rt(n-i,:)<=t(i,j)); % Use rt to choose a transition 
sx = rs(n-i,:).^(1/i); % Use rs to compute next simplex coord. 
sm = sm + (1-sx).*pr.*s/(i+1); % Update sum 
pr = sx.*pr; % Update product 
x(n-i,:) = sm + pr.*e; % Calculate x using simplex coords. 
s = s - e; j = j - e; % Transition adjustment 
end 
x(n,:) = sm + pr.*s; % Compute the last x 

% Randomly permute the order in the columns of x and rescale. 
rp = rand(n,m); % Use rp to carry out a matrix 'randperm' 
[ig,p] = sort(rp); % The values placed in ig are ignored 
x = (b-a)*x(p+repmat([0:n:n*(m-1)],n,1))+a; % Permute & rescale x 


return 

но я хочу сгенерировать матрицу, которая имеет каждый элемент строки sum 1, и каждая строка имеет равномерную вероятность в матлабе. как это сделать с помощью вышеуказанной программы. для запуска выше программы я вызываю его в другой файл и настройки параметров т.е.

m=4;n=4; a=0; b=1.5;s=1; 
[x,v] = randfixedsum(n,m,s,a,b) 
+0

Создать столбец один и транспонировать его –

+0

нет, столбец не имеет равномерной вероятности, поэтому я не могу использовать, принимая транспонирование –

+0

Применить [this] (http://stackoverflow.com/q/8064629/2586922) (с '1' вместо '100') для каждой строки? –

ответ

0

Я бы сказал, что проще всего было, чтобы сгенерировать массив с данной функцией.

[x,v] = randfixedsum(n,m,s,a,b); 

Затем просто передайте результаты.

x = x'; 
+0

нет, столбец не имеет равномерной вероятности, поэтому я не могу использовать, принимая транспонирование –

0

создать случайную матрицу и разделить каждую строку на сумму элементов этой строки:

function result = randrowsum(m ,n) 
    rnd = rand(m,n); 
    rowsums = sum(rnd,2); 
    result = bsxfun(@rdivide, rnd, rowsums); 
end 

создать m * n случайную матрицу:

a=randrowsum(3,4) 

чек, если сумма каждой строки 1:

sum(a,2) 
+0

Я запускаю программу ниже, но это ошибка посева –

+0

m = 3; n = 3; a = randrowsum (m, n) Результат функции = a rnd = rand (m, n); rowsums = sum (rnd, 2); result = bsxfun (@rdivide, rnd, rowsums) end –

+0

@JayRam сохранить функцию в файле с именем 'randrowsum.m', а затем вызвать функцию – rahnema1

 Смежные вопросы

  • Нет связанных вопросов^_^