Подтвердить что ты не робот

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

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

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

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

200px-2D-simplex.svg.png

(Изображение из http://en.wikipedia.org/wiki/Simplex.)

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

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

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

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

Связанные вопросы

4b9b3361

Ответ 1

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

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

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

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

Ответ 2

После небольшого копания я нашел эту страницу, которая дает приятную реализацию Дирихле-дистрибутива. Оттуда кажется, что было бы довольно просто следовать методу Википедии 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

Ответ 3

Здесь хорошая краткая реализация второго алгоритма из Wikipedia:

SimplexSample[n_] := [email protected]# - [email protected]# &[[email protected][{0,1}, RandomReal[{0,1}, n-1]]]

Это адаптировано здесь: http://www.mofeel.net/1164-comp-soft-sys-math-mathematica/14968.aspx (Первоначально у него был Union вместо Sort @Join - последний немного быстрее.)

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

Ответ 4

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

Времени, это немного накладные расходы, чтобы сначала выполнить полный дистрибутив Дирихле, так как все, что вам действительно нужно, это n random Gamma[1,1] samples. Сравнить ниже
Простая реализация

SimplexSample[n_, opts:OptionsPattern[RandomReal]] :=
  (#/Total[#])& @ RandomReal[GammaDistribution[1,1],n,opts]

Полная реализация Dirichlet

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

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}].

Ответ 5

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

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

<я > х <суб > 1суб >= 1- <я > R <суб > 1суб > 1/п-1

<я > х <суб > <я > ксуб >= (1- & Sigma; <суб >= 1суб > <я > к <я > х <суб > суб > ) (1- <я > R <суб > <я > к 1/nk), k = 2,..., n-1

<я > х <суб > <я > псуб >= 1- & Sigma; <суб >= 1суб > п-1 <я > х <суб > суб >

Где R_i - случайное число от 0 до 1.

Теперь я пытаюсь создать алгоритм для генерации случайных однородных выборок из ограниченного симплекса. Это пересечение между симплексом и выпуклым телом.