Квадратура для аппроксимации преобразованного бета-распределения в R

Я использую R для запуска моделирования, в котором я использую тест отношения правдоподобия для сравнения двух моделей ответа вложенных элементов. Одна версия LRT использует совместную функцию правдоподобия L (θ, ρ), а другая - функцию маргинального правдоподобия L (ρ). Я хочу проинтегрировать L (θ, ρ) по f (θ), чтобы получить маргинальное правдоподобие L (ρ). У меня есть два условия: в одном f (θ) является стандартной нормой (μ = 0, σ = 1), и я понимаю, что я могу просто выбрать количество точек абсцисс, скажем 20 или 30, и использовать Гаусса-Эрмита. квадратурная аппроксимация этой плотности. Но в другом условии f (θ) является линейно преобразованным бета-распределением (a = 1,25, b = 10), где линейное преобразование B '= 11,14 * (B-0,11) таково, что B' также имеет (приблизительно) μ = 0, σ = 1.

Я достаточно запутался в том, как реализовать квадратуру для бета-распределения, но линейное преобразование сбивает меня с толку еще больше. У меня тройной вопрос: (1) могу ли я использовать некоторые вариации квадратуры для аппроксимации f (θ), когда θ распределяется как это линейно преобразованное бета-распределение, (2) как бы я реализовал это в R, и (3) это нелепая трата времени на то, что существует явно более быстрый и лучший способ выполнить эту задачу? (Я попытался написать свою собственную функцию численного приближения, но обнаружил, что моя реализация, ограниченная языком R, была слишком медленной, чтобы ее хватило.)

Спасибо!


r numerical-integration
person psychometriko    schedule 04.06.2013    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Во-первых, я предполагаю, что вы можете выразить свои 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")

Изменить: не уверен, что вы также спрашиваете, как определить f (θ) для преобразованного бета-распределения; это кажется довольно элементарным для человека, работающего с совместной и предельной вероятностью. Но если да, то плотность B '= a * B + b, заданная f (B), равна

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

Итак, в вашем случае f(theta) это dbeta(theta/11.14 + 0.11, 1.25, 10) / 11.14

person Hong Ooi    schedule 04.06.2013
comment
Спасибо, это очень информативно! У меня действительно есть L (θ, ρ) и f (θ), выраженные в терминах реального кода. Я получил все необходимые результаты на основе L (θ, ρ), и мне просто было трудно понять, как перейти от L (θ, ρ) к L (ρ). Так что мне даже не нужно идти по пути квадратуры; это хорошо знать! Я подожду, чтобы увидеть, почувствует ли кто-нибудь другое, а если нет, то отмечу это как принятый ответ. - person psychometriko; 04.06.2013