I tried to build a quantile regression myself and compare the result with the sklearn: https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.QuantileRegressor.html#sklearn.linear_model.QuantileRegressor
This is the code of the quantile regression
class QuantileRegression:
def __init__(self, theta, verbose=False):
self.theta = theta # theta-quantile
self.verbose = verbose # to print the loss
def fit(self, X, y, learning_rate=0.001, num_iter=1000):
self.dim = X.shape[1]
self.num_iter = num_iter
self.learning_rate = learning_rate
# initialize the params for gradient descent
self.params = np.random.uniform(0, 1, self.dim)
self.gradient_descent(X, y)
def sum_abs_loss(self, y, y_hat):
# sum of absolute error
dev = y - y_hat
return np.sum(np.maximum(self.theta * dev, (self.theta - 1) * dev))
def gradient_descent(self, X, y):
self.count = 0
y_hat = X @ self.params
best_loss = self.sum_abs_loss(y, y_hat)
if self.verbose:
print(f'{self.count}/{self.num_iter}:', best_loss)
self.best_params = np.copy(self.params)
for i in range(self.num_iter):
# update params
grads = self.gradient(X, y, self.params)
self.params -= grads * self.learning_rate
# print loss
y_hat = X @ self.params
self.count += 1
loss = self.sum_abs_loss(y, y_hat)
if self.verbose:
print(f'{self.count}/{self.num_iter}:', loss)
if loss < best_loss:
self.best_params = np.copy(self.params)
best_loss = loss
def gradient(self, X, y, params):
y_hat = X @ self.params
dev = y - y_hat
grads = np.zeros(self.dim)
# this gradient for sum(|y - xb|) only
for i in range(len(grads)):
grads[i] = np.sum(np.where(dev > 0,
self.theta * np.sign(- X[:, i] * self.params[i]),
(self.theta-1) * np.sign(- X[:, i] * self.params[i])))
return grads
this is the test case given that the true beta is [10, 1, 2, 3, 4, 5]
X = np.c_[np.ones(100), np.random.normal(0, 5, size=(100, 5))]
beta = np.array([10, 1, 2, 3, 4, 5])
y = X @ beta + np.random.normal(0, 1, 100)
for q in [0.01, 0.05, 0.2, 0.5, 0.7, 0.95, 0.99]:
qr = QuantileRegression(q)
qr.fit(X, y)
print(qr.params)
qr_sklearn1 = QuantileRegressor(quantile=q, alpha=0, fit_intercept=False)
qr_sklearn1.fit(X,y)
print(qr_sklearn1.coef_)
when I conducted the experiment and print the coefficients of my QR and sklearn's QR:
when theta = 0.01
my qr: [-0.6498561 1.14255861 1.96963312 2.42596852 3.87759425 4.07229887]
sklearn qr [8.27991656 1.02353173 2.04483203 3.00242634 4.03772444 4.96552982]
when theta = 0.05
my qr: [3.1335077 1.0661197 1.72773243 2.75970899 3.91690744 4.31994194]
sklearn qr [8.39770131 1.04236782 2.02841619 3.00370016 4.01664672 4.93823634]
when theta = 0.2
my qr: [9.15718701 1.02003304 1.97818333 2.94035991 4.00303393 4.99738457]
sklearn qr [9.17828948 1.03784577 1.96222859 2.97330991 4.00220429 4.99989697]
when theta = 0.5
my qr: [10.10550684 1.03086459 2.01128614 2.93384029 3.96102259 5.00320016]
sklearn qr [10.02735652 1.03415381 1.9980452 2.96871738 3.99386424 4.99261563]
when theta = 0.7
my qr: [10.6616965 0.99940346 2.01492713 2.92801633 3.99833129 4.93782001]
sklearn qr [10.6246958 1.00520561 2.01437818 2.92953241 3.99722636 4.94740772]
when theta = 0.95
my qr: [11.98306491 0.8544847 1.92519793 2.87262039 4.08435037 4.8680982 ]
sklearn qr [11.63871321 0.94709746 2.08566667 2.94908355 3.98721936 5.01007832]
when theta = 0.99
my qr: [12.32849148 0.99429387 1.77393081 2.94237375 4.06992124 4.89181963]
sklearn qr [11.97795949 0.91029108 1.99793332 2.94873895 3.98622897 5.02365461]
I realised when theta (quantile) < 0.2, the performance is terribly bad.
May I know if I have made any mistake in my code building the quantile regression?
