Может ли эмпирический гессиан М-оценки быть неопределенным?

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

Это кажется мне неправильным, поскольку (помимо численных проблем) гессиан должен быть положительно полуопределенным в результате определения М-оценки как значения параметра, который минимизирует целевую функцию для данного образца, и общеизвестный факт, что в (локальном) минимуме гессиан положительно полуопределен.

Прав ли мой аргумент?

[РЕДАКТИРОВАТЬ: заявление было удалено во 2-е изд. книги. Смотрите комментарий.]

Предпосылки Предположим , что θ N является оценка получается путем минимизации 1$\widehat \theta_N$

Обозначим гессиан $q$ как $H$ ,

Асимптотическая ковариация & thetas п включает Е [ Н ( д , & thetas ; 0 ) ] , где θ 0 является истинным значением параметра. Один из способов оценить это - использовать эмпирический гессенский$\widehat \theta_n$$E[H(q,\theta_0)]$$\theta_0$

Это определенность Н , которая находится под вопросом.$\widehat H$

Я думаю ты прав. Давайте разберем ваш аргумент в его сущности:

1. минимизирует функциюQопределяется какQ(thetas)=$\widehat \theta_N$$Q$$Q(\theta) = {1 \over N}\sum_{i=1}^N q(w_i,\theta).$

2. Пусть - гессиан Q , откуда H ( θ ) = ∂ 2$H$$Q$ , по определениюи этов свою очередь, по линейности дифференциации, равна$H(\theta) = \frac{\partial^2 Q}{\partial \theta_i \partial \theta_j}$.$\frac{1}{N}\sum_{i=1}^N H(w_i, \theta_n)$

3. Предполагая , что & thetas ; N лежит во внутренней части области Q , то Н ( θ N ) должен быть положительным полуопределенным.$\widehat \theta_N$$Q$$H(\widehat \theta_N)$

Это просто утверждение о функции : как она определяется всего лишь отвлечение, за исключением тех случаев, когда предполагается второго порядка дифференцируемость ц относительно второго аргумента ( & thetas ) обеспечивает второй порядок дифференцируемости Q .$Q$$q$$\theta$$Q$

Найти М-оценки может быть сложно. Рассмотрим эти данные, предоставленные @mpiktas:

{1.168042, 0.3998378}, {1.807516, 0.5939584}, {1.384942, 3.6700205}, {1.327734, -3.3390724}, {1.602101, 4.1317608}, {1.604394, -1.9045958}, {1.124633, -3.0865249}, {1.294601, -1.8331763},{1.577610, 1.0865977}, { 1.630979, 0.7869717}

Процедура R для нахождения M-оценки с дает решение ( c 1 , c 2 ) = ( - 114,91316 , - 32,54386 ) . Значение целевой функции (среднее значение q ) в этой точке равно 62,3542. Вот сюжет подгонки:$q((x,y),\theta)=(y-c_1x^{c_2})^4$$(c_1, c_2)$$(-114.91316, -32.54386)$$q$

Вот график (log) целевой функции в окрестности этого соответствия:

Здесь что-то подозрительно: параметры подгонки чрезвычайно далеки от параметров, используемых для моделирования данных (около ), и мы, кажется, не находимся на минимуме: мы находимся в чрезвычайно мелкой долине с наклоном в сторону больших значений обоих параметров:$(0.3, 0.2)$

Отрицательный определитель гессиана в этой точке подтверждает, что это не локальный минимум! Тем не менее, когда вы посмотрите на метки оси z, вы увидите, что эта функция имеет точность от пяти до пяти цифр во всем регионе, поскольку она равна константе 4,1329 (логарифм 62,354). Это, вероятно, привело к тому, что минимизатор функции R (с допусками по умолчанию) пришел к выводу, что он близок к минимуму.

