Выборка положительных весов с единицей суммы и ограничением равенства

Для равномерной выборки по набору {(a1, a2, a3, a4 | a2=a3, a1+a2+a3+a4=1, a1>0, a2>0, a3>0, a4>0} я бы сначала выбрал значение для a2 (которое равно a3). Для этого нам нужно знать распределение этого значения. Если a2 = a3 = r, то имеем a1+a4 = 1-2r; для положительных a1 и a4 существует отрезок длины (1-2k)*sqrt(2), содержащий все допустимые значения a1 и a4. При интегрировании вероятность того, что a2 равна k или меньше, равна 4(k - k^2). Более подробно:

Prob (a2 <= k) = Integral(0 to k) (1-2r)*sqrt(2) dr / Integral(0 to 0.5) (1-2r)*sqrt(2) dr
               = ((k-k^2)*sqrt(2)) / (sqrt(2)/4)
               = 4k - 4k^2

Таким образом, мы можем выбрать значения для a2, выбрав равномерно распределенное значение u~U(0, 1) и установив a2 равным значению k, для которого 4k - 4k^2 = u. Решая по квадратичной формуле, получаем:

a2 = 0.5 * (1 - sqrt(1-u))

В R мы можем выбрать 1000 значений для a2 с помощью:

set.seed(144)
a2 <- 0.5 * (1 - sqrt(1 - runif(1000)))
a3 <- a2

Учитывая фиксированное значение a2 = a3 = k, значение a1 равномерно распределяется в [0, 1-2k]:

a1 <- runif(1000) * (1 - 2*a2)

После указания a1, a2 и a3 существует только одно возможное значение для a4:

a4 <- 1 - a1 - a2 - a3

Мы можем взглянуть на некоторые из наших выборочных значений:

head(cbind(a1, a2, a2, a4))
#              a1         a2         a2         a4
# [1,] 0.83455239 0.01251016 0.01251016 0.14042729
# [2,] 0.02744599 0.22932773 0.22932773 0.51389856
# [3,] 0.45835472 0.23860119 0.23860119 0.06444291
# [4,] 0.36843649 0.14679703 0.14679703 0.33796946
# [5,] 0.35109881 0.08702039 0.08702039 0.47486041
# [6,] 0.02916818 0.19942616 0.19942616 0.57197949

Вот распределение значений a1 (обратите внимание, что по симметрии оно идентично распределению значений a4). Поскольку мы выбираем a1 равномерно в диапазоне [0, 1-2*a2], более низкие значения встречаются чаще, чем более высокие:

Вот распределение значений a2 (по определению это то же самое, что и распределение значений a3). Форма распределения аналогична a1, но максимальное значение равно 0,5:

Я пытался подумать об использовании распределения Дирихле,

Ну, для меня это похоже на распределение Дирихле.

но он не дает механизма, чтобы заставить две переменные быть равными.

но вам не нужно. На самом деле у вас есть три варианта распределения Дирихле — A, B, C, все >= 0, равномерно распределенные U (0,1), так что A + B + C = 1

После выборки (A, B, C) вы просто назначаете

a1 = A;
a2 = B/2.0;
a3 = B/2.0;
a4 = C;

Пожалуйста, посмотрите, как сэмплировать (ну, на Python)

Генерация N однородных случайных чисел, сумма которых равна M