Выберите однородно наугад от n-мерного симплекса единицы
Выборка однородно наугад от n-мерного симплекса единицы является необычным способом сказать желание n случайных чисел, таким образом что
- они являются все неотрицательными,
- они суммируют одному, и
- каждый возможный вектор n неотрицательных чисел, которые суммируют, каждый одинаково вероятен.
В n=2 случае Вы хотите выбрать однородно от сегмента строки x+y=1 (т.е., y=1-x), который находится в положительном квадранте. В n=3 случае Вы выбираете от части треугольной формы плоскости x+y+z=1, который находится в положительном октанте R3:
(Изображение из http://en.wikipedia.org/wiki/Simplex.)
Обратите внимание, что, выбирая n универсальные случайные числа и затем нормализуя их так они суммируют, каждый не работает. Вы заканчиваете с предвзятостью к менее экстремальным числам.
Точно так же выбор n-1 универсальные случайные числа и затем взятие энного, чтобы быть один минус сумма их также представляют предвзятость.
Википедия дает два алгоритма, чтобы сделать это правильно: http://en.wikipedia.org/wiki/Simplex#Random_sampling (Хотя второй в настоящее время утверждает, что только был корректен на практике, не в теории. Я надеюсь очистить это или разъяснить его, когда я понимаю это лучше. Я первоначально всунул "ПРЕДУПРЕЖДЕНИЕ: такая-то и такая-то бумага утверждает, что следующее является неправильным" на той странице Wikipedia, и кто-то еще превратил его в "работы только в практике" протест.)
Наконец, вопрос: Что Вы рассматриваете лучшей реализацией выборки симплекса в Mathematica (предпочтительно с эмпирическим подтверждением, что это корректно)?
Связанные вопросы
4 ответа
Этот код может работать:
samples[n_] := Differences[Join[{0}, Sort[RandomReal[Range[0, 1], n - 1]], {1}]]
Обычно вы просто выбираете n - 1 мест в интервале [0,1] , чтобы разделить его, затем возьмите размер каждой части, используя Различия .
Быстрый просмотр Timing показывает, что это немного быстрее, чем первый ответ Януса.
Я согласен с zdav: распределение Дирихле кажется самым простым способом, а алгоритм выборки распределения Дирихле, на который ссылается zdav, также представлен на странице Распределение Дирихле в Википедии.
С точки зрения реализации, делать сначала полное распределение Дирихле немного накладно, поскольку все, что вам действительно нужно - это n случайных Gamma[1,1] выборок. Сравните ниже
Простая реализация
SimplexSample[n_, opts:OptionsPattern[RandomReal]] :=
(#/Total[#])& @ RandomReal[GammaDistribution[1,1],n,opts]
Полная реализация Дирихле
DirichletDistribution/:Random`DistributionVector[
DirichletDistribution[alpha_?(VectorQ[#,Positive]&)],n_Integer,prec_?Positive]:=
Block[{gammas}, gammas =
Map[RandomReal[GammaDistribution[#,1],n,WorkingPrecision->prec]&,alpha];
Transpose[gammas]/Total[gammas]]
SimplexSample2[n_, opts:OptionsPattern[RandomReal]] :=
(#/Total[#])& @ RandomReal[DirichletDistribution[ConstantArray[1,{n}]],opts]
Время
Timing[Table[SimplexSample[10,WorkingPrecision-> 20],{10000}];]
Timing[Table[SimplexSample2[10,WorkingPrecision-> 20],{10000}];]
Out[159]= {1.30249,Null}
Out[160]= {3.52216,Null}
Таким образом, полное Дирихле в 3 раза медленнее. Если вам нужно m>1 точек выборки за раз, вы, вероятно, можете выиграть еще больше, делая (#/Total[#]&)/@RandomReal[GammaDistribution[1,1],{m,n}].
Вот хорошая краткая реализация второго алгоритма из Википедии :
SimplexSample[n_] := Rest@# - Most@# &[Sort@Join[{0,1}, RandomReal[{0,1}, n-1]]]
Адаптировано отсюда: http: // www. mofeel.net/1164-comp-soft-sys-math-mathematica/14968.aspx (Первоначально в нем было Union вместо Sort @ Join - последнее работает немного быстрее.)
(Некоторые доказательства того, что это правильно, см. В комментариях!)
Немного покопавшись, я нашел эту страницу, где дана хорошая реализация распределения Дирихле. Отсюда кажется, что будет довольно просто следовать методу 1 Википедии. Кажется, что это лучший способ сделать это.
В качестве теста:
In[14]:= RandomReal[DirichletDistribution[{1,1}],WorkingPrecision->25]
Out[14]= {0.8428995243540368880268079,0.1571004756459631119731921}
In[15]:= Total[%]
Out[15]= 1.000000000000000000000000
График из 100 образцов:
alt text http://www.public.iastate.edu/~zdavkeos/simplex-sample.png
