Среднее вычисление матрицы в моделировании Монте-Карло в R

В следующем коде я создаю имитацию падения снаряда методом Монте-Карло. Я вполне уверен, что симуляционная часть кода верна, но у меня возникают проблемы с вычислением среднего значения переменной g_estimate (гравитация) в симуляции. Я вычисляю g_estimate на каждой итерации, но когда я беру mean, он возвращает то же значение, что и g_estimate. Я думаю, это потому, что он вычисляет mean для каждой итерации, когда мне нужно общее mean. Я хочу создать матрицу или вектор для хранения всех значений g_estimate, но пока мне это не удалось. Вот код:

#Monte Carlo simulation of distance formula equation
set.seed(1)
N <- 100000
g = 9.8
h0 = 56.67
v0 = 0
n = 25 #Number of devisions for time variable
tt = seq(0, 3.4, len = n) #Split time tt into 25 values between 0 and 3.4
y = h0 + v0 * tt - 0.5 * g * tt ^ 2 + rnorm(n, sd = 1) #Traces the path of the projectile
X = cbind(1, tt, tt ^ 2) #Define the matrix of variables
A = solve(crossprod(X)) %*% t(X) # A = (X^T X)^-1 X^T

replicate(N, {
    y = h0 + v0 * tt - 0.5 * g * tt ^ 2 + rnorm(n, sd = 1) #Traces the path of the projectile
    Exercise_Beta <- A %*% y #Returns a matrix of beta values
    g_estimate <- c(-2 * Exercise_Beta[3]) #Calculates g constant each iteration
    return(mean(g_estimate))
})

Любая помощь приветствуется! Пожалуйста, дайте мне знать, если вам понадобится дополнительная информация. Заранее спасибо!


r mean matrix montecarlo
person William    schedule 04.09.2019    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Я переместил replicate... в переменную g_estimate, и этот код работает:

set.seed(1)
N <- 100000
g = 9.8
h0 = 56.67
v0 = 0
n = 25 #Number of devisions for time variable
tt = seq(0, 3.4, len = n) #Split time tt into 25 values between 0 and 3.4
y = h0 + v0 * tt - 0.5 * g * tt ^ 2 + rnorm(n, sd = 1) #Traces the path of the projectile
X = cbind(1, tt, tt ^ 2) #Define the matrix of variables
A = solve(crossprod(X)) %*% t(X) # A = (X^T X)^-1 X^T
g_estimate <- matrix(replicate(N, {
    y = h0 + v0 * tt - 0.5 * g * tt ^ 2 + rnorm(n, sd = 1) #Traces the path of the projectile
    Exercise_Beta <- A %*% y #Returns a matrix of beta values
    return(-2 * Exercise_Beta[3])
    }))

mean(g_estimate)
person William    schedule 04.09.2019