На самом деле, решение далеко от этой точки. Чтобы быть уверенным в его нахождении, я использовал в Mathematica дорогой в вычислительном отношении, но очень эффективный метод « Основной оси » , используя точность в 50 цифр (основание 10), чтобы избежать возможных численных проблем. Он находит минимум вблизи ( c 1 , c 2 ) = ( 0,02506 , 7,55973 ), где целевая функция имеет значение 58.292655: примерно на 6% меньше, чем «минимум», найденный R. Этот минимум имеет место в чрезвычайно плоском участке , но я могу заставить его выглядеть (едва ли) как истинный минимум с эллиптическими контурами, преувеличивая c 2$(c_1, c_2) = (0.02506, 7.55973)$$c_2$ Направление на участке:

Контуры варьируются от 58.29266 в середине до 58.29284 в углах (!). Вот трехмерный вид (опять же цель журнала):

Здесь гессиан положительно определен: его собственные значения 55062.02 и 0.430978. Таким образом, эта точка является локальным минимумом (и, вероятно, глобальным минимумом). Вот подгонка, которой он соответствует:

Я думаю, что это лучше, чем другой. Значения параметров, безусловно, более реалистичны, и ясно, что мы не сможем добиться большего успеха с этим семейством кривых.

Из этого примера мы можем извлечь полезные уроки:

1. Численная оптимизация может быть трудной, особенно с нелинейной подгонкой и неквадратичными функциями потерь. Следовательно:
2. Дважды проверьте результаты как можно больше способов, в том числе:
3. График целевой функции, когда вы можете.
4. Если числовые результаты нарушают математические теоремы, будьте крайне подозрительны.
5. Когда статистические результаты являются неожиданными - например, неожиданные значения параметров, возвращаемые кодом R - быть особенно подозрительными.

$\hat{\theta}_N$

$\hat{\theta}_N$$\Theta$$\hat{H}$

$N^{-1}\sum_{i=1}^Nq(w_i,\theta)$$\theta_0$

$N^{-1}\sum_{i=1}^Nq(w_i,\theta)$$\Theta$

Далее в своей книге Вулдридж приводит примеры оценок Гессиана, которые гарантированно являются численно положительно определенными. На практике неположительная определенность гессиана должна указывать на то, что решение либо находится на граничной точке, либо алгоритм не может найти решение. Что обычно является дополнительным указанием на то, что установленная модель может быть неподходящей для данных.

Вот числовой пример. Я генерирую нелинейную задачу наименьших квадратов:

$X$$[1,2]$$\varepsilon$$\sigma^2$set.seed(3)$x_i$$y_i$

Я выбрал квадрат целевой функции обычной нелинейной целевой функции наименьших квадратов:

Вот код в R для оптимизации функции, ее градиент и гессиан.

##First set-up the epxressions for optimising function, its gradient and hessian.
##I use symbolic derivation of R to guard against human error    
mt <- expression((y-c1*x^c2)^4)

gradmt <- c(D(mt,"c1"),D(mt,"c2"))

hessmt <- lapply(gradmt,function(l)c(D(l,"c1"),D(l,"c2")))

##Evaluate the expressions on data to get the empirical values. 
##Note there was a bug in previous version of the answer res should not be squared.
optf <- function(p) {
    res <- eval(mt,list(y=y,x=x,c1=p[1],c2=p[2]))
    mean(res)
}

gf <- function(p) {
    evl <- list(y=y,x=x,c1=p[1],c2=p[2]) 
    res <- sapply(gradmt,function(l)eval(l,evl))
    apply(res,2,mean)
}

hesf <- function(p) {
    evl <- list(y=y,x=x,c1=p[1],c2=p[2]) 
    res1 <- lapply(hessmt,function(l)sapply(l,function(ll)eval(ll,evl)))
    res <- sapply(res1,function(l)apply(l,2,mean))
    res
}

Сначала проверьте, что градиент и гессиан работают так, как рекламируется.

set.seed(3)
x <- runif(10,1,2)
y <- 0.3*x^0.2

> optf(c(0.3,0.2))
[1] 0
> gf(c(0.3,0.2))
[1] 0 0
> hesf(c(0.3,0.2))
     [,1] [,2]
[1,]    0    0
[2,]    0    0
> eigen(hesf(c(0.3,0.2)))$values
[1] 0 0

$x$$y$

> df <- read.csv("badhessian.csv")
> df
          x          y
