Генерация случайных чисел в R, удовлетворяющих ограничениям - PullRequest
Новая
       18

Генерация случайных чисел в R, удовлетворяющих ограничениям

2 голосов
/ 21 апреля 2020

Мне нужна помощь с кодом для генерации случайных чисел в соответствии с ограничениями. В частности, я пытаюсь смоделировать случайные числа ALFA и BETA, соответственно, из нормального и гамма-распределения, чтобы ALFA - BETA <1. </p>

Вот то, что я написал, но оно не работает вообще. 

set.seed(42) 
n <- 0 
repeat {
  n <- n + 1
  a <- rnorm(1, 10, 2)
  b <- rgamma(1, 8, 1)
  d <- a - b
  if (d < 1) 
  alfa[n] <- a
  beta[n] <- b
  l = length(alfa)
  if (l == 10000) break
}

Ответы [ 2 ]

3 голосов
/ 21 апреля 2020

Из-за векторизации будет быстрее генерировать числа «все сразу», а не как все oop: 

set.seed(42)
N = 1e5
a = rnorm(N, 10, 2)
b = rgamma(N, 8, 1)
d = a - b
alfa = a[d < 1]
beta = b[d < 1]
length(alfa)
# [1] 36436

В результате получено 100 000 кандидатов, 36 436 из которых соответствуют вашим критериям. Если вы хотите сгенерировать n семплов, попробуйте установить N = 4 * n, и вы, вероятно, сгенерируете более чем достаточно, оставьте первый n.

У вашего l oop есть 2 проблемы: (а) вам нужны фигурные скобки, чтобы заключить несколько строк после оператора if. (б) вы используете n в качестве счетчика попыток, но он должен быть счетчиком успеха. Как написано, ваш l oop остановится только в случае успеха 10000-й попытки. Переместите n <- n + 1 внутри оператора if, чтобы исправить: 

set.seed(42) 
n <- 0 
alfa = numeric(0)
beta = numeric(0)
repeat {
  a <- rnorm(1, 10, 2)
  b <- rgamma(1, 8, 1)
  d <- a - b
  if (d < 1) {
    n <- n + 1
    alfa[n] <- a
    beta[n] <- b
    l = length(alfa)
    if (l == 500) break
  }
}

Но первый путь лучше ... из-за "роста" alfa и beta в l oop, и генерируя числа по одному, этот метод генерирует 500 чисел дольше, чем приведенный выше код для генерации 30000.

0 голосов
/ 21 апреля 2020

Как прокомментировал @Gregor Thomas, ваша попытка провалилась из-за отсутствия фигурных скобок для включения оператора if. Если вы хотите пропустить {} для if управления, возможно, вы можете попробовать код ниже 

set.seed(42) 
r <- list()
repeat {
  a <- rnorm(1, 10, 2)
  b <- rgamma(1, 8, 1)
  d <- a - b
  if (d < 1) r[[length(r)+1]] <- cbind(alfa = a, beta = b)
  if (length(r) == 100000) break
}
r <- do.call(rbind,r)

, такой что 

> head(r)
          alfa      beta
[1,]  9.787751 12.210648
[2,]  9.810682 14.046190
[3,]  9.874572 11.499204
[4,]  6.473674  8.812951
[5,]  8.720010  8.799160
[6,] 11.409675 10.602608
...