quantile regression --python实现
前言
不管是quantile regression,还是在其他方法上把loss function改为quantile loss function,都是根据设定的quantile, 返回一个值。并不是直接返回一个范围。即一个分位数quantile,对应一个模型。
分位数回归可调用的库
1. scikit-learn
官网: https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.QuantileRegressor.html#s%20%20klearn.linear_model.QuantileRegressor
实例code:
from sklearn.linear_model import QuantileRegressor
import numpy as np
n_samples, n_features = 10, 1
rng = np.random.RandomState(0)
y = rng.randn(n_samples)
x = rng.randn(n_samples, n_features)
reg = QuantileRegressor(quantile=0.8).fit(x, y)
2. statsmodels
官网: https://www.statsmodels.org/stable/examples/notebooks/generated/quantile_regression.html
实例code:
import statsmodels.api as sm
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
data = sm.datasets.engel.load_pandas().data
mod = smf.quantreg("foodexp ~ income", data)
res = mod.fit(q=0.5)
print(res.summary())
quantile loss function - python实现
这个要注意,不同的基础函数所对应的loss function的输入输出格式是不同的。
1. 在neural network 中添加分位数损失函数
代码
def tilted_loss(quantile, y, f):
"""
the final overall loss is the mean of the sum.
y: y_true
f: y_pred
"""
e = (y - f)
return K.mean(K.maximum(quantile * e, (quantile - 1) * e), axis=-1)
def mcycle_model(length):
model = Sequential()
model.add(LSTM(units=128, input_shape=(1, length)))
model.add(Dense(units=10, input_dim=1, activation='relu'))
model.add(Dense(1))
return model
if __name__ == '__main__':
x_train, x_test, y_train, y_test = train_test_split(X, Y, random_state=2)
x_train = np.array(x_train)[:, np.newaxis, :]
x_test = np.array(x_test)[:, np.newaxis, :]
qs = [0.1, 0.5, 0.9]
plt.plot(y_test.values, 'h', label='Origin')
len_x = len(cols_x)
model = mcycle_model(len_x)
for q in qs:
model.compile(loss=lambda y, f: tilted_loss(q, y, f), optimizer='adadelta')
model.fit(x_train, y_train, epochs=20, batch_size=32, verbose=0)
y_pred = model.predict(x_test)
plt.plot(y_pred, '.--', label=q, alpha=0.3)
plt.legend()
plt.show()
2. 在xgboost上添加分位数损失函数
我们使用构造函数构造 XGBRegressor 的时候,里边的 objective 参数才是真正的损失函数(loss function)。xgb 使用 sklearn api 的时候需要用到的损失函数,其返回值是一阶导和二阶导。 这里首先要参考xgboost官网(https://xgboost.readthedocs.io/en/latest/tutorials/custom_metric_obj.html)来看他的loss function的格式,输入输出是什么。 如果直接把上面的loss function带入的话,会报错。
File "E:/ml_regression.py", line 77, in xgb_reg
reg.fit(x_train, y_train, eval_set=[(x_test, y_test)])
File "C:\Users\AppData\Local\Programs\Python\Python37\lib\site-packages\xgboost\core.py", line 575, in inner_f
return f(**kwargs)
File "C:\Users\AppData\Local\Programs\Python\Python37\lib\site-packages\xgboost\sklearn.py", line 972, in fit
callbacks=callbacks,
File "C:\Users\AppData\Local\Programs\Python\Python37\lib\site-packages\xgboost\core.py", line 575, in inner_f
return f(**kwargs)
File "C:\Users\AppData\Local\Programs\Python\Python37\lib\site-packages\xgboost\training.py", line 181, in train
bst.update(dtrain, i, obj)
File "C:\Users\AppData\Local\Programs\Python\Python37\lib\site-packages\xgboost\core.py", line 1783, in update
grad, hess = fobj(pred, dtrain)
File "C:\Users\AppData\Local\Programs\Python\Python37\lib\site-packages\tensorflow\python\framework\ops.py", line 511, in __iter__
raise TypeError("Cannot iterate over a scalar tensor.")
TypeError: Cannot iterate over a scalar tensor.
代码
def quantile_loss(y_true, y_pred, alpha, delta, threshold, var):
"""
loss function:
x is the error between y_true and y_pred
if 0<x<alpha: grad1 = -x, hess1 = 1; grad2 = grad1, hess2 = hess1
elif x >= max(alpha, threshold): grad1 = 0, hess1 = 0; grad2 = (2 * np.random.randint(2, size=len(y_true)) - 1.0) * var, hess2 = 1
:param y_true:
:param y_pred:
:param alpha: quantile
:param delta:
:param threshold:
:param var: i think this is variance
:return: grad: i think that this is negtive gradient
"""
x = y_true - y_pred
grad = (x < (alpha - 1.0) * delta) * (1.0 - alpha) - (
(x >= (alpha - 1.0) * delta) & (x < alpha * delta)) * x / delta - alpha * (x > alpha * delta)
hess = ((x >= (alpha - 1.0) * delta) & (x < alpha * delta)) / delta
grad = (np.abs(x) < threshold) * grad - (np.abs(x) >= threshold) * (
2 * np.random.randint(2, size=len(y_true)) - 1.0) * var
hess = (np.abs(x) < threshold) * hess + (np.abs(x) >= threshold)
return grad, hess
q=0.1
reg = xgb.XGBRegressor(
objective=lambda y,f:quantile_loss(y,f,q))
补充
上面我们提到,xgboost的loss function要求return grad, hess。 其中grad即一阶梯度,hess为二阶梯度,即hessian。在数学中, 海森矩阵(Hessian matrix或Hessian)是一个自变量为向量的实值函数的二阶偏导数组成的方块矩阵。
reference
1.杜克大学的材料,包括公式、实例。 https://www2.stat.duke.edu/~fl35/teaching/440-19F/decks/cs01_5_deck.html#/example-engel-curves 2. University of Minho大学课件 http://www.econ.uiuc.edu/~roger/courses/NIPE/lectures/L1.pdf
|