1  1.168042  0.3998378
2  1.807516  0.5939584
3  1.384942  3.6700205
4  1.327734 -3.3390724
5  1.602101  4.1317608
6  1.604394 -1.9045958
7  1.124633 -3.0865249
8  1.294601 -1.8331763
9  1.577610  1.0865977
10 1.630979  0.7869717
> x <- df$x
> y <- df$y
> opt <- optim(c(1,1),optf,gr=gf,method="BFGS")  
> opt$par
[1] -114.91316  -32.54386
> gf(opt$par)
[1] -0.0005795979 -0.0002399711
> hesf(opt$par)
              [,1]         [,2]
[1,]  0.0002514806 -0.003670634
[2,] -0.0036706345  0.050998404
> eigen(hesf(opt$par))$values
[1]  5.126253e-02 -1.264959e-05

Градиент равен нулю, но гессиан не является положительным.

Примечание: это моя третья попытка дать ответ. Надеюсь, мне наконец удалось дать точные математические утверждения, которые ускользали от меня в предыдущих версиях.

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

Обновление: позвольте мне уточнить. Во-первых, давайте предположим, что эмпирический гессиан существует везде.

$\hat{\theta}_n$$\sum_i q(w_i, \cdot)$$(1/N) \sum_i H(w_i, \hat{\theta}_n)$$\hat{\theta}_n$$\sum_i q(w_i, \cdot)$$\hat{\theta}_n$

$\arg\min_\theta \sum_i q(w_i, \theta)$

Практически говоря, даже положительно определенный гессиан, который является почти единственным или плохо обусловленным, предполагает, что оценщик является плохим, и вам нужно больше беспокоиться, чем оценивать его дисперсию.

В этой теме было много ударов вокруг куста о том, должен ли гессиан быть положительным (полу) определенным при локальном минимуме. Поэтому я сделаю четкое заявление по этому поводу.

Если предположить, что целевая функция и все функции ограничения дважды непрерывно дифференцируемы, то при любом локальном минимуме гессиан лагранжиана, проецируемый в нулевое пространство якобиана активных ограничений, должен быть положительно полуопределенным. То есть если$Z$ является базисом для нулевого пространства якобиана активных ограничений, то $Z^T*(\text{Hessian of Lagrangian})*Z$должен быть положительным полуопределенным. Это должно быть положительно определенным для строгого локального минимума.

Таким образом, гессиан целевой функции в ограниченной задаче, имеющей активные ограничения, не обязательно должен быть положительным полуопределенным, если существуют активные ограничения.

Примечания:

1) Активные ограничения состоят из всех ограничений равенства, плюс ограничения неравенства, которые удовлетворяются равенством.

2) См. Определение лагранжиана на странице https://www.encyclopediaofmath.org/index.php/Karush-Kuhn-Tucker_conditions .

3) Если все ограничения линейны, то гессиан лагранжиана = гессиан целевой функции, поскольку 2-ые производные линейных функций равны нулю. Но вам все равно нужно выполнить джаз проекции, если какое-либо из этих ограничений активно. Обратите внимание, что ограничения нижней или верхней границ являются частными случаями ограничений линейного неравенства. Если единственными активными ограничениями являются связанные ограничения, то проекция гессиана в пустое пространство якобиана активных ограничений сводится к удалению строк и столбцов гессиана, соответствующих этим компонентам на их границах.

4) Поскольку множители Лагранжа неактивных ограничений равны нулю, если нет активных ограничений, гессиан лагранжиана = гессиан целевой функции, а матрица тождества является основой для нулевого пространства якобиана активных ограничений, которое приводит к упрощению критерия, являющегося знакомым условием, чтобы гессиан целевой функции был положительно полуопределенным на локальном минимуме (положительно определенным, если строгий локальный минимум).

Вышеуказанные положительные ответы верны, но они не учитывают основополагающего предположения об идентификации - если ваша модель не идентифицирована (или если она только идентифицирована как набор), вы, действительно, как правильно указал Вулдридж, окажетесь с эмпирическим гессианом без PSD. Просто запустите какую-нибудь не игрушечную психометрическую / эконометрическую модель и убедитесь сами.