分位数曲线相互交叉的问题

机器算法验证 线性的 分位数回归 曲线
2022-03-22 19:39:48

在尝试对某些数据执行一些分位数回归时,我遇到了一个经典问题:0.25、0.50、0.75、0.95 分位数曲线相互交叉(我正在处理线性曲线)。

我目前正在使用sqreg命令与 Stata 合作。

我已阅读 Howard Bondell 的论文“非交叉分位数回归曲线估计”(2010 年),我知道他使用 R 编写了解决方案。

我还阅读了 Andriyana 和 Chernozhukov 的论文(关于分位数回归中的非单调性问题)。

为了清楚我们在说什么,下面是一个不交叉的分位数曲线示例(我的问题是当它们交叉时该怎么办):

图片

我想知道你们中是否有人可以推荐工作实践中最常用的软件解决方案,以处理交叉分位数曲线(=避免它们交叉),我知道有多种方法,所以问题是:什么方法你更喜欢什么软件来解决这个问题

注意:我没有发布我的数据样本,因为它无关紧要,我只是对问题的一般理论以及如何专门使用软件(理想情况下是 STATA 或 Mathematica,或 R 或​​ Python 或 Perl)处理它感兴趣。

1个回答

一般情况下,非交叉线性分位数回归问题没有解

实际上,您的分位数线要么是平行的(并且微不足道),要么它们确实在某处交叉。

通过设置限制,您可以确保训练样本中的线不会交叉,但绝不会保证在训练模型后看到的下一次观察中不会发生交叉。如果分位数线在训练样本中相交,则很可能意味着您的模型指定不正确:均值或标准差非线性变化,或者您在拟合模型时应用了错误的成本函数。

如果您仍想估计这样的模型,我建议将其拟合为神经网络。

你的模型应该接受输入X,乘以矩阵(!)β得到一个预测矩阵f=Xβ有大小n×k, 在哪里n是观察次数和k是估计的分位数。我假设分位数百分比q是递增的顺序。您应该最小化该功能

L=i=1n(j=1kmax(qj(yifij),(qj1)(yifij))+j=1k1αmax(0,δ(fi,j+1fij)))

内部和中的第一项只是普通的分位数回归损失。第二项是如果两个连续的分位数预测相差小于时应用的惩罚。δ

通过梯度下降最小化此函数将为您提供不交叉的分位数线(如果足够大)。但我仍然警告您,如果“自然”分位数线相交,那么您的模型的功能形式可能会出现问题。也许您更喜欢随机森林的分位数估计(如在 R 中),它们总是一致的。αquantregForest

有一个 Python 示例。至于我,受限制和不受限制的版本看起来都很丑。

# import everything
import keras
from keras import backend as K
from keras import Sequential
from keras.layers import Dense
import numpy as np
import matplotlib.pyplot as plt

# create the dataset
n = 200
np.random.seed(1)
X = np.sort(np.random.normal(size=n))[:, np.newaxis] + 4
y = 5 + X.ravel()*(1 + np.random.normal(size=n)*0.2)
quantiles = np.array([0.1, 0.25, 0.5, 0.75, 0.9])

# define loss function
def quantile_ensemble_loss(q, y, f, alpha=100, margin=0.1):
    error = (y - f)
    quantile_loss = K.mean(K.maximum(q*error, (q-1)*error))
    diff = f[:, 1:] - f[:, :-1]
    penalty = K.mean(K.maximum(0.0, margin - diff)) * alpha
    return quantile_loss + penalty

# fit two models
for i, alpha, name in [(1, 0, 'w/o penalty'), (2, 10, 'with_penalty')]:
    model = Sequential()
    model.add(Dense(len(quantiles), input_dim=X.shape[1]))
    model.compile(loss=lambda y,f: quantile_ensemble_loss(quantiles,y,f,alpha), optimizer=keras.optimizers.RMSprop(lr=0.003))
    model.fit(X, y, epochs=3000, verbose=0, batch_size=100);
    plt.subplot(1,2,i)
    plt.scatter(X.ravel(), y, s=1)
    plt.plot(X.ravel(), model.predict(X))
    plt.title(name)

plt.show()

在此处输入图像描述