Я пытаюсь вычислить логарифмическую вероятность для обобщенной нелинейной регрессии наименьших квадратов для функции оптимизированной с помощью функция в пакете R , используя ковариационную матрицу дисперсии, генерируемую расстояниями на филогенетическом дереве, предполагающем броуновское движение ( из пакета). Следующий воспроизводимый код R подходит для модели gnls с использованием данных x, y и случайного дерева с 9 таксонами:gnls
nlme
corBrownian(phy=tree)
ape
require(ape)
require(nlme)
require(expm)
tree <- rtree(9)
x <- c(0,14.51,32.9,44.41,86.18,136.28,178.21,262.3,521.94)
y <- c(100,93.69,82.09,62.24,32.71,48.4,35.98,15.73,9.71)
data <- data.frame(x,y,row.names=tree$tip.label)
model <- y~beta1/((1+(x/beta2))^beta3)
f=function(beta,x) beta[1]/((1+(x/beta[2]))^beta[3])
start <- c(beta1=103.651004,beta2=119.55067,beta3=1.370105)
correlation <- corBrownian(phy=tree)
fit <- gnls(model=model,data=data,start=start,correlation=correlation)
logLik(fit)
Я хотел бы рассчитать логарифмическую вероятность "вручную" (в R, но без использования logLik
функции) на основе оценочных параметров, полученных из, gnls
чтобы он соответствовал выходным данным из logLik(fit)
. ПРИМЕЧАНИЕ: я не пытаюсь оценить параметры; Я просто хочу вычислить логарифмическую вероятность параметров, оцениваемых gnls
функцией (хотя, если у кого-то есть воспроизводимый пример того, как оценивать параметры без него gnls
, мне было бы очень интересно увидеть его!).
Я не совсем уверен, как сделать это в R. Линейная алгебраическая запись, описанная в моделях со смешанными эффектами в S и S-Plus (Pinheiro и Bates), очень важна для меня, и ни одна из моих попыток не совпала logLik(fit)
. Вот подробности, описанные Пинейро и Бейтсом:
Логарифмическая вероятность для обобщенной нелинейной модели наименьших квадратов где вычисляется следующим образом:
где - количество наблюдений, а .f ∗ i ( β ) = f ∗ i ( ϕ i , v i )
положительно определен, и
Для фиксированных и оценка ML равна
и профилированная логарифмическая вероятность
который используется с алгоритмом Гаусса-Зейделя для нахождения ML-оценок и . Используется менее смещенная оценка :
где представляет длину .
Я составил список конкретных вопросов, с которыми я сталкиваюсь:
- Что такое ? Это матрица расстояний, создаваемая in , или она должна быть каким-то образом преобразована или параметризована , или что-то еще целиком?
big_lambda <- vcv.phylo(tree)
ape
- Would В , или уравнение для менее пристрастной оценки (последнего уравнения в этой должности)?
fit$sigma^2
- Нужно ли использовать для вычисления логарифмического правдоподобия или это просто промежуточный шаг для оценки параметров? Кроме того, как используется ? Это одно значение или вектор, и оно умножается на все или просто недиагональные элементы и т. Д.?
- Что такое? Это будет в пакете ? Если это так, я запутался в том, как рассчитать сумму , потому что возвращает одно значение, а не вектор.M ∑ i = 1 | | y ∗ i - f ∗ i ( β ) | | 2
norm(y-f(fit$coefficients,x),"F")
Matrix
norm()
- Как рассчитать? Является ли это , где находится , или это из пакета ? Если да , то как взять сумму матрицы (или подразумевается, что это только диагональные элементы)?Λ я
log(diag(abs(big_lambda)))
big_lambda
logm(abs(big_lambda))
expm
logm()
- Просто чтобы подтвердить, является рассчитывается следующим образом: ?
t(solve(sqrtm(big_lambda)))
- Как рассчитываются и ? Это одно из следующего: f ∗ i ( β )
y_star <- t(solve(sqrtm(big_lambda))) %*% y
и
f_star <- t(solve(sqrtm(big_lambda))) %*% f(fit$coefficients,x)
или это будет
y_star <- t(solve(sqrtm(big_lambda))) * y
и
f_star <- t(solve(sqrtm(big_lambda))) * f(fit$coefficients,x)
?
Если ответить на все эти вопросы, теоретически, я думаю, что логарифмическая вероятность должна быть вычисляемой, чтобы соответствовать выходным данным logLik(fit)
. Любая помощь по любому из этих вопросов будет принята с благодарностью. Если что-то нуждается в разъяснении, пожалуйста, дайте мне знать. Благодарность!
ОБНОВЛЕНИЕ : я экспериментировал с различными возможностями для вычисления логарифмической вероятности, и вот лучшее, что я придумал до сих пор. logLik_calc
последовательно от 1 до 3 от значения, возвращаемого logLik(fit)
. Либо я близок к реальному решению, либо это чисто случайно. Есть предположения?
C <- vcv.phylo(tree) # variance-covariance matrix
tC <- t(solve(sqrtm(C))) # C^(-T/2)
log_C <- log(diag(abs(C))) # log|C|
N <- length(y)
y_star <- tC%*%y
f_star <- tC%*%f(fit$coefficients,x)
dif <- y_star-f_star
sigma_squared <- sum(abs(y_star-f_star)^2)/N
# using fit$sigma^2 also produces a slightly different answer than logLik(fit)
logLik_calc <- -((N*log(2*pi*(sigma_squared)))+
sum(((abs(dif)^2)/(sigma_squared))+log_C))/2
Ответы:
Давайте начнем с более простого случая, когда нет корреляционной структуры для остатков:
Вероятность логарифма может быть легко вычислена вручную с помощью:
Поскольку остатки независимы, мы можем просто использовать,
dnorm(..., log=TRUE)
чтобы получить отдельные логарифмические вероятностные термины (а затем суммировать их). В качестве альтернативы мы могли бы использовать:Обратите внимание, чтоσ2
fit$sigma
это не «менее предвзятая оценка » - поэтому нам нужно сначала выполнить коррекцию вручную.Теперь для более сложного случая, когда остатки коррелируют:
Здесь нам нужно использовать многомерное нормальное распределение. Я уверен, что где-то есть функция для этого, но давайте просто сделаем это вручную:
источник
vcv
функцией) - но вам нужно получить матрицу корреляции и затем использовать чтобы превратить ее в матрицу var-cov.