如何将约束添加到对数似然函数中?
How to add constraint into loglikelihood function?
我有一个时间序列模型 (INGARCH):
lambda_t = alpha0 + alpha1*(x_(t-1)) + beta1*(lambda_(t-1))
X_t ~ poisson (lambda_t)
其中 t 是观察或数据的长度,alpha0、alpha1 和 beta1 是参数。
X_t
是数据序列,lambda_t是均值序列。
该模型的条件为alpha1 + beta1 < 1
。
据我估计,我想在我的代码中加入alpha1 + beta1 <
1的条件,我在对数似然函数中加入了一个while循环,但是循环无法停止。
我该怎么做才能解决这个问题?有没有其他方法可以在不使用 while 循环的情况下添加 alpha1 + beta1 < 1
的约束?
下面是我的代码:
ll <- function(par) {
h.new = rep(0,n)
#par[1] is alpha0
#par[2] is alpha1
#par[3] is beta1
while(par[2] + par[3] < 1){
for (i in 2:n) {
h.new[i] <- par[1] + par[2] * dat[i-1] + par[3] * h.new[i-1]
}
-sum(dpois(dat, h.new, log=TRUE))
}
}
#simply generate a dataset as I have not found a suitable real dataset to fit in
set.seed(77)
n=400
dat <- rpois(n,36)
nlminb(start = c(0.1,0.1,0.1), lower = 1e-6, ll)
在此期间您根本不更改 par
。特别是,如果你在 while 中打印了 par[1]
和 par[2]
,你会看到你在无休止地打印原始值 0.1 - 因此你永远被困在 while
中。
par
是来自 nlminb
的每次调用中的单个不变对象。你只需要确定 par 是否不好,你 return 不是最小值,所以 nlminb
不会继续朝那个方向搜索:
ll <- function(par) {
#If alpha + beta > 1, this is terrible and return an infinite score
#It may be better to throw an error if you get NaN values! The if will
#fail anyway, but if you want to power through add checks:
if( is.nan(par[2]) || is.nan(par[3]) || par[2]+par[3]>1) return(Inf)
h.new = rep(0,n)
#remove while
for (i in 2:n) {
h.new[i] <- par[1] + par[2] * dat[i-1] + par[3] * h.new[i-1]
}
-sum(dpois(dat, h.new, log=TRUE))
}
算法nlminb
(或任何最小化函数)非常粗略:
- 将参数设置为初始猜测
- 将参数发送到 objective 函数
猜新参数:
一个。如果分数没有太大提高,return 最小化猜测
b。如果分数不错,继续往这个方向搜索
c。否则,在其他方向搜索
使用新参数返回(2)
请注意,您必须 return 每组参数的分数,您不要在 objective 函数中重复它们。
我有一个时间序列模型 (INGARCH):
lambda_t = alpha0 + alpha1*(x_(t-1)) + beta1*(lambda_(t-1))
X_t ~ poisson (lambda_t)
其中 t 是观察或数据的长度,alpha0、alpha1 和 beta1 是参数。
X_t
是数据序列,lambda_t是均值序列。
该模型的条件为alpha1 + beta1 < 1
。
据我估计,我想在我的代码中加入alpha1 + beta1 <
1的条件,我在对数似然函数中加入了一个while循环,但是循环无法停止。
我该怎么做才能解决这个问题?有没有其他方法可以在不使用 while 循环的情况下添加 alpha1 + beta1 < 1
的约束?
下面是我的代码:
ll <- function(par) {
h.new = rep(0,n)
#par[1] is alpha0
#par[2] is alpha1
#par[3] is beta1
while(par[2] + par[3] < 1){
for (i in 2:n) {
h.new[i] <- par[1] + par[2] * dat[i-1] + par[3] * h.new[i-1]
}
-sum(dpois(dat, h.new, log=TRUE))
}
}
#simply generate a dataset as I have not found a suitable real dataset to fit in
set.seed(77)
n=400
dat <- rpois(n,36)
nlminb(start = c(0.1,0.1,0.1), lower = 1e-6, ll)
在此期间您根本不更改 par
。特别是,如果你在 while 中打印了 par[1]
和 par[2]
,你会看到你在无休止地打印原始值 0.1 - 因此你永远被困在 while
中。
par
是来自 nlminb
的每次调用中的单个不变对象。你只需要确定 par 是否不好,你 return 不是最小值,所以 nlminb
不会继续朝那个方向搜索:
ll <- function(par) {
#If alpha + beta > 1, this is terrible and return an infinite score
#It may be better to throw an error if you get NaN values! The if will
#fail anyway, but if you want to power through add checks:
if( is.nan(par[2]) || is.nan(par[3]) || par[2]+par[3]>1) return(Inf)
h.new = rep(0,n)
#remove while
for (i in 2:n) {
h.new[i] <- par[1] + par[2] * dat[i-1] + par[3] * h.new[i-1]
}
-sum(dpois(dat, h.new, log=TRUE))
}
算法nlminb
(或任何最小化函数)非常粗略:
- 将参数设置为初始猜测
- 将参数发送到 objective 函数
猜新参数:
一个。如果分数没有太大提高,return 最小化猜测
b。如果分数不错,继续往这个方向搜索
c。否则,在其他方向搜索
使用新参数返回(2)
请注意,您必须 return 每组参数的分数,您不要在 objective 函数中重复它们。