«оптимизировать», не находя переменные внутри вызова функции - PullRequest
2 голосов
/ 17 ноября 2010

У меня есть следующая функция, которая извлекает некоторые данные из распределения хи-квадрат и сравнивает распределение X с известным распределением хи-квадрат, используя максимальную вероятность.Эта процедура моделируется nSims раз.(Я сравниваю эти результаты с результатами теста перестановки, но этот код исключен.)

chi2c <- function(xdf=2, yObs=100, xObs=100, nSims=1000, nPerm=500, alpha=0.05){
  simResults <- sapply(1:nSims, function(x){
    # Draw variables
    x  <- rchisq(xObs, df=xdf)
    # Other variables not relevant here
    # [[snip]]

    # Permutation test
    # [[snip]]

    # Calculate the statistics necessary for maximum likelihood
    n       <<- length(x)
    sumlx   <<- sum(log(x))
    sumx    <<- sum(x)
    # Calculate the maximum likelihood estimate
    dfhat   <- optimize(f=c2ll, interval=c(1, 10), maximum=TRUE)$maximum
    # Calculate the test statistic: -2 times the log likelihood ratio
    llr     <- -2 * (c2ll(2) - c2ll(dfhat))
    # Compare the test statistic to its asymptotic dist: chi-squared
    lReject <- llr > qchisq(1 - alpha, df=1)

    # Provide the results
    # [[snip]]
  })
  # Calcuate means across simulations
  rowMeans(simResults)
}

Эта функция вызывает c2ll, функцию вероятности хи-квадрат

c2ll <- function(dfHat){
  -n * log(gamma(dfHat/2)) - n * (dfHat/2) * log(2) + 
  (dfHat/2 - 1) * sumlx - sumx/2
}

Эта функция делает то, что мне хотелось бы, и точна, но я не понимаю, почему мне нужно установить статистику максимального правдоподобия (n, sumlx и sumx) в глобальном масштабе, чтобы заставить ее работать;optimize не находит их, если я устанавливаю их только внутри функции, используя <-.Я попытался установить их внутри optimize, но это тоже не сработало.Спасибо за вашу помощь.

Чарли

Ответы [ 3 ]

2 голосов
/ 17 ноября 2010

R имеет лексическую область видимости, что означает, что функции ищут переменные в среде, в которой они были определены. c2ll определен в глобальной среде, поэтому он не видит ваши определения n, sumx и sum1x внутри функции. S, с другой стороны, использует динамическую область видимости, которая ведет себя как , которую вы ожидаете (то есть ищет переменные в области действия, в которой она была вызвана). Компьютерные ученые, как правило, считают, что динамическое определение области было тупиковой плохой идеей, и что лексическое определение области - это путь.

На практике, что вы можете с этим поделать?

Ну, есть пара вариантов ...

Сначала , вы можете определить свою функцию локально:

n       <- length(x)
sumlx   <- sum(log(x))
sumx    <- sum(x)
c2ll <- function(dfHat){
    -n * log(gamma(dfHat/2)) - n * (dfHat/2) * log(2) + (dfHat/2 - 1) * sumlx - sumx/2
}
dfhat   <- optimize(f=c2ll, interval=c(1, 10), maximum=TRUE)$maximum

Второй , вы можете сделать так, чтобы c2ll принял дополнительные параметры и передал эти параметры через оптимизацию.

#in global env
c2ll <- function(dfHat,n,sum1x,sumx){
  -n * log(gamma(dfHat/2)) - n * (dfHat/2) * log(2) + 
  (dfHat/2 - 1) * sumlx - sumx/2
}

#...
#in function
n       <- length(x)
sumlx   <- sum(log(x))
sumx    <- sum(x)
dfhat   <- optimize(f=c2ll, interval=c(1, 10), n=n,sum1x=sum1x,sumx=sumx, maximum=TRUE)$maximum

Оба варианта являются чистыми и сохраняют инкапсуляцию ваших функций.

1 голос
/ 17 ноября 2010

Ваша проблема проистекает из лексических правил определения области применения, которые использует R. Подробнее см. В руководстве language определений . Короче говоря, ваша функция c2ll ищет переменные в среде, в которой она была определена.

Чтобы избежать этой проблемы, вы должны явным образом передать n, sum1x и sum2x в качестве аргументов вашей функции или определить вашу функцию локально в функции ch2c напрямую.

Это довольно распространенный вопрос, на SO .

есть много интересных примеров.
1 голос
/ 17 ноября 2010

Ваша функция simResults возвращает логический вектор.Вместо того чтобы использовать rowMeans, просто используйте среднее значение (simResults), и результаты выглядят разумно разумными ... по крайней мере, в той степени, в которой:

> chi2c(alpha=0.05)
[1] 0.057
> chi2c(alpha=0.5)
[1] 0.503
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...