Генерация случайного мултьтиномиального вектора

* Опубликован: 29.12.2005.

Solution ID:	 	100000320	
Product:	 	SPSS Base 	
Version:	 		
O/S:	 		
Тип вопроса:		 	Статистика	
Подтип вопроса:		 	Статистические распределения	
Заголовок:
Генерация мультиномиальных случайных переменных в SPSS 

Описание:

Вопрос. 
Как можно в SPSS сгенерировать вектор с мультиномиальным распределением с заданным числом наблюдений? 

Ответ. 
Если осуществлять выборки объёма n с возвращениями из генеральной совокупности,
имеющей k классов объектов, где k>2, то вектор из k чисел, означающих количество объектов,
выбранных из соответствующих классов, будет иметь мультиномиальное распределение.

Следующий макрос генерирует наблюдения и переменные с таким распределением.
Вы указываете число наблюдений случайной величины с мультиномиальным распределением,
которое необходимо сгенерировать (ncases), число классов объектов (classes),
объём выборки n (samsize) и доли или соотношение в абсолютных частотах классов
в генеральной совокупности (prop).

Если вы вводите соотношение в абсолютных частотах вместо долей,
то следует указать это, поставив значение параметра "norp = 1" при вызове макроса.
(Если вы введете соотношение, но не укажете, "norp = 1", макрос "определит" это
самостоятельно, если сумма долей всех классов превысит 1.001 и переведёт соотношение в доли.
Значение 0.001 как составляющая порога, приведённого выше, используется для того, чтобы 
предусмотреть возможные погрешности при складывании долей.

Используемый алгоритм генерации работает по аналогии извлечения шаров из урны. Он
описан в книге Johnson, Kotz, & Balakrishnan 
[(1997). "Discrete Multivariate Distributions", Wiley. pp. 68-69]. 

Как отмечают авторы, существуют и другие, более быстрые алгоритмы для генерации
мультиномиальных переменных. Однако, эти методоы не вполне оправдывают свою эффективность
из-за сложности их реализации в рамках структуры данных внутри SPSS. Алгорим следующий:

1. Доли классов в генеральной совокупности (pop1- popk) инициализируются из соответствующих
значений параметра prop. Количество объектов каждого класса в конкретной выборке (sam1-samk)
полагается равным 0. Сумма долей определяется как сумма pop1-popk и записывается в poptot. 
Если пропорции заданы корректно через параметр prop, poptot будет равен 1. Если было
задано соотношение классов в генеральной совокупности в абсолютных частотах, то poptot будет
больше 1.001 и значения pop1-popk пересчитаются в доли путём деления на poptot.


2. Для каждого из samsize наблюдений, подлежащих выборке в каждом случае:
а) вычисляется случайна величина, равномерно распределённая на интервале (0, 1) и её
значение записывается в Y.
б) Для каждого из k классов вычисляется значение psum как сумма долей классов в генеральной совокупности,
выбранных на данный момент. Если Y меньше или равна psum, но больше, чем psum для предыдущего класса, 
наблюдение считается выбранным из текущего класса. Количество выбранных объектов из этого класса
увеличивается на 1. 

NB! Алгоритм имеет следующее свойство. Если сумма заданных долей классов в генеральной совокупности будет
меньше 1, poptot будет также меньше 1 на ту же величину. Это будет означать, что на самом деле имеется
не k, а k+1 класс, где доля (k+1)-го класса определяется как 1 минус сумма заданных долей. Таким образом, 
количетсво выбранных объектов из (k+1)-го класса впоследствии можно определить как
samsize-(sum1+sum2+...+sumk) - А.Б.


* Макрос для генерации мультиномиального распределения. 
* Первый пример вызова задаёт 3 класса с долями .5, .3 и .2. 
* Делается выборка по 125 объектов с возвращением. 
* sam1 - sam3 содержат соответствующие частоты. 
* Генерируется 200 наблюдений случайного вектора, имеющего мультиномиальное распределение. 
* Второй пример задаёт 4 класса с соотношением частот как 20 к 10 к 30 к 20. 
* Делаются выборки по 30 объектов с возвращением. 
* sam1 - sam4 содержат частоты. 
* Параметр norp = 1 при вызове, что указывает на то, что заданы соотношения в абсолютных частотах.
* Генерируется 300 наблюдений вектора. 
* Третий пример вызова похож на второй, но параметр norp в явном виде не установлен. 
* Макрос сам определяет, что заданы соотношения в частотах на основании того, что 
* сумма элементов prop превосходит 1.001 . 
* . 
*************************************************************. 
* ОПРЕДЕЛЕНИЕ МАКРОСА.
define multnomg 
(ncases = !tokens(1) 
/classes = !tokens(1) 
/samsize = !tokens(1) 
/norp= !default (0) !tokens(1) 
/prop = !enclose('[',']') ). 
new file. 
input program . 
loop id = 1 to !ncases . 
+ compute np = !norp . 
vector pop (!classes , F8.5) sam (!classes , F8) . 
+ do repeat popn = pop1 to !concat('pop',!classes) 
/samn = sam1 to !concat('sam',!classes) 
/pc = !prop . 
+ compute popn = pc. 
+ compute samn = 0. 
+ end repeat. 
+ compute poptot = sum(pop1 to !concat('pop',!classes)). 
+ do if (np = 1 or poptot > 1.001). 
+ do repeat propn = pop1 to !concat('pop',!classes). 
compute propn = propn/poptot. 
+ end repeat. 
+ end if. 
+ loop #j = 1 to !samsize . 
+ compute y = uniform(1) . 
+ compute psum = 0. 
+ loop #k = 1 to !classes . 
+ compute psum = psum + pop(#k). 
+ if (y le psum and y gt (psum - pop(#k))) 
sam(#k) = sam(#k) + 1. 
+ end loop. 
+ end loop. 
+ end case. 
end loop. 
end file. 
end input program. 
execute. 
!enddefine . 
*.
*ПРИМЕРЫ ВЫЗОВОВ МАКРОСА.
multnomg ncases = 200 classes = 3 samsize = 125 
prop = [ 0.5 0.3 0.1 ] . 
multnomg ncases = 300 classes = 4 samsize = 30 norp = 1 
prop = [ 20 10 30 20 ] . 
multnomg ncases = 300 classes = 4 samsize = 30 
prop = [ 20 10 30 20 ] .