Тест Вуонга и Кларка, использующий вывод mle в R? - PullRequest
Новая
       91

Тест Вуонга и Кларка, использующий вывод mle в R?

0 голосов
/ 03 марта 2020

Я хочу проверить, какая из двух не вложенных моделей, которые мне подходят, используя stats4 :: mle в R, обеспечивает лучшее соответствие с использованием теста Вуонга и Кларка.

Вот небольшая часть данных, которые я подгоняю, две разные модели (функция "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))

1 Ответ

0 голосов
/ 11 марта 2020

С помощью автора пакета "games", который включает в себя тесты Vuong и Clarke для ультимативных игр, я смог закодировать тесты для своей оптимизации MLE.

Важно отметить, что следует записать логарифмические вероятности каждой точки данных при оптимальном наборе параметров, что можно сделать, сохранив «ll» в глобальной переменной. Значения, которые этот вектор принимает на последней итерации процесса оптимизации, являются именно тем необходимым вектором, который можно легко проверить, суммируя его и сравнивая его с logLike (fit).

С подгонками и индивидуальными плотностями Vuong и Затем можно определить критерий Кларка. 

### 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 = pnorm(z,mean = 0,sd = 1)}
  if(alternative == "g"){pval = 1-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 = pbinom(q = z,size=n.obs,prob = 0.5)}
  if(alternative == "g"){pval = 1-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 = вектор логарифмических правдоподобий при оптимальном наборе параметров второй модели. - alternative = проверить с помощью «t» / «l» / «g», можно ли отклонить, что модель 1 равноудалена / дальше / ближе к истинному процессу генерирования данных, чем модель 2.

Например, отклонение не может быть отклонено Эквидистантность между моделью 1 и моделью 2 в соответствии с Vuong: 

> vuong.test(fit1 = fit.model1, ldens1 = ll.fin1,
+            fit2 = fit.model2, ldens2 = ll.fin2,
+            alternative = "t")
[1] 0.5842877

Учитывая третью модель: 

##### 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))

Например, не может отклонить эквидистантность на основе теста Vuong, но отклонить на основе Кларка: 

> vuong.test(fit1 = fit.model1, ldens1 = ll.fin1,
+             fit2 = fit.model3, ldens2 = ll.fin3,
+             alternative = "t")
[1] 0.1823776
> 
> 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

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.

Похожие темы

...