Равномерно случайная выборка из n-мерного единичного симплекса - PullRequest
Новая
       18

Равномерно случайная выборка из n-мерного единичного симплекса

19 голосов
/ 10 июня 2010

Равномерная выборка случайным образом из n-мерного единичного симплекса - причудливый способ сказать, что вы хотите n случайных чисел, таких что

  • все они неотрицательны,
  • они составляют одну, а
  • одинаково вероятны все возможные векторы из n неотрицательных чисел, сумма которых равна единице.

В случае n = 2 вы хотите произвести равномерную выборку из сегмента линии x + y = 1 (т. Е. Y = 1-x), который находится в положительном квадранте. В случае n = 3 вы производите выборку из треугольной части плоскости x + y + z = 1, которая находится в положительном октанте R3:

image

(Image from http://en.wikipedia.org/wiki/Simplex.)

Обратите внимание, что сбор n равномерных случайных чисел и затем их нормализация, чтобы они суммировались в одно, не работает. Вы в конечном итоге смещены в сторону менее экстремальных чисел.

Аналогичным образом, выбор n-1 равномерных случайных чисел, а затем принятие n-го равным одному минус сумма их также приводит к смещению.

Википедия дает два алгоритма, чтобы сделать это правильно: http://en.wikipedia.org/wiki/Simplex#Random_sampling (Хотя второй в настоящее время утверждает, что он верен только на практике, а не в теории. Я надеюсь это исправить или уточнить, когда лучше пойму это. Сначала я застрял в «ПРЕДУПРЕЖДЕНИИ: такие-то и такие-то бумажные заявления на этой странице Википедии неверно следующее, и кто-то еще превратил это в предупреждение «работает только на практике».)

Наконец, вопрос: Что вы считаете лучшей реализацией симплексной выборки в Mathematica (желательно с эмпирическим подтверждением того, что это правильно)?

Похожие вопросы

Ответы [ 5 ]

10 голосов
/ 21 июня 2010

Этот код может работать: 

samples[n_] := Differences[Join[{0}, Sort[RandomReal[Range[0, 1], n - 1]], {1}]]

В основном вы просто выбираете n - 1 мест на интервале [0,1], чтобы разделить его, а затем берете размер каждой части, используя Differences.

Быстрый показ Timing показывает, что это немного быстрее, чем первый ответ Януса.

8 голосов
/ 10 июня 2010

Немного покопавшись, я нашел эту страницу , которая дает хорошую реализацию Dirichlet Distribution. Оттуда кажется, что было бы довольно просто следовать методу Википедии 1. Это кажется лучшим способом сделать это.

В качестве теста: 

In[14]:= RandomReal[DirichletDistribution[{1,1}],WorkingPrecision->25]
Out[14]= {0.8428995243540368880268079,0.1571004756459631119731921}
In[15]:= Total[%]
Out[15]= 1.000000000000000000000000

Сюжет из 100 образцов:

альтернативный текст http://www.public.iastate.edu/~zdavkeos/simplex-sample.png

6 голосов
/ 10 июня 2010

Вот хорошая краткая реализация второго алгоритма из Wikipedia : 

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 (Первоначально вместо Sort @ Join у него был Union - последний немного быстрее.)

(См. Комментарии для некоторых доказательств того, что это правильно!)

6 голосов
/ 10 июня 2010

Я с zdav: распределение Dirichlet кажется наиболее простым, а алгоритм выборки распределения Dirichlet, на который ссылается zdav, также представлен на странице Википедии в Dirichlet .

С точки зрения реализации, сначала стоит выполнить полный дистрибутив Дирихле, так как все, что вам действительно нужно, это 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}].

2 голосов
/ 12 июня 2014

Я создал алгоритм для равномерной генерации случайных чисел по симплексу. Вы можете найти подробности в документе по следующей ссылке: http://www.tandfonline.com/doi/abs/10.1080/03610918.2010.551012#.U5q7inJdVNY

Вкратце, вы можете использовать следующие формулы рекурсии, чтобы найти случайные точки над n-мерным симплексом:

х 1 = 1- R 1 1 / п-1

х * * K одна тысяча двадцать-один * ** +1023 тысяча двадцать две * = (1- & Sigma; = 1 * тысяча двадцать-пять * к * ** тысяча двадцать восемь одна тысяча двадцать-девять * х * ** 1032 тысяча тридцать-одна * г ) (* 1- +1036 * R * +1039 * к 1 / nk ), k = 2, ..., n-1

* +1045 ** * х одна тысяча сорок-шесть * ** +1047 +1048 * N * * = тысяче пятьдесят-один 1- & Sigma; * 1 052 * = 1 п-1 * ** одна тысяча пятьдесят пять тысяча пятьдесят-шесть * х я Где R_i - это случайное число от 0 до 1. Теперь я пытаюсь создать алгоритм для генерации случайных равномерных выборок из ограниченного симплекса. Это пересечение симплекса и выпуклого тела.
...