Как создать образцы из модели MVN?

Question

Я пытаюсь запустить некоторый код на R на основе этой статьи здесь в примере 5.1. Я хочу смоделировать следующее:

Мой фон на R не очень хороший, поэтому у меня есть следующий код ниже, как я могу сгенерировать гистограмму а образцы из этого?

xseq<-seq(0, 100, 1)  
n<-100
Z<- pnorm(xseq,0,1)
U<- pbern(xseq, 0.4, lower.tail = TRUE, log.p = FALSE)
Beta <- (-1)^U*(4*log(n)/(sqrt(n)) + abs(Z))

Answer 1

Некоторые демонстрации инструментов, которые будут полезны:

rnorm(1)                 # generates one standard normal variable
rnorm(10)                # generates 10 standard normal variables
rnorm(1, 5, 6)           # generates 1 normal variable with mu = 5, sigma = 6
                         # not needed for this problem, but perhaps worth saying anyway

rbinom(5, 1, 0.4)        # generates 5 Bernoulli variables that are 1 w/ prob. 0.4

Итак, чтобы сгенерировать один экземпляр беты:

n <- 100                 # using the value you gave; I have no idea what n means here
u <- rbinom(1, 1, 0.4)   # make one Bernoulli variable
z <- rnorm(1)            # make one standard normal variable
beta <- (-1)^u * (4 * log(n) / sqrt(n) + abs(z))

Но теперь вы хотели бы сделать это много раз для моделирования Монте-Карло. Один из способов сделать это - создать функцию с выводом beta и использовать функцию replicate(), например так:

n <- 100                    # putting this here because I assume it doesn't change
genbeta <- function(){      # output of this function will be one copy of beta
  u <- rbinom(1, 1, 0.4)
  z <- rnorm(1)
  return((-1)^u * (4 * log(n) / sqrt(n) + abs(z)))
}

# note that we don't need to store beta anywhere directly; 
# rather, it is just the return()ed value of the function we defined

betadraws <- replicate(5000, genbeta())
hist(betadraws)

Это даст эффект создания 5000 копий ваша бета-переменная и положить их в гистограмму.

Есть и другие способы сделать это - например, можно просто создать большую матрицу из случайных переменных и работать с ней напрямую - но я подумал, что это будет самый ясный подход для начала.

---

РЕДАКТИРОВАТЬ : я понял, что полностью проигнорировал второе уравнение, которое вы, вероятно, не хотели.

Теперь мы создали вектор beta значения, и вы можете контролировать длину вектора в первом параметре функции replicate() выше. Я оставлю это как 5000 в моем следующем примере ниже.

Чтобы получить случайные выборки вектора Y, вы можете использовать что-то вроде:

x <- replicate(5000, rnorm(17))     
  # makes a 17 x 5000 matrix of independent standard normal variables
epsilon <- rnorm(17)
  # vector of 17 standard normals
y <- x %*% betadraws + epsilon
  # y is now a 17 x 1 matrix (morally equivalent to a vector of length 17)

, и если вы хотите получить много из них вы можете обернуть , что , в другую функцию и replicate() it.

В качестве альтернативы, если вам не нужен вектор Y, а только один компонент Y_i:

x <- rnorm(5000)    
  # x is a vector of 5000 iid standard normal variables 
epsilon <- rnorm(1)
  # epsilon_i is a single standard normal variable
y <- t(x) %*% betadraws + epsilon
  # t() is the transpose function; y is now a 1 x 1 matrix