0

我想测试我使用 R 中的 stats4::mle 拟合的两个非嵌套模型中的哪一个使用 Vuong 和 Clarke 测试提供了更好的拟合。

这是我正在拟合的数据的一小部分，两个不同的模型（函数“w”不同）和相应的 mle()：

library(stats4)
### Data 
z1 <- c(0.1111111, 0.1037037, 0.1222222, 0.1111111, 0.1074074, 0.1666667, 0.1333333, 0.2000000, 0.1333333, 0.1074074,
        0.1037037, 0.1111111, 0.1333333, 0.2000000, 0.1222222, 0.1111111, 0.1666667, 0.1333333, 0.1111111, 0.1333333,
        0.1111111, 0.1666667, 0.1074074, 0.1333333, 0.1222222, 0.2000000, 0.1037037)

z2 <- c(0.08888889, 0.06666667, 0.07777778, 0.00000000, 0.03333333, 0.09259259, 0.09629630, 0.08888889, 0.06666667,
        0.03333333, 0.06666667, 0.08888889, 0.06666667, 0.08888889, 0.07777778, 0.00000000, 0.09259259, 0.09629630,
        0.00000000, 0.09629630, 0.08888889, 0.09259259, 0.03333333, 0.06666667, 0.07777778, 0.08888889, 0.06666667)

p <-  c(0.5, 0.9, 0.5, 0.9, 0.9, 0.1, 0.1, 0.1, 0.5, 0.9, 0.9, 0.5, 0.5, 0.1, 0.5, 0.9, 0.1, 0.1, 0.9, 0.1, 0.5, 0.1, 0.9, 0.5, 0.5, 0.1, 0.9)

zce <- c(0.11055556, 0.10277778, 0.11000000, 0.10833333, 0.10185185, 0.11666667, 0.13240741, 0.14166667, 0.13166667,
         0.07222222, 0.08796296, 0.09944444, 0.09500000,0.10833333, 0.09444444, 0.05277778, 0.10925926, 0.11759259,
         0.05833333, 0.10277778, 0.09277778, 0.10925926, 0.06111111, 0.08833333, 0.09222222, 0.12500000, 0.09166667)


### Functions:
u <- function(x,n) 
{
  ifelse(n!=1,util <- x^(1-n)/(1-n), util <- log(x))
  return(util)
}
u.inv <- function(x,n)
{
  ifelse(n !=1, inv.util <- ((1-n)*(x))^(1/(1-n)), inv.util <- exp(x))
  return(inv.util)
}

v = function(x,n){return(1/(u(maxz,n)-u(minz,n))*(u(x,n)-u(minz,n)))}
v.inv = function(x,n){return(u.inv(x*(u(maxz,n)-u(minz,n))+u(minz,n),n))}

maxz = 135
minz = 0


### model 1
w <- function(p,a,b){return(exp(-b*(-log(p))^(1-a)))}

### Loglikelihood 1
LL1 <- function(n,a,b,s)
{
  V = (v(z1,n)-v(z2,n))*w(p,a,b) + v(z2,n) 
  res = zce - v.inv(V,n)
  ll = dnorm(res, 0, s,log=T)
  ll.fin1 <<- ll ### record ll per datapoint given optimal parameters
  return(-sum(ll))
}

### mle 1
fit.model1 <- mle(LL1,
           start = list(n = 0.1,a=0.1,b=0.1,s=0.1),
           method = "L-BFGS-B",
           lower = list(n=-Inf,a = -Inf, b = 0.0001, s=0.0001),
           upper = list(n=0.9999,a = 0.9999, b = Inf, s=Inf),
           control = list(maxit = 500, ndeps = rep(0.000001,4)),
           nobs=length(z1))


######################

### model 2
w <- function(p,a,b){return((b*p^a)/(b*p^a+(1-p)^a))}

### Loglikelihood 2
LL2 <- function(n,a,b,s)
{
  V = (v(z1,n)-v(z2,n))*w(p,a,b) + v(z2,n) 
  res = zce - v.inv(V,n)
  ll = dnorm(res, 0, s,log=T)
  ll.fin2 <<- ll ### record ll per datapoint given optimal parameters
  return(-sum(ll))
}

### mle 2
fit.model2 <- mle(LL2,
                  start = list(n = 0.1,a=0.1,b=0.1,s=0.1),
                  method = "L-BFGS-B",
                  lower = list(n=-Inf,a = 0.0001, b = 0.0001, s=0.0001),
                  upper = list(n=0.9999,a = Inf, b = Inf, s=Inf),
                  control = list(maxit = 500, ndeps = rep(0.000001,4)),
                  nobs=length(z1))
