R中蒙特卡罗模拟中矩阵的平均计算答案

【问题标题】：Mean Calculation of a matrix within a Monte Carlo simulation in RR中蒙特卡罗模拟中矩阵的平均计算
【发布时间】：2020-01-07 12:12:56
【问题描述】：

在下面的代码中，我正在创建一个弹丸下落的蒙特卡罗模拟。我相当有信心代码的模拟部分是正确的，但是我在计算模拟中变量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 matrix mean montecarlo

【解决方案1】：

我将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)

【讨论】：