0

我正在使用 R 运行模拟，其中我使用似然比检验来比较两个嵌套项目响应模型。LRT 的一个版本使用联合似然函数 L(θ,ρ)，而另一个版本使用边际似然函数 L(ρ)。我想在 f(θ) 上积分 L(θ,ρ) 以获得边际似然 L(ρ)。我有两个条件：其中一个，f(θ) 是标准正态（μ=0，σ=1），我的理解是我可以选择一些横坐标点，比如 20 或 30，并使用 Gauss-Hermite正交来近似这个密度。但在另一种情况下，f(θ) 是线性变换的 beta 分布 (a=1.25,b=10)，其中线性变换 B'=11.14*(B-0.11) 使得 B' 也具有（大约） μ=0，σ=1。

我对如何为 beta 分布实现正交感到很困惑，但是线性变换让我更加困惑。我的问题有三个：（1）当 θ 分布为这种线性变换的 beta 分布时，我可以使用正交的一些变化来近似 f(θ)，（2）我将如何在 R 中实现它，以及（3）这是一个荒谬的时间浪费，以至于有明显更快更好的方法来完成这项任务？（我尝试编写自己的数值逼近函数，但发现我的实现仅限于 R 语言，速度太慢而无法满足。）

谢谢！

4

1 回答 1

1

首先，我假设你可以用实际代码来表达你的 L(θ,ρ) 和 f(θ)；否则你有点搞砸了。鉴于该假设，您可以使用它integrate来执行必要的计算。这样的事情应该让你开始；只需插入 L 和 f 的表达式。

marglik <- function(rho) {
    integrand <- function(theta, rho) L(theta, rho) * f(theta)
    # set your lower/upper integration limits as appropriate
    integrate(integrand, lower=-5, upper=5, rho=rho)
}

为此，您的被积函数必须被矢量化；即，给定一个向量输入theta，它必须返回一个输出向量。如果您的代码不符合要求，您可以Vectorize在将其传递给之前在被积函数上使用integrate：

integrand <- Vectorize(integrand, "theta")

编辑：不确定您是否还询问如何为转换后的 beta 分布定义 f(θ)；对于处理联合和边际可能性的人来说，这似乎相当基本。但如果你是，那么在给定 f(B) 的情况下，B' = a*B + b 的密度是

f'(B') = f(B)/a = f((B' - b)/a) / a

所以在你的情况下，f(theta)是dbeta(theta/11.14 + 0.11, 1.25, 10) / 11.14

于 2013-06-04T12:52:28.727 回答