Генерация вектора «случайных» пропорций заданной длины в заданных пределах c

Question

Я хочу сгенерировать вектор заданной длины, например, n = 5. Каждое значение в векторе должно быть пропорцией (т. Е. Значением от 0 до 1), чтобы по n элементам они суммировались до 1.

К сожалению, у меня есть два вектора: один (mymins) определяет допустимые нижние границы каждой пропорции, а другой (mymaxs) определяет разрешенные верхние границы каждой пропорции.

В моем примере ниже желаемая пропорция для первого элемента может падать где-то между 0,3 и 0,9. А для последнего элемента желаемая пропорция может упасть между 0,05 и 0,7.

mymins <- c(0.3, 0.1, 0, 0.2, 0.05)
mymaxs <- c(0.9, 1, 1, 1, 0.7)

Давайте предположим, что mymins всегда «законны» (т. Е. Их сумма никогда не превышает 1).

Как я могу найти набор из 5 пропорций, чтобы все они составляли 1, но l ie в границах?

Вот что я попробовал:

n = 5
mydif <- mymaxs - mymins    # possible range for each proportion
myorder <- rank(mydif)      # order those differences from smallest to largest
mytarget <- sum(mydif)      # sum up the 5 ranges
x <- sort(runif(n))[myorder] # generate 5 random values an sort them in the order of mydif
x2 <- mymins + x / sum(x) * mytarget  # rescale random values to sum up to mytarget and add them to mymins
x3 <- x2/sum(x2)             # rescale x2 to sum up to 1

Как вы можете видеть, я не очень далеко - потому что после перемасштабирования некоторые значения выходят за допустимые границы.

Я, вероятно, должен также упомянуть, что мне нужна эта операция, чтобы быть быстрой - потому что я использую ее в Оптимизация l oop.

Я также пытался найти решение с помощью Optim, однако проблема в том, что он всегда находит одно и то же решение - и мне нужно генерировать РАЗНЫЕ решения каждый раз, когда я нахожу пропорциональность:

    myfun <- function(x) {
      x <- round(x, 4)
      abovemins <- x - mymins
      n_belowmins <- sum(abovemins < 0)
      if (n_belowmins > 0) return(100000)
      belowmax <- x - mymaxs
      n_abovemax <- sum(belowmax > 0)
      if (n_abovemax > 0) return(100000)
      mydist <- abs(sum(x) - 1)
      return(mydist)
    }

    myopt <- optim(par = mymins + 0.01, fn = myfun)
    myopt$par
    sum(round(myopt$par, 4))

Большое спасибо за ваши предложения!

Answer 1

Возможно, лучше подумать об этом по-другому. Ваши сэмплы на самом деле должны суммироваться до 0,35 (то есть 1 - сумма (mymins)), а затем добавляться к минимальным значениям

constrained_sample <- function(mymins, mymaxs)
{
 sizes <- mymaxs - mymins
 samp <- (runif(5) * sizes)
 samp/sum(samp) * (1 - sum(mymins)) + mymins
}

Это работает так:

constrained_sample(mymins, mymaxs)
#> [1] 0.31728333 0.17839397 0.07196067 0.29146744 0.14089459

Мы можем проверить это, запустив следующую команду l oop, которая выведет сообщение на консоль, если какой-либо из критериев не будет соблюден:

for(i in 1:1000)
{
  test <- constrained_sample(mymins, mymaxs)
  if(!all(test > mymins) | !all(test < mymaxs) | abs(sum(test) - 1) > 1e6) cat("failure")
}

Это не вызывает ошибок, так как критерии всегда встречал Однако, как указывает @GregorThomas, в этом случае границы не являются реалистичными c. Мы можем увидеть ряд решений, ограниченных вашими условиями, используя блокпост:

samp <- constrained_sample(mymins, mymaxs)
for(i in 1:999) samp <- rbind(samp, constrained_sample(mymins, mymaxs))
df <- data.frame(val = c(samp[,1], samp[,2], samp[,3], samp[,4], samp[,5]), 
                 index = factor(rep(1:5, each = 1000)))
ggplot(df, aes(x = index, y = val)) + geom_boxplot()

Answer 2

Если ваши примерные границы реалистичны c, мы можем немного их уточнить, сузив диапазон возможностей. Для текущей версии вопрос с:

mymins = c(0.3, 0.1, 0, 0.2, 0.05)
mymaxs = c(0.9, 1, 1, 1, 0.7)