4

1 回答 1

1

在“游戏”包作者的帮助下，其中包括最后通牒游戏的 Vuong 和 Clarke 测试，我能够为我的 MLE 优化编写测试代码。

重要的是，应该在最佳参数集记录每个数据点的对数似然，这可以通过将“ll”保存在全局变量中来完成。该向量在优化过程的最终迭代中所取的值正是所需的向量，可以通过对其求和并将其与 logLike(fit) 进行比较来轻松检查。

有了拟合和个体密度，就可以定义 Vuong 和 Clarke 检验。

### Define tests
vuong.test <- function(fit1,fit2,ldens1,ldens2,alternative) 
{
  if(nobs(fit1) != nobs(fit2)){return("Error: number of observations used to fit the models unequal")}
  n.obs <- nobs(fit1)
  n.par1 <- nrow(coef(summary(fit1)))
  n.par2 <- nrow(coef(summary(fit2)))
  
  LR <- sum(ldens1) - sum(ldens2) - (n.par1-n.par2)/2*log(n.obs)
  LDdif <- sqrt(sum((ldens1 - ldens2)^2)/n.obs)
  z = LR/(sqrt(n.obs)*LDdif)
  
  alternative <- substr(alternative,1,1) ## alternative specifies whether you test whether fit1 is lesser, greater or equal to fit2.
  if(alternative == "l"){pval = 1-pnorm(z,mean = 0,sd = 1)}
  if(alternative == "g"){pval = pnorm(z,mean = 0,sd = 1)}
  if(alternative == "t"){pval = 2*pnorm(-abs(z),mean = 0,sd = 1)}
  
  return(pval)
}

clarke.test <- function(fit1,fit2,ldens1,ldens2,alternative)
{
  if(nobs(fit1) != nobs(fit2)){return("Error: number of observations used to fit the models unequal")}
  n.obs <- nobs(fit1)
  n.par1 <- nrow(coef(summary(fit1)))
  n.par2 <- nrow(coef(summary(fit2)))
  
  correction <- (n.par1-n.par2)*(log(n.obs)/(2*n.obs))
  z = sum(ldens1 - ldens2 > correction)
  
  alternative <- substr(alternative,1,1) ## alternative specifies whether you test whether fit1 is lesser, greater or equal to fit2.
  if(alternative == "l"){pval = 1-pbinom(q = z,size=n.obs,prob = 0.5)}
  if(alternative == "g"){pval = pbinom(q = z,size=n.obs,prob = 0.5)}
  if(alternative == "t"){zz <- min(z,n.obs-z); pval = 2*pbinom(q = zz,size=n.obs,prob = 0.5)}
  
  return(pval)
}

在哪里：

  • fit1 = 使用 stats4::mle 保存的第一个模型的拟合
  • fit2 = 使用 stats4::mle 保存的第二个模型的拟合
  • ldens1 = 第一个模型的最优参数集的对数似然向量。
  • ldens2 = 第二个模型的最佳参数集的对数似然向量。
  • 替代方案 = 用“t”/“l”/“g”测试是否可以拒绝模型 1 与模型 2 相比是等距/更远/更接近真实数据生成过程。

给定第三种模型：

##### Model 3
w3 <- function(p,a,b){return(p)}
w <- w3

### LL3
LL3 <- function(n,s)
{
  V = (v(z1,n)-v(z2,n))*w(p,a,b) + v(z2,n) 
  res = zce - v.inv(V,n)
  ll = dnorm(res, 0, s,log=T)
  ll.fin3 <<- ll
  -sum(ll)
}

### mle 3
fit.model3 <- mle(LL3,
                  start = list(n = 0.1,s=0.1),
                  method = "L-BFGS-B",
                  lower = list(n=-Inf,s=0.0001),
                  upper = list(n=0.9999, s=Inf),
                  control = list(maxit = 500, ndeps = rep(0.000001,2)),
                  nobs=length(z1))

将拒绝基于 Clarke 的等距离：

> clarke.test(fit1 = fit.model1, ldens1 = ll.fin1,
+             fit2 = fit.model3, ldens2 = ll.fin3,
+             alternative = "t")
[1] 0.01915729

游戏包：https ://github.com/brentonk/games/blob/master/games/R/nonnest.r

于 2020-03-11T12:02:03.173 回答