Какое максимальное значение для x[1]? Что ж, если x[2:5] примет минимальные значения, они составят в сумме 0.1 + 0 + 0.2 + 0.05 = 0.35, поэтому, исходя из только других минут , мы знаем, что максимальное значение для x[1] равно 1 - 0.35 = 0.65. 0.9 в mymaxs слишком велико.

Мы можем вычислить фактические максимальные значения, взяв минимум максимальных значений на основе минимумов и вектора mymaxs:

new_max = pmin(mymaxs, 1 - (sum(mymins) - mymins))
new_max
# [1] 0.65 0.45 0.35 0.55 0.40

Мы можем аналогичным образом пересмотреть минимальные границы, хотя в этом случае даже пересмотренные максимальные границы new_max достаточно высоки, чтобы это оказало какое-либо влияние на минимумы.

new_min = pmax(mymins, 1 - (sum(new_max) - new_max))
new_min
# [1] 0.30 0.10 0.00 0.20 0.05

С этими корректировками мы должен легко видеть, возможны ли какие-либо решения (all(new_min < new_max)). И затем генерация случайных чисел, как в ответе r2evans, должна go намного быстрее с использованием новых границ.

Answer 3

Поскольку вам нужно 5 случайных чисел для суммирования с 1, у вас действительно есть только 4 независимых числа и одно зависимое число.

mymins <- c(0.3, 0.1, 0, 0.2, 0.05)
mymaxs <- c(0.9, 1, 1, 1, 0.7)

set.seed(42)
iter <- 1000
while(iter > 0 &&
        (
          (1 - sum(x <- runif(4, mymins[-5], mymaxs[-5]))) < mymins[5] ||
            (1 - sum(x)) > mymaxs[5]
        )
      ) iter <- iter - 1
if (iter < 1) {
  # failed
  stop("unable to find something within 1000 iterations")
} else {
  x <- c(x, 1-sum(x))
}

sum(x)
# [1] 1
all(mymins <= x & x <= mymaxs)
# [1] TRUE
x
# [1] 0.37732330 0.21618036 0.07225311 0.24250359 0.09173965

Причина, по которой я использую iter, заключается в том, чтобы убедиться, что вы не тратить "бесконечное" количество времени, чтобы найти что-то. Если ваша комбинация mymins и mymaxs делает это математически неосуществимым (как в первом примере), то вам не нужно вращаться вечно. Если математически маловероятно найти его за разумное время, вам нужно взвесить, сколько времени вы хотите сделать это.

Одна из причин, по которой это занимает так много времени, заключается в том, что мы многократно увеличиваем энтропию. Если вы ожидаете, что это значение будет go в течение длительного времени, тогда, как правило, лучше предварительно рассчитать столько, сколько вы считаете нужным (в целом), и запустить все как матрицу.

set.seed(42)
n <- 10000
m <- matrix(runif(prod(n, length(mymins)-1)), nrow = n)
m <- t(t(m) * (mymaxs[-5] - mymins[-5]) + mymins[-5])
remainders <- (1 - rowSums(m))
ind <- mymins[5] <= remainders & remainders <= mymaxs[5]
table(ind)
# ind
# FALSE  TRUE 
#  9981    19 
m <- cbind(m[ind,,drop=FALSE], remainders[ind])
nrow(m)
# [1] 19
rowSums(m)
#  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
head(m)
#           [,1]      [,2]       [,3]      [,4]       [,5]
# [1,] 0.3405821 0.1306152 0.05931363 0.2199362 0.24955282
# [2,] 0.3601376 0.1367465 0.20235704 0.2477507 0.05300821
# [3,] 0.4469526 0.1279795 0.02265618 0.2881733 0.11423845
# [4,] 0.5450527 0.1029903 0.07503371 0.2052423 0.07168103
# [5,] 0.3161519 0.1469783 0.15290720 0.3268470 0.05711557
# [6,] 0.4782448 0.1185735 0.01664063 0.2178225 0.16871845
all(
  mymins[1] <= m[,1] & m[,1] <= mymaxs[1],
  mymins[2] <= m[,2] & m[,2] <= mymaxs[2],
  mymins[3] <= m[,3] & m[,3] <= mymaxs[3],
  mymins[4] <= m[,4] & m[,4] <= mymaxs[4],
  mymins[5] <= m[,5] & m[,5] <= mymaxs[5]
)
# [1] TRUE

На этот раз потребовалось 10000 попыток собрать 19 правильных комбинаций. Это может занять больше или меньше попыток, основанных на случайности, так что, ymmv, в отношении того, сколько вам нужно предварительно сгенерировать.