精华内容
下载资源
问答
  • Python 分位数回归
    更多相关内容
  • 分位数回归

    2018-11-15 13:58:57
    实现分布滞后分位数回归,包括方差标准差,t检验各系数是否显著
  • 基于神经网络分位数回归及核密度估计的概率密度预测方法,闻才喜,何耀耀,本文引入神经网络分位数回归和核密度估计方法,把神经网络强大的非线性自适应能力及分位数回归能更加细致刻画解释变量的优点结合
  • o(* ̄︶ ̄*)o我们要探测自变量 与因变量 的关系,最简单的方法是线性回归,即假设:我们通过最小二乘方法 (OLS: ordinary least squares) 的无偏估计 , 。为了解决, 的可靠性问题,我们同时对残差 做了假设,即:...

    天朗气清,惠风和畅。赋闲在家,正宜读书。前人文章,不得其解。代码开源,无人注释。你们不来,我行我上。废话少说,直入主题。o(* ̄︶ ̄*)o

    我们要探测自变量

    equation?tex=y 与因变量

    equation?tex=x 的关系,最简单的方法是线性回归,即假设:

    equation?tex=y_%7Bi%7D+%3D+%5Cbeta_%7B0%7D+%2B+%5Cbeta_%7B1%7D+x_%7Bi%7D+%2B+e_%7Bi%7D%2C+i%3D1%2C...%2Cn

    我们通过最小二乘方法 (OLS: ordinary least squares)

    equation?tex=%5Cbeta_%7B0%7D%E4%B8%8E%5Cbeta_%7B1%7D 的无偏估计

    equation?tex=%5Ctilde%7B%5Cbeta_%7B0%7D%7D

    equation?tex=%5Ctilde%7B%5Cbeta_%7B1%7D%7D 。为了解决

    equation?tex=%5Ctilde%7B%5Cbeta_%7B0%7D%7D

    equation?tex=%5Ctilde%7B%5Cbeta_%7B1%7D%7D的可靠性问题,我们同时对残差

    equation?tex=e_%7Bi%7D 做了假设,即:

    equation?tex=e_%7Bi%7D为均值为0,方差恒定的独立随机变量。

    equation?tex=%5Ctilde%7By_%7Bi%7D%7D+%3D%5Ctilde%7B%5Cbeta_%7B0%7D%7D+%2B+%5Ctilde%7B%5Cbeta_%7B1%7D%7Dx_%7Bi%7D 即为给定自变量

    equation?tex=x_%7Bi%7D 下,因变量

    equation?tex=y_%7Bi%7D 的条件均值。

    假如残差

    equation?tex=e_%7Bi%7D 不满足我们的假设,或者更重要地,我们不仅仅想要知道

    equation?tex=y_%7Bi%7D 的在给定

    equation?tex=x_%7Bi%7D下的条件均值,而且想知道条件中位数(更一般地,条件分位数),那么OLS下的线性回归就不能满足我们的需求。分位数回归(Quantile Regression)

    1. 一个例子:收入与食品消费

    这个例子出自statasmodels:Quantile Regression.

    1.1 预处理

    %matplotlib inline

    import numpy as np

    import pandas as pd

    import statsmodels.api as sm

    import statsmodels.formula.api as smf

    import matplotlib.pyplot as plt

    data = sm.datasets.engel.load_pandas().data

    data.head()

    income foodexp

    0420.157651255.839425

    1541.411707310.958667

    2901.157457485.680014

    3639.080229402.997356

    4750.875606495.560775

    1.2 中位数回归 (分位数回归的特例,q=0.5)

    mod = smf.quantreg('foodexp ~ income', data)

    res = mod.fit(q=.5)

    print(res.summary())

    QuantReg Regression Results

    ==============================================================================

    Dep. Variable: foodexp Pseudo R-squared: 0.6206

    Model: QuantReg Bandwidth: 64.51

    Method: Least Squares Sparsity: 209.3

    Date: Mon, 21 Oct 2019 No. Observations: 235

    Time: 17:46:59 Df Residuals: 233

    Df Model: 1

    ==============================================================================

    coef std err t P>|t| [0.025 0.975]

    ------------------------------------------------------------------------------

    Intercept 81.4823 14.634 5.568 0.000 52.649 110.315

    income 0.5602 0.013 42.516 0.000 0.534 0.586

    ==============================================================================

    The condition number is large, 2.38e+03. This might indicate that there are

    strong multicollinearity or other numerical problems.

    由结果可以知道

    equation?tex=foodexp_%7Bmedian%7D%3D81.4623+%2B+0.5602income+ ,如何得到回归系数的估计?结果中的std err, t, Pseudo R-squared等是什么?我会在稍后解释。

    1.3 数据可视化

    我们先拟合10个分位数回归,分位数q分别在0.05到0.95之间。

    quantiles = np.arange(.05, .96, .1)

    def fit_model(q):

    res = mod.fit(q=q)

    return [q, res.params['Intercept'], res.params['income']] + \

    res.conf_int().loc['income'].tolist()

    models = [fit_model(x) for x in quantiles]

    models = pd.DataFrame(models, columns=['q', 'a', 'b', 'lb', 'ub'])

    ols = smf.ols('foodexp ~ income', data).fit()

    ols_ci = ols.conf_int().loc['income'].tolist()

    ols = dict(a = ols.params['Intercept'],

    b = ols.params['income'],

    lb = ols_ci[0],

    ub = ols_ci[1])

    print(models)

    print(ols)

    q a b lb ub

    0 0.05 124.880096 0.343361 0.268632 0.418090

    1 0.15 111.693660 0.423708 0.382780 0.464636

    2 0.25 95.483539 0.474103 0.439900 0.508306

    3 0.35 105.841294 0.488901 0.457759 0.520043

    4 0.45 81.083647 0.552428 0.525021 0.579835

    5 0.55 89.661370 0.565601 0.540955 0.590247

    6 0.65 74.033435 0.604576 0.582169 0.626982

    7 0.75 62.396584 0.644014 0.622411 0.665617

    8 0.85 52.272216 0.677603 0.657383 0.697823

    9 0.95 64.103964 0.709069 0.687831 0.730306

    {'a': 147.47538852370562, 'b': 0.48517842367692354, 'lb': 0.4568738130184233,

    这里拟合了10个回归,其中q是对应的分位数,a是斜率,b是回归系数。lb和ub分别是b的95%置信区间的下界与上界。

    现在来画出这10条回归线:

    x = np.arange(data.income.min(), data.income.max(), 50)

    get_y = lambda a, b: a + b * x

    fig, ax = plt.subplots(figsize=(8, 6))

    for i in range(models.shape[0]):

    y = get_y(models.a[i], models.b[i])

    ax.plot(x, y, linestyle='dotted', color='grey')

    y = get_y(ols['a'], ols['b'])

    ax.plot(x, y, color='red', label='OLS')

    ax.scatter(data.income, data.foodexp, alpha=.2)

    ax.set_xlim((240, 3000))

    ax.set_ylim((240, 2000))

    legend = ax.legend()

    ax.set_xlabel('Income', fontsize=16)

    ax.set_ylabel('Food expenditure', fontsize=16);

    上图中虚线是分位数回归线,红线是线性最小二乘(OLS)的回归线。通过观察,我们可以发现3个现象:随着收入提高,食品消费也在提高。

    随着收入提高,家庭间食品消费的差别拉大。穷人别无选择,富人能选择生活方式,有喜欢吃贵的,也有喜欢吃便宜的。然而我们无法通过OLS发现这个现象,因为它只给了我们一个均值。

    对于穷人来说,OLS预测值过高。这是因为少数的富人拉高了整体的均值,可见OLS对异常点敏感,不是一个稳健的模型。

    2.分位数回归的原理

    这部分是数理统计的内容,只关心如何实现的朋友可以略过。我们要解决以下这几个问题:什么是分位数?

    如何求分位数?

    什么是分位数回归?

    分位数回归的回归系数如何求得?

    回归系数的检验如何进行?

    如何评估回归拟合优度?

    2.1 分位数的定义

    equation?tex=Y 是随机变量,

    equation?tex=Y 的累积密度函数是

    equation?tex=F_%7By%7D%28y%29+%3D+P%28Y%5Cleq+y%29 .

    equation?tex=Y

    equation?tex=q 分位数为:

    equation?tex=Y_%7Bq%7D+%3DF_%7BY%7D%5E%7B-1%7D%28q%29%3Dinf%5Cleft%5C%7B+y%3AF_%7BY%7D+%28y%29%5Cgeq+q%5Cright%5C%7D ,

    equation?tex=q%5Cin%280%2C1%29

    假设有100个人,95%的人身高少于1.9m, 1.9m就是身高的95%分位数。

    2.2 分位数的求法

    equation?tex=%5Ctilde%7BY_%7Bq%7D%7D+%3D+argmin_%7By_%7Bq%7D%7D%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%7B%5Crho_%7Bq%7D%28y_%7Bi%7D-y_%7Bq%7D%7D%29+%3Dargmin_%7By_%7Bq%7D%7D%5B%28q-1%29%5Csum_%7By_%7Bi%7D%3Cy%7Bq%7D%7D%5E%7B%7D%7B%28y_%7Bi%7D-y_%7Bq%7D%29%2Bq%5Csum_%7By_%7Bi%7D%5Cgeq+y_%7Bq%7D%7D%5E%7B%7D%7B%28y_%7Bi%7D-y_%7Bq%29%7D%7D%5D%7D+

    通过选择不同的

    equation?tex=y_%7Bq%7D 值,使

    equation?tex=%28q-1%29%5Csum_%7By_%7Bi%7D%3Cy%7Bq%7D%7D%5E%7B%7D%7B%28y_%7Bi%7D-y_%7Bq%7D%29%2Bq%5Csum_%7By_%7Bi%7D%5Cgeq+y_%7Bq%7D%7D%5E%7B%7D%7B%28y_%7Bi%7D-y_%7Bq%29%7D%7D%7D 最小,对应的

    equation?tex=y_%7Bq%7D 值即为

    equation?tex=Y

    equation?tex=q 分位数的估计

    equation?tex=%5Ctilde%7BY_%7Bq%7D%7D .

    2.3 分位数回归

    对于OLS, 我们有:

    equation?tex=%5Ctilde%7BY%7D%7CX%3DX%5Ctilde%7B%5Cbeta%7D

    equation?tex=%5Ctilde%7B%5Cbeta%7D

    equation?tex=%7C%7CY-X%5Cbeta%7C%7C%5E%7B2%7D 最小化所对应的

    equation?tex=%5Cbeta ,类比地,对于

    equation?tex=q 分位数回归,我们有:

    equation?tex=%5Ctilde%7BY_%7Bq%7D%7D%7CX%3DX%5Ctilde%7B%5Cbeta_%7Bq%7D%7D

    equation?tex=%5Ctilde%7B%5Cbeta_%7Bq%7D%7D 为最小化:

    equation?tex=%7C%7C%5Crho_%7Bq%7D%28Y-X%5Cbeta%29%7C%7C 即最小化

    equation?tex=%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%7B%5Crho_%7Bq%7D%28y_%7Bi%7D-x_%7Bi%7D%27%5Cbeta%7D%29%3D%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%5B%7B%28q-1%29%5Csum_%7By_%7Bi%7D%3Cx_%7Bi%7D%E2%80%99%5Cbeta%7D%5E%7B%7D%28%7By_%7Bi%7D-x_%7Bi%7D%27%5Cbeta%29%2Bq%5Csum_%7By_%7Bi%7D%5Cgeq+x_%7Bi%7D%5Cbeta%7D%5E%7B%7D%7B%28y_%7Bi%7D-x_%7Bi%7D%27%5Cbeta%29%5D%7D%7D%7D

    所对应的

    equation?tex=%5Cbeta

    2.4 系数估计

    由于

    equation?tex=%7C%7C%5Crho_%7Bq%7D%28Y-X%5Cbeta%29%7C%7C 不能直接对

    equation?tex=%5Cbeta 求导,我们只能用迭代的方法来逼近

    equation?tex=%7C%7C%5Crho_%7Bq%7D%28Y-X%5Cbeta%29%7C%7C 最小时对应的

    equation?tex=%5Cbeta 值。statsmodels采用了Iteratively reweighted least squares (IRLS)的方法。

    假设我们要求

    equation?tex=%5Cbeta 最小化形如下的

    equation?tex=L%5E%7Bp%7D 范数:

    equation?tex=argmin_%7B%5Cbeta%7D%7C%7CY-X%5Cbeta%7C%7C_%7Bp%7D%3Dargmin_%7B%5Cbeta%7D%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%7B%7Cy_%7Bi%7D-x_%7Bi%7D%27%5Cbeta%7D%7C%5E%7Bp%7D

    则第t+1步迭代的

    equation?tex=%5Cbeta 值为:

    equation?tex=%5Cbeta%5E%7Bt%2B1%7D%3Dargmin_%7B%5Cbeta%7D%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%7Bw_%7Bi%7D%5E%7B%28t%29%7D%7D%7Cy_%7Bi%7D-x_%7Bi%7D%27%5Cbeta%7C%5E%7B2%7D%3D%28X%5E%7BT%7DW%5E%7B%28t%29%7DX%29%5E%7B-1%7DX%5E%7BT%7DW%5E%7B%28t%29%7Dy

    equation?tex=W%5E%7Bt%7D 是对角矩阵且初始值为

    equation?tex=w_%7Bi%7D%5E%7B%280%29%7D%3D1

    第t次迭代

    equation?tex=w_%7Bi%7D%5E%7B%28t%29%7D%3D%7Cy_%7Bi%7D-x_%7Bi%7D%27%5Cbeta%7C%5E%7Bp-2%7D

    以中位数回归为例子(q=0.5),我们求:

    equation?tex=argmin_%7B%5Cbeta%7D%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%5B%7B%28q-1%29%5Csum_%7By_%7Bi%7D%3Cx_%7Bi%7D%27%5Cbeta%7D%5E%7B%7D%28%7By_%7Bi%7D-x_%7Bi%7D%27%5Cbeta%29%2Bq%5Csum_%7By_%7Bi%7D%5Cgeq+x_%7Bi%7D%27%5Cbeta%7D%5E%7B%7D%7B%28y_%7Bi%7D-x_%7Bi%7D%27%5Cbeta%29%5D%7D%7D%7D

    equation?tex=argmin_%7B%5Cbeta%7D%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%5B%7B-0.5%5Csum_%7By_%7Bi%7D%3Cx_%7Bi%7D%27%5Cbeta%7D%5E%7B%7D%28%7By_%7Bi%7D-x_%7Bi%7D%27%5Cbeta%29%2B0.5%5Csum_%7By_%7Bi%7D%5Cgeq+x_%7Bi%7D%27%5Cbeta%7D%5E%7B%7D%7B%28y_%7Bi%7D-x_%7Bi%7D%27%5Cbeta%29%5D%7D%7D%7D%3Dargmin_%7B%5Cbeta%7D%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%7B%7Cy_%7Bi%7D-x_%7Bi%7D%27%5Cbeta%7C%7D

    即最小化形如上的

    equation?tex=L%5E%7B1%7D 范数,

    equation?tex=w_%7Bi%7D%5E%7B%28t%29%7D%3D%7Cy_%7Bi%7D-x_%7Bi%7D%27%5Cbeta%7C%5E%7B-1%7D

    为避免分母为0,我们取

    equation?tex=w_%7Bi%7D%5E%7B%28t%29%7D%3D1%2Fmax%28%7B%7B%5Cdelta%2C%7Cy_%7Bi%7D-x_%7Bi%7D%27%5Cbeta%7C%7D%7D%29 ,

    equation?tex=%5Cdelta 是一个很小的数,例如0.0001.

    2.5 回归系数的检验

    我们通过2.4,多次迭代得出

    equation?tex=%5Cbeta 的估计值,为了得到假设检验的t统计量,我们只需得到

    equation?tex=%5Cbeta 的方差的估计。

    equation?tex=q分位数回归

    equation?tex=%5Cbeta_%7Bq%7D 的协方差矩阵的渐近估计为:

    equation?tex=Est.Asy.Var%5B%5Cbeta_%7Bq%7D%5D%3D%28X%27X%29%5E%7B-1%7DX%27DX%28X%27X%29%5E%7B-1%7D

    其中

    equation?tex=D 是对角矩阵,

    equation?tex=y_%7Bi%7D%3Ex_%7Bi%7D%27%5Cbeta ,

    equation?tex=d_%7Bi%7D%3D%5B%5Cfrac%7Bq%7D%7Bf%280%29%7D%5D%5E%7B2%7D , 当

    equation?tex=y_%7Bi%7D%5Cleq+x_%7Bi%7D%27%5Cbeta

    equation?tex=d_%7Bi%7D%3D%5B%5Cfrac%7B1-q%7D%7Bf%280%29%7D%5D%5E%7B2%7D

    equation?tex=f%280%29 的估计为

    equation?tex=%5Ctilde%7Bf%280%29%7D%3D%5Cfrac%7B1%7D%7Bn%7D%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%7B%5Cfrac%7B1%7D%7Bh%7DK%5B%5Cfrac%7Be_%7Bi%7D%7D%7Bh%7D%7D%5D

    其中

    equation?tex=e_%7Bi%7D%3Dy_%7Bi%7D-x_%7Bi%7D%27%5Cbeta

    equation?tex=K%28%29 为核函数(Kernel),可取Epa,Gaussian等.

    equation?tex=h 为根据Stata 12所选的窗宽(bandwidth)

    回归结果中的std err即由

    equation?tex=Est.Asy.Var%5B%5Cbeta_%7Bq%7D%5D 获得,t统计量等于

    equation?tex=%5Ctilde%7B%5Cbeta%7D%2F%5Ctilde%7BVar%28%5Cbeta%29%7D

    2.6 拟合优度

    对于OLS,我们用

    equation?tex=R%5E%7B2%7D%3D1-%5Cfrac%7B%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%7B%28y_%7Bi%7D-%5Ctilde%7By_%7Bi%7D%7D%29%5E%7B2%7D%7D%7D%7B%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%7B%28y_%7Bi%7D-%5Cbar%7By_%7Bi%7D%7D%7D%29%5E%7B2%7D%7D 来衡量拟合优度。对于

    equation?tex=q 分位数回归,我们类比得到:

    equation?tex=R_%7Bq%7D%5E%7B2%7D%3D1-%5Cfrac%7B%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%7B%5Crho_%7Bq%7D%28y_%7Bi%7D-x_%7Bi%7D%27%5Cbeta%7D%29%7D%7B%5Csum_%7Bi%3D1%7D%5E%7Bn%7D%7B%5Crho_%7Bq%7D%28y_%7Bi%7D-y_%7Bq%7D%7D%29%7D ,其中

    equation?tex=y_%7Bq%7D 为所有

    equation?tex=y 观察值的

    equation?tex=q 分位数。

    equation?tex=R_%7Bq%7D%5E%7B2%7D 即为回归结果中的Pseudo R-squared。

    3.Python源码分析

    实现分位数回归的完整源码在

    3.1 QuantReg类

    #QuantReg是包中RegressionModel的一个子类

    class QuantReg(RegressionModel):

    #计算回归系数及其协方差矩阵。q是分位数,vcov是协方差矩阵,默认robust即2.5的方法。核函数kernel默认

    #epa,窗宽bandwidth默认hsheather.IRLS最大迭代次数默认1000,差值默认小于1e-6时停止迭代

    def fit(self, q=.5, vcov='robust', kernel='epa', bandwidth='hsheather',

    max_iter=1000, p_tol=1e-6, **kwargs):

    """Solve by Iterative Weighted Least SquaresParameters----------q : floatQuantile must be between 0 and 1vcov : str, method used to calculate the variance-covariance matrixof the parameters. Default is ``robust``:- robust : heteroskedasticity robust standard errors (as suggestedin Greene 6th edition)- iid : iid errors (as in Stata 12)kernel : str, kernel to use in the kernel density estimation for theasymptotic covariance matrix:- epa: Epanechnikov- cos: Cosine- gau: Gaussian- par: Parzenebandwidth : str, Bandwidth selection method in kernel densityestimation for asymptotic covariance estimate (fullreferences in QuantReg docstring):- hsheather: Hall-Sheather (1988)- bofinger: Bofinger (1975)- chamberlain: Chamberlain (1994)"""

    if q < 0 or q > 1:

    raise Exception('p must be between 0 and 1')

    kern_names = ['biw', 'cos', 'epa', 'gau', 'par']

    if kernel not in kern_names:

    raise Exception("kernel must be one of " + ', '.join(kern_names))

    else:

    kernel = kernels[kernel]

    if bandwidth == 'hsheather':

    bandwidth = hall_sheather

    elif bandwidth == 'bofinger':

    bandwidth = bofinger

    elif bandwidth == 'chamberlain':

    bandwidth = chamberlain

    else:

    raise Exception("bandwidth must be in 'hsheather', 'bofinger', 'chamberlain'")

    #endog样本因变量,exog样本自变量

    endog = self.endog

    exog = self.exog

    nobs = self.nobs

    exog_rank = np.linalg.matrix_rank(self.exog)

    self.rank = exog_rank

    self.df_model = float(self.rank - self.k_constant)

    self.df_resid = self.nobs - self.rank

    #IRLS初始化

    n_iter = 0

    xstar = exog

    beta = np.ones(exog_rank)

    # TODO: better start, initial beta is used only for convergence check

    # Note the following does not work yet,

    # the iteration loop always starts with OLS as initial beta

    # if start_params is not None:

    # if len(start_params) != rank:

    # raise ValueError('start_params has wrong length')

    # beta = start_params

    # else:

    # # start with OLS

    # beta = np.dot(np.linalg.pinv(exog), endog)

    diff = 10

    cycle = False

    history = dict(params = [], mse=[])

    #IRLS迭代

    while n_iter < max_iter and diff > p_tol and not cycle:

    n_iter += 1

    beta0 = beta

    xtx = np.dot(xstar.T, exog)

    xty = np.dot(xstar.T, endog)

    beta = np.dot(pinv(xtx), xty)

    resid = endog - np.dot(exog, beta)

    mask = np.abs(resid) < .000001

    resid[mask] = ((resid[mask] >= 0) * 2 - 1) * .000001

    resid = np.where(resid < 0, q * resid, (1-q) * resid)

    resid = np.abs(resid)

    #1/resid[:, np.newaxis]为更新权重W

    xstar = exog / resid[:, np.newaxis]

    diff = np.max(np.abs(beta - beta0))

    history['params'].append(beta)

    history['mse'].append(np.mean(resid*resid))

    #检查是否收敛,若收敛则提前停止迭代

    if (n_iter >= 300) and (n_iter % 100 == 0):

    # check for convergence circle, should not happen

    for ii in range(2, 10):

    if np.all(beta == history['params'][-ii]):

    cycle = True

    warnings.warn("Convergence cycle detected", ConvergenceWarning)

    break

    if n_iter == max_iter:

    warnings.warn("Maximum number of iterations (" + str(max_iter) +

    ") reached.", IterationLimitWarning)

    #计算协方差矩阵

    e = endog - np.dot(exog, beta)

    # Greene (2008, p.407) writes that Stata 6 uses this bandwidth:

    # h = 0.9 * np.std(e) / (nobs**0.2)

    # Instead, we calculate bandwidth as in Stata 12

    iqre = stats.scoreatpercentile(e, 75) - stats.scoreatpercentile(e, 25)

    h = bandwidth(nobs, q)

    h = min(np.std(endog),

    iqre / 1.34) * (norm.ppf(q + h) - norm.ppf(q - h))

    fhat0 = 1. / (nobs * h) * np.sum(kernel(e / h))

    if vcov == 'robust':

    d = np.where(e > 0, (q/fhat0)**2, ((1-q)/fhat0)**2)

    xtxi = pinv(np.dot(exog.T, exog))

    xtdx = np.dot(exog.T * d[np.newaxis, :], exog)

    vcov = chain_dot(xtxi, xtdx, xtxi)

    elif vcov == 'iid':

    vcov = (1. / fhat0)**2 * q * (1 - q) * pinv(np.dot(exog.T, exog))

    else:

    raise Exception("vcov must be 'robust' or 'iid'")

    #用系数估计值和协方差矩阵创建一个QuantResults对象,并输出结果

    lfit = QuantRegResults(self, beta, normalized_cov_params=vcov)

    lfit.q = q

    lfit.iterations = n_iter

    lfit.sparsity = 1. / fhat0

    lfit.bandwidth = h

    lfit.history = history

    return RegressionResultsWrapper(lfit)

    #核函数表达式

    def _parzen(u):

    z = np.where(np.abs(u) <= .5, 4./3 - 8. * u**2 + 8. * np.abs(u)**3,

    8. * (1 - np.abs(u))**3 / 3.)

    z[np.abs(u) > 1] = 0

    return z

    kernels = {}

    kernels['biw'] = lambda u: 15. / 16 * (1 - u**2)**2 * np.where(np.abs(u) <= 1, 1, 0)

    kernels['cos'] = lambda u: np.where(np.abs(u) <= .5, 1 + np.cos(2 * np.pi * u), 0)

    kernels['epa'] = lambda u: 3. / 4 * (1-u**2) * np.where(np.abs(u) <= 1, 1, 0)

    kernels['gau'] = lambda u: norm.pdf(u)

    kernels['par'] = _parzen

    #kernels['bet'] = lambda u: np.where(np.abs(u) <= 1, .75 * (1 - u) * (1 + u), 0)

    #kernels['log'] = lambda u: logistic.pdf(u) * (1 - logistic.pdf(u))

    #kernels['tri'] = lambda u: np.where(np.abs(u) <= 1, 1 - np.abs(u), 0)

    #kernels['trw'] = lambda u: 35. / 32 * (1 - u**2)**3 * np.where(np.abs(u) <= 1, 1, 0)

    #kernels['uni'] = lambda u: 1. / 2 * np.where(np.abs(u) <= 1, 1, 0)

    #窗宽计算

    def hall_sheather(n, q, alpha=.05):

    z = norm.ppf(q)

    num = 1.5 * norm.pdf(z)**2.

    den = 2. * z**2. + 1.

    h = n**(-1. / 3) * norm.ppf(1. - alpha / 2.)**(2./3) * (num / den)**(1./3)

    return h

    def bofinger(n, q):

    num = 9. / 2 * norm.pdf(2 * norm.ppf(q))**4

    den = (2 * norm.ppf(q)**2 + 1)**2

    h = n**(-1. / 5) * (num / den)**(1. / 5)

    return h

    def chamberlain(n, q, alpha=.05):

    return norm.ppf(1 - alpha / 2) * np.sqrt(q*(1 - q) / n)

    3.2 QuantRegResults类

    这里我只给出计算拟合优度的代码。

    class QuantRegResults(RegressionResults):

    '''Results instance for the QuantReg model'''

    @cache_readonly

    def prsquared(self):

    q = self.q

    endog = self.model.endog

    #e为残差

    e = self.resid

    e = np.where(e < 0, (1 - q) * e, q * e)

    e = np.abs(e)

    ered = endog - stats.scoreatpercentile(endog, q * 100)

    ered = np.where(ered < 0, (1 - q) * ered, q * ered)

    ered = np.abs(ered)

    return 1 - np.sum(e) / np.sum(ered)

    4.总结

    上文我先给出了一个分位数回归的应用例子,进而叙述了分位数回归的原理,最后再分析了Python实现的源码。

    分位数回归对比起OLS回归,虽然较为复杂,但它有三个主要优势:能反映因变量分位数与自变量的关系,而不仅仅反映因变量均值与自变量的关系。

    分位数回归对残差不作任何假设。

    分位数回归受异常点的影响较小。

    (欢迎转载,但请注明出处)

    参考

    展开全文
  • 分位数回归 参考文献 Python statsmodels 介绍 - 树懒学堂 (shulanxt.com) Quantile Regression - IBM Documentation https://www.cnblogs.com/TMesh/p/11737368.html 传统的线性回归模型 其的求解方式是一个最小...

    分位数回归

    参考文献

    Python statsmodels 介绍 - 树懒学堂 (shulanxt.com)

    Quantile Regression - IBM Documentation

    https://www.cnblogs.com/TMesh/p/11737368.html

    传统的线性回归模型

    其的求解方式是一个最小二乘法,保证观测值与你的被估值的差的平方和应该保持最小,
    M S E   =   1 n ∑ i = 1 n ( y i − f ^ ( x i ) ) 2   =   E ( y − f ^ ( x ) ) 2 MSE\ =\ \frac{1}{n}\sum_{i=1}^n{\left( y_i-\widehat{f}\left( x_i \right) \right) ^2\ =\ E\left( y-\widehat{f}\left( x \right) \right)}^2 MSE = n1i=1n(yif (xi))2 = E(yf (x))2

    • 因变量的条件均值分布受自变量x的影响过程,因此我们拟合出来的曲线是在给定x的情况下,y的条件均值
    • 随机误差项来均值为0、同方差,因此估计的斜率在现有的基础上是最好的

    分位数回归

    • 首先提出中位数回归(最小绝对偏差估计)
    • 改进出分位数回归
    • 描述自变量X对因变量Y的变化范围,以及其不受分布形状的影响。即其不止可以描述条件均值的影响,还可以描述中位数的影响

    因此我们能够得到如下的一个损失函数
    Q Y ^ ( τ ) = a r g min ⁡ ξ r ∈ R ( ∑ i : Y i ≥ ξ t a u τ ∣ Y i − ξ r ∣ + ∑ i : Y i < ξ t a u ( 1 − τ ) ∣ Y i − ξ r ∣ ) \widehat{Q_Y}\left( \tau \right) =arg\min _{\xi _{r\in R}}\left( \sum_{i:Y_i\ge \xi _{t^{au}}}{\tau \left| Y_i-\xi _r \right|}+\sum_{i:Y_i<\xi _{t^{au}}}{\left( 1-\tau \right) \left| Y_i-\xi _r \right|} \right) QY (τ)=argξrRmini:YiξtauτYiξr+i:Yi<ξtau(1τ)Yiξr
    参数 τ \tau τ的估计算法有:

    • 单纯形算法
    • 内点算法
    • 平滑算法

    总结来说,在我心目中,分位数回归是对传统回归的一种改进,它不在局限于原来最小二乘法,使得数据可以更多影响其他的点或者类似于中位数的影响。

    接下来我们将采用python语言进行实现,采用的数据集是我们之前的文章中cpu—time_tamp的数据

    class QuantileRegression:
        def __init__(self,data):
            # self.data = pd.DataFrame(data=np.hstack([time_stamp,cpu_util_percent]),columns=["time_stamp","cpu_util_percent"])
            self.data = data
            # self.num = len(time_stamp)
            pass
        def __QuantileReq_1__(self):
            # 主义这里,前面是Y轴,后面是X轴
            mod = smf.quantreg('cpu_util_percent~time_stamp',self.data)
            print(mod)
            res = mod.fit()
            print(res)
            fig = plt.subplots(figsize=(8, 6))
            # x = np.arange(self.data.time_stamp.min(),self.data.time_stamp.max(),1000)
    
            print(res.summary())
    

    在这里插入图片描述

    数据解释:

    • Dep. Variable :因变量

    • Model:方法模块

    • Method:方法(最小二乘法)默认使用迭代加权最小二乘法(IRLS)

    • Date:日期

    • Time:时间

    • Pseudo R-squared: 拟合优度

      采用的公式为:
      R q 2 = 1 − ∑ i = 1 n ρ q ( y i − x i ′ β ) ∑ i = 1 n ρ q ( y i − y q ) R_{q}^{2}=1-\frac{\sum_{i=1}^{n} \rho_{q}\left(y_{i}-x_{i}^{\prime} \beta\right)}{\sum_{i=1}^{n} \rho_{q}\left(y_{i}-y_{q}\right)} Rq2=1i=1nρq(yiyq)i=1nρq(yixiβ)

    • Bandwidth:窗宽h

      公式来源于:
      当 y i > x i ′ β , d i = [ q f ( 0 ) ] 2 , 当 y i ≤ x i ′ β , d i = [ 1 − q f ( 0 ) ] 2 f ( 0 ) 的估计为  f ( 0 ) ~ = 1 n ∑ i = 1 n 1 h K [ e i h ] 当 y_{i}>x_{i}^{\prime} \beta , d_{i}=\left[\frac{q}{f(0)}\right]^{2} , 当 y_{i} \leq x_{i}^{\prime} \beta , d_{i}=\left[\frac{1-q}{f(0)}\right]^{2} f(0)_{\text {的估计为 }} \tilde{f(0)}=\frac{1}{n} \sum_{i=1}^{n} \frac{1}{h} K\left[\frac{e_{i}}{h}\right] yi>xiβ,di=[f(0)q]2,yixiβdi=[f(0)1q]2f(0)的估计为 f(0)~=n1i=1nh1K[hei]

      其 中 , f ( 0 ) 的估计为  f ( 0 ) ~ = 1 n ∑ i = 1 n 1 h K [ e i h ] 其 中 e i = y i − x i ′ β ​ , K [ ] 表 示 为 核 函 数 其中, f(0)_{\text {的估计为 }} \tilde{f(0)}=\frac{1}{n} \sum_{i=1}^{n} \frac{1}{h} K\left[\frac{e_{i}}{h}\right] 其中e_i=y_i-x_{i}^{'}\beta ​ ,K[]表示为核函数 ,f(0)的估计为 f(0)~=n1i=1nh1K[hei]ei=yixiβ,K[]

    • Sparsity

    • No. Observations:

    • Df Residuals :Df残差

    • Df Model

    • coef:系数

    • std err:协方差(标准差)

      采用以下公式得到:
      E s t . A s y . Var ⁡ [ β q ] = ( X ′ X ) − 1 X ′ D X ( X ′ X ) − 1 ​ Est. Asy.\operatorname{Var}\left[\beta_{q}\right]=\left(X^{\prime} X\right)^{-1} X^{\prime} D X\left(X^{\prime} X\right)^{-1} ​ Est.Asy.Var[βq]=(XX)1XDX(XX)1 其中D为对角矩阵,

    • t:统计量,表示为 β ~ V a r ~ ( β ) \dfrac{\widetilde{\beta }}{Va\tilde{r}}\left( \beta \right) Var~β (β)

    • P>|t|

    • [0.025 0.975]

    • Intercept:截距

    • cpu_util_percent : 斜率

    但在多次实验的过程中,发现一直报过时未收敛的警告,所以我查看了源代码,最终我们怀疑python的分位数回归可能不太适用于曲线回归,可能只能分段式线性回归比较合适,以下是源代码的部分

    #!/usr/bin/env python
    
    '''
    Quantile regression model
    
    Model parameters are estimated using iterated reweighted least squares. The
    asymptotic covariance matrix estimated using kernel density estimation.
    
    Author: Vincent Arel-Bundock
    License: BSD-3
    Created: 2013-03-19
    
    The original IRLS function was written for Matlab by Shapour Mohammadi,
    University of Tehran, 2008 (shmohammadi@gmail.com), with some lines based on
    code written by James P. Lesage in Applied Econometrics Using MATLAB(1999).PP.
    73-4.  Translated to python with permission from original author by Christian
    Prinoth (christian at prinoth dot name).
    '''
    
    import numpy as np
    import warnings
    import scipy.stats as stats
    from numpy.linalg import pinv
    from scipy.stats import norm
    from statsmodels.tools.decorators import cache_readonly
    from statsmodels.regression.linear_model import (RegressionModel,
                                                     RegressionResults,
                                                     RegressionResultsWrapper)
    from statsmodels.tools.sm_exceptions import (ConvergenceWarning,
                                                 IterationLimitWarning)
    
    [docs]class QuantReg(RegressionModel):
        # 计算回归系数及其协方差矩阵。q是分位数,vcov是协方差矩阵,默认robust即2.5的方法,核函数kernel默认
        # epa,窗宽bandwidth默认hsheather.IRLS最大迭代次数默认1000,差值默认小于1e-6时停止迭代
        '''Quantile Regression
    
    	使用迭代加权最小二乘法估计分位数回归模型。
        Estimate a quantile regression model using iterative reweighted least
        squares.
    
        Parameters
        ----------
        endog : array or dataframe   数据/数据帧
            endogenous/response variable  内源性/响应变量
        exog : array or dataframe
            exogenous/explanatory variable(s) 外生/解释变量(s)
    
        Notes
        -----
        The Least Absolute Deviation (LAD) estimator is a special case where
        quantile is set to 0.5 (q argument of the fit method).
        最小绝对偏差(LAD)估计量是一种特殊情况
    	Quantile被设置为0.5 (fit方法的q参数)。
    
        The asymptotic covariance matrix is estimated following the procedure in
        Greene (2008, p.407-408), using either the logistic or gaussian kernels
        (kernel argument of the fit method).
        在此基础上,对渐近协方差矩阵进行了估计
    	格林(2008,p.407-408),使用logistic或高斯核
    	(拟合方法的核心参数)。
    
        References
        ----------
        General:
    
        * Birkes, D. and Y. Dodge(1993). Alternative Methods of Regression, John Wiley and Sons.
        * Green,W. H. (2008). Econometric Analysis. Sixth Edition. International Student Edition.
        * Koenker, R. (2005). Quantile Regression. New York: Cambridge University Press.
        * LeSage, J. P.(1999). Applied Econometrics Using MATLAB,
        * Birkes, D.和Y. Dodge(1993)。回归的可选方法,约翰·威利和儿子。
    	* Green,W。h(2008)。计量经济学分析。第六版。国际学生版。
    	* Koenker, R.(2005)。分位数回归。纽约:剑桥大学出版社。
    	* LeSage J. P.(1999)。应用计量经济学
    
        Kernels (used by the fit method):
    
        * Green (2008) Table 14.2
    
        Bandwidth selection (used by the fit method):
        
    
        * Bofinger, E. (1975). Estimation of a density function using order statistics. Australian Journal of Statistics 17: 1-17.
        * Chamberlain, G. (1994). Quantile regression, censoring, and the structure of wages. In Advances in Econometrics, Vol. 1: Sixth World Congress, ed. C. A. Sims, 171-209. Cambridge: Cambridge University Press.
        * Hall, P., and S. Sheather. (1988). On the distribution of the Studentized quantile. Journal of the Royal Statistical Society, Series B 50: 381-391.
    
        Keywords: Least Absolute Deviation(LAD) Regression, Quantile Regression,
        Regression, Robust Estimation.
        * Bofinger E.(1975)。使用顺序统计量估计密度函数。澳大利亚统计杂志17:1-17。
    
    	*张伯伦,G.(1994)。分位数回归、审查和工资结构。《计量经济学进展》,第1卷:第六届世界大会,c.a.西姆斯编,171-209。剑桥:剑桥大学出版社。
    
    	* Hall, P.和S. Sheather。(1988)。研究分位数的分布。皇家统计学会学报,B辑50:381-391。
    
    
    
    	关键词:最小绝对偏差回归分位数回归
    
    	回归,稳健估计。
        '''
    
        # 初始化
        def __init__(self, endog, exog, **kwargs):
            self._check_kwargs(kwargs)
            super(QuantReg, self).__init__(endog, exog, **kwargs)
    
    [docs]    def whiten(self, data):
            """
            QuantReg model whitener does nothing: returns data.
            QuantReg模型增白器什么也不做:返回数据。
            """
            return data
    
    
    [docs]    def fit(self, q=.5, vcov='robust', kernel='epa', bandwidth='hsheather',
                max_iter=1000, p_tol=1e-6, **kwargs):
            """
            Solve by Iterative Weighted Least Squares
            用迭代加权最小二乘法求解
    
            Parameters
            ----------
            q : float
                Quantile must be strictly between 0 and 1
            vcov : str, method used to calculate the variance-covariance matrix
                of the parameters. Default is ``robust``:
    
                - robust : heteroskedasticity robust standard errors (as suggested
                  in Greene 6th edition)
                - iid : iid errors (as in Stata 12)
                q:浮动型小数
    
    			分位数必须严格在0和1之间
    
    		vcoc:str,用于计算方差-协方差矩阵的参数方法。默认是“robust”:
    
    		-robust:异方差鲁棒性标准误差(如在格林第六版中的建议)
    
    		- iid: iid错误(如Stata 12)
    
            kernel : str, kernel to use in the kernel density estimation for the
                asymptotic covariance matrix:
            Kernel: str,用于核密度估计的渐近协方差矩阵的核:
    
                - epa: Epanechnikov
                - cos: Cosine  余旋
                - gau: Gaussian 高斯
                - par: Parzene
    
            bandwidth : str, Bandwidth selection method in kernel density
                estimation for asymptotic covariance estimate (full
                references in QuantReg docstring):
            bandwidth: str,渐近协方差估计核密度估计中的带宽选择方法(完整参考QuantReg文档字符串):
    
                - hsheather: Hall-Sheather (1988)
                - bofinger: Bofinger (1975)
                - chamberlain: Chamberlain (1994)
            """
    
            if q <= 0 or q >= 1:
                raise Exception('q must be strictly between 0 and 1')
    
            kern_names = ['biw', 'cos', 'epa', 'gau', 'par']
            if kernel not in kern_names:
                raise Exception("kernel must be one of " + ', '.join(kern_names))
            else:
                kernel = kernels[kernel]
    
            if bandwidth == 'hsheather':
                bandwidth = hall_sheather
            elif bandwidth == 'bofinger':
                bandwidth = bofinger
            elif bandwidth == 'chamberlain':
                bandwidth = chamberlain
            else:
                raise Exception("bandwidth must be in 'hsheather', 'bofinger', 'chamberlain'")
    
            #endog样本因变量,exog样本自变量
            endog = self.endog
            exog = self.exog
            nobs = self.nobs
            exog_rank = np.linalg.matrix_rank(self.exog)
            self.rank = exog_rank
            self.df_model = float(self.rank - self.k_constant)
            self.df_resid = self.nobs - self.rank
            #IRLS初始化
            n_iter = 0
            xstar = exog
    
            beta = np.ones(exog.shape[1])
            # TODO: better start, initial beta is used only for convergence check
            # 待办事项:更好的开始,初始测试版仅用于收敛检查
    
            # Note the following does not work yet,
            # the iteration loop always starts with OLS as initial beta
            # if start_params is not None:
            #    if len(start_params) != rank:
            #       raise ValueError('start_params has wrong length')
            #       beta = start_params
            #    else:
            #       # start with OLS
            #       beta = np.dot(np.linalg.pinv(exog), endog)
            """
            #注意以下内容还不能使用,
    
    		迭代循环总是以OLS作为初始测试开始
    
    		#如果start_params不是None:
    
    		# if len(start_params) != rank:
    
    		#引发ValueError('start_params has wrong length')
    
    		# beta = start_params
    
    		其他:
    
    		# #从OLS开始
    
    		# beta = np.dot(np. linalgr .pinv(exog), endog)
            """
    
            diff = 10
            cycle = False
    
            history = dict(params = [], mse=[])
            #IRLS迭代
            while n_iter < max_iter and diff > p_tol and not cycle:
                n_iter += 1
                beta0 = beta
                xtx = np.dot(xstar.T, exog)
                xty = np.dot(xstar.T, endog)
                beta = np.dot(pinv(xtx), xty)
                resid = endog - np.dot(exog, beta)
    
                mask = np.abs(resid) < .000001
                resid[mask] = ((resid[mask] >= 0) * 2 - 1) * .000001
                resid = np.where(resid < 0, q * resid, (1-q) * resid)
                resid = np.abs(resid)
                #1/resid[:, np.newaxis]为更新权重W
                xstar = exog / resid[:, np.newaxis]
                diff = np.max(np.abs(beta - beta0))
                history['params'].append(beta)
                history['mse'].append(np.mean(resid*resid))
    
                #检查是否收敛,若收敛则提前停止迭代
                if (n_iter >= 300) and (n_iter % 100 == 0):
                    # check for convergence circle, should not happen
                    for ii in range(2, 10):
                        if np.all(beta == history['params'][-ii]):
                            cycle = True
                            warnings.warn("Convergence cycle detected", ConvergenceWarning)
                            break
    		# 超出迭代次数,发出警告并结束,迭代次数默认为1000
            if n_iter == max_iter:
                warnings.warn("Maximum number of iterations (" + str(max_iter) +
                              ") reached.", IterationLimitWarning)
    
            #计算协方差矩阵
            e = endog - np.dot(exog, beta)
            # Greene (2008, p.407) writes that Stata 6 uses this bandwidth:
            # h = 0.9 * np.std(e) / (nobs**0.2)
            # Instead, we calculate bandwidth as in Stata 12
            iqre = stats.scoreatpercentile(e, 75) - stats.scoreatpercentile(e, 25)
            h = bandwidth(nobs, q)
            h = min(np.std(endog),
                    iqre / 1.34) * (norm.ppf(q + h) - norm.ppf(q - h))
    
            fhat0 = 1. / (nobs * h) * np.sum(kernel(e / h))
    
            if vcov == 'robust':
                d = np.where(e > 0, (q/fhat0)**2, ((1-q)/fhat0)**2)
                xtxi = pinv(np.dot(exog.T, exog))
                xtdx = np.dot(exog.T * d[np.newaxis, :], exog)
                vcov = xtxi @ xtdx @ xtxi
            elif vcov == 'iid':
                vcov = (1. / fhat0)**2 * q * (1 - q) * pinv(np.dot(exog.T, exog))
            else:
                raise Exception("vcov must be 'robust' or 'iid'")
    
                #用系数估计值和协方差矩阵创建一个QuantResults对象,并输出结果
            lfit = QuantRegResults(self, beta, normalized_cov_params=vcov)
    
            lfit.q = q
            lfit.iterations = n_iter
            lfit.sparsity = 1. / fhat0
            lfit.bandwidth = h
            lfit.history = history
    
            return RegressionResultsWrapper(lfit)
    
    
    
        #核函数表达式
    def _parzen(u):
        z = np.where(np.abs(u) <= .5, 4./3 - 8. * u**2 + 8. * np.abs(u)**3,
                     8. * (1 - np.abs(u))**3 / 3.)
        z[np.abs(u) > 1] = 0
        return z
    
    
    kernels = {}
    kernels['biw'] = lambda u: 15. / 16 * (1 - u**2)**2 * np.where(np.abs(u) <= 1, 1, 0)
    kernels['cos'] = lambda u: np.where(np.abs(u) <= .5, 1 + np.cos(2 * np.pi * u), 0)
    kernels['epa'] = lambda u: 3. / 4 * (1-u**2) * np.where(np.abs(u) <= 1, 1, 0)
    kernels['gau'] = lambda u: norm.pdf(u)
    kernels['par'] = _parzen
    #kernels['bet'] = lambda u: np.where(np.abs(u) <= 1, .75 * (1 - u) * (1 + u), 0)
    #kernels['log'] = lambda u: logistic.pdf(u) * (1 - logistic.pdf(u))
    #kernels['tri'] = lambda u: np.where(np.abs(u) <= 1, 1 - np.abs(u), 0)
    #kernels['trw'] = lambda u: 35. / 32 * (1 - u**2)**3 * np.where(np.abs(u) <= 1, 1, 0)
    #kernels['uni'] = lambda u: 1. / 2 * np.where(np.abs(u) <= 1, 1, 0)
    
    
    #窗宽计算
    def hall_sheather(n, q, alpha=.05):
        z = norm.ppf(q)
        num = 1.5 * norm.pdf(z)**2.
        den = 2. * z**2. + 1.
        h = n**(-1. / 3) * norm.ppf(1. - alpha / 2.)**(2./3) * (num / den)**(1./3)
        return h
    
    
    def bofinger(n, q):
        num = 9. / 2 * norm.pdf(2 * norm.ppf(q))**4
        den = (2 * norm.ppf(q)**2 + 1)**2
        h = n**(-1. / 5) * (num / den)**(1. / 5)
        return h
    
    
    def chamberlain(n, q, alpha=.05):
        return norm.ppf(1 - alpha / 2) * np.sqrt(q*(1 - q) / n)
    
    
    [docs]class QuantRegResults(RegressionResults):
        '''Results instance for the QuantReg model'''
    
        @cache_readonly
        def prsquared(self):
            q = self.q
            endog = self.model.endog
            e = self.resid
            e = np.where(e < 0, (1 - q) * e, q * e)
            e = np.abs(e)
            ered = endog - stats.scoreatpercentile(endog, q * 100)
            ered = np.where(ered < 0, (1 - q) * ered, q * ered)
            ered = np.abs(ered)
            return 1 - np.sum(e) / np.sum(ered)
    
        #@cache_readonly
    [docs]    def scale(self):
            return 1.
    
    
        @cache_readonly
        def bic(self):
            return np.nan
    
        @cache_readonly
        def aic(self):
            return np.nan
    
        @cache_readonly
        def llf(self):
            return np.nan
    
        @cache_readonly
        def rsquared(self):
            return np.nan
    
        @cache_readonly
        def rsquared_adj(self):
            return np.nan
    
        @cache_readonly
        def mse(self):
            return np.nan
    
        @cache_readonly
        def mse_model(self):
            return np.nan
    
        @cache_readonly
        def mse_total(self):
            return np.nan
    
        @cache_readonly
        def centered_tss(self):
            return np.nan
    
        @cache_readonly
        def uncentered_tss(self):
            return np.nan
    
        @cache_readonly
        def HC0_se(self):
            raise NotImplementedError
    
        @cache_readonly
        def HC1_se(self):
            raise NotImplementedError
    
        @cache_readonly
        def HC2_se(self):
            raise NotImplementedError
    
        @cache_readonly
        def HC3_se(self):
            raise NotImplementedError
    
    [docs]    def summary(self, yname=None, xname=None, title=None, alpha=.05):
            """Summarize the Regression Results
    
            Parameters
            ----------
            yname : str, optional
                Default is `y`
            xname : list[str], optional
                Names for the exogenous variables. Default is `var_##` for ## in
                the number of regressors. Must match the number of parameters
                in the model
            title : str, optional
                Title for the top table. If not None, then this replaces the
                default title
            alpha : float
                significance level for the confidence intervals
    
            Returns
            -------
            smry : Summary instance
                this holds the summary tables and text, which can be printed or
                converted to various output formats.
    
            See Also
            --------
            statsmodels.iolib.summary.Summary : class to hold summary results
            """
            """
            总结回归结果
    
    		参数
    		----------
    		Yname: str,可选
    		默认是“y”
    		list[str],可选
    		外生变量的名称。默认是' var_## '的## in
    		回归量的数量。必须匹配的参数个数
    		在模型中
    		标题:str,可选
    		顶级的头衔。如果不是None,则替换
    		默认的标题
    		α:浮动
    		置信区间的显著性水平
    
    		返回
    		-------
    		smry:概要实例
    		这包含了汇总表和文本,可以打印或
    		转换为各种输出格式。
    
    		另请参阅
    		--------
    		summary:保存汇总结果的类
            """
            
            eigvals = self.eigenvals
            condno = self.condition_number
    
            top_left = [('Dep. Variable:', None),
                        ('Model:', None),
                        ('Method:', ['Least Squares']),
                        ('Date:', None),
                        ('Time:', None)
                        ]
    
            top_right = [('Pseudo R-squared:', ["%#8.4g" % self.prsquared]),
                         ('Bandwidth:', ["%#8.4g" % self.bandwidth]),
                         ('Sparsity:', ["%#8.4g" % self.sparsity]),
                         ('No. Observations:', None),
                         ('Df Residuals:', None),
                         ('Df Model:', None)
                         ]
    
            if title is None:
                title = self.model.__class__.__name__ + ' ' + "Regression Results"
    
            # create summary table instance
            from statsmodels.iolib.summary import Summary
            smry = Summary()
            smry.add_table_2cols(self, gleft=top_left, gright=top_right,
                                 yname=yname, xname=xname, title=title)
            smry.add_table_params(self, yname=yname, xname=xname, alpha=alpha,
                                  use_t=self.use_t)
    
            # add warnings/notes, added to text format only
            etext = []
            if eigvals[-1] < 1e-10:
                wstr = "The smallest eigenvalue is %6.3g. This might indicate "
                wstr += "that there are\n"
                wstr += "strong multicollinearity problems or that the design "
                wstr += "matrix is singular."
                wstr = wstr % eigvals[-1]
                etext.append(wstr)
            elif condno > 1000:  # TODO: what is recommended
                wstr = "The condition number is large, %6.3g. This might "
                wstr += "indicate that there are\n"
                wstr += "strong multicollinearity or other numerical "
                wstr += "problems."
                wstr = wstr % condno
                etext.append(wstr)
    
            if etext:
                smry.add_extra_txt(etext)
    
            return smry
    
    展开全文
  • 多种方法实现分位数回归
  • 偶尔在机器学习的论文中了解到了分位数回归,发现这个方法应用也满广的。 一般的回归方法是最小二乘法,即最小化误差的平方和: min⁡∑(yi−y^i)2\min\quad \sum(y_i-\hat{y}_i)^2min∑(yi​−y^​i​)2 其中,yiy_...

    偶尔在机器学习的论文中了解到了分位数回归,发现这个方法应用也满广的。

    1. 分位数回归的数学原理

    一般的回归方法是最小二乘法,即最小化误差的平方和:

    min ⁡ ∑ ( y i − y ^ i ) 2 \min\quad \sum(y_i-\hat{y}_i)^2 min(yiy^i)2

    其中, y i y_i yi 是真实值,而 y ^ i \hat{y}_i y^i 是预测值。而分位数的目标是最小化加权的误差绝对值和:

    min ⁡ y ^ i ∑ y i ≥ y ^ i τ ∣ y i − y ^ i ∣ + ∑ y i < y ^ i ( 1 − τ ) ∣ y i − y ^ i ∣ \min_{\hat{y}_i}\quad \sum_{y_i\geq \hat{y}_i} \tau|y_i-\hat{y}_i|+\sum_{y_i<\hat{y}_i}(1-\tau)|y_i-\hat{y}_i| y^iminyiy^iτyiy^i+yi<y^i(1τ)yiy^i

    其中, τ \tau τ 是给定的分位数。决策变量是 y ^ i \hat{y}_i y^i,可以证明,使上面表达式最小化的 y ^ i \hat{y}_i y^i 就是给定分位数 τ \tau τ 对应的分位点(将上面式子转化为连续密度函数的积分,然后求一阶导数即可证明)。

    上式也可以简写成:
    min ⁡ y ^ i ∑ [ τ ( y i − y ^ i ) + + ( 1 − τ ) ( y ^ i − y i ) + ] \min_{\hat{y}_i}\quad \sum \left[\tau(y_i-\hat{y}_i)^++(1-\tau)(\hat{y}_i-y_i)^+\right] y^imin[τ(yiy^i)++(1τ)(y^iyi)+]

    2. 分位数回归的求解原理

    为了求出分位数的回归方程,假设 y ^ i = X i β \hat{y}_i=\bm{X_i \beta} y^i=Xiβ,那么求解的目标函数转化为:

    arg min ⁡ β ∈ R k ∑ [ τ ( y i − X i β ) + + ( 1 − τ ) ( X i β − y i ) + ] \argmin_{\bm{\beta}\in\mathbb{R}^k} \sum \left[\tau(y_i-\bm{X_i \beta})^++(1-\tau)(\bm{X_i \beta}-y_i)^+\right] βRkargmin[τ(yiXiβ)++(1τ)(Xiβyi)+]

    决策变量为 k k k 维回归方程的参数向量 β \bm{\beta} β。在实际的求解中,将上式转化为一个线性规划问题,引入两个虚拟变量 u i + u_i^+ ui+ u i − u_i^- ui

    arg min ⁡ β ∈ R k ∑ [ τ u i + + ( 1 − τ ) u i − ] s . t . i = 1 , 2 , … , n y i = X i β + u i + − u i − u i + ≥ 0 , u i − ≥ 0 \begin{aligned} &\argmin_{\bm{\beta}\in\mathbb{R}^k}\quad && \sum \left[\tau u_i^++(1-\tau)u_i^-\right]\\ &s.t. && i=1,2,\dots, n\\ &&& y_i=\bm{X_i \beta}+u_i^+-u_i^-\\ &&& u_i^+\geq 0, u_i^-\geq 0 \end{aligned} βRkargmins.t.[τui++(1τ)ui]i=1,2,,nyi=Xiβ+ui+uiui+0,ui0

    然后用单纯形法或内点法求解,就能得出分位数回归方程(python 与 R 软件求出的分位数回归方程可能略微不同,因为求解方法不一样, python 使用了迭代的加权最小二乘法求解)。

    3 python 分位数回归

    分位数回归要用到 statsmodels,下面的代码得到分位数为 0.6 的回归方程,并画图:

    import statsmodels.formula.api as smf
    import pandas as pd
    import numpy as np
    import statsmodels.api as sm
    import matplotlib.pyplot as plt
    
    
    data = sm.datasets.engel.load_pandas().data
    mod = smf.quantreg('foodexp ~ income', data)
    res = mod.fit(q=0.6)
    print(res.summary())
    
    plt.scatter(data['income'], data['foodexp'])
    xx = np.arange(min(data['income']), max(data['income']))
    yy = [i*res.params['income'] + res.params['Intercept'] for i in xx]
    plt.plot(xx, yy, color='red')
    plt.show()
    

    输出结果:

    在这里插入图片描述
    在这里插入图片描述
    回归方程就是上面的红线,它将 40% 的数据分割在上面,60% 的数据分割在下面。

    转载于个人公众号:Python 统计分析与数据科学

    在这里插入图片描述

    展开全文
  • 分位数回归及其Python源码

    千次阅读 2020-12-06 14:51:03
    分位数回归及其Python源码天朗气清,惠风和畅。赋闲在家,正宜读书。前人文章,不得其解。代码开源,无人注释。你们不来,我行我上。废话少说,直入主题。o( ̄︶ ̄)o我们要探测自变量 与因变量 的关系,最简单的方法...
  • 分位数回归森林

    2020-12-06 14:51:01
    分位数回归森林(Quantile Regression Forests),一般回归模型预测均值,但该算法预测数据的分布。它可以用来预测给定输入的价格分布,例如,给定一些属性,汽车价格分布的第25和75百分位是多少。大多数预测器在预测...
  • 分位数回归,分位数回归模型,Python源码.zip
  • Python中statsmodels包中的QuantReg使用如下代码中所示的数据,得到了与R中截然不同的结果。在我分别用Python和R测试了STACKLOSS数据,结果是一样的。我想知道数据本身是否在Python中引起了一些问题,或者可能两种...
  • 我遵循StatsModels示例here来绘制分位数回归线.只需对我的数据稍作修改,该示例效果很好,生成此绘图(请注意,我已修改代码以仅绘制0.05,0.25,0.5,0.75和0.95分位数):但是,我想绘制OLS拟合和相应分位数的二阶多项式...
  • 假设房地产分析师想要根据家庭、...这称为预测区间,产生它们的一般方法称为分位数回归。在这篇文章中,我将描述这个问题是如何正式化,如何采用六种线性,基于决策树和深度学习的方法中实现它(在Python中,这是Ju...
  • 你需要弄清楚某个点是在95%分位数线上还是在5%分位数线之下。这可以使用叉积来实现,请参见this answer以获得简单的实现。在在您的例子中,您需要将分位数线上下的点组合起来,可能是在一个遮罩中。在下面是一个例子...
  • 这种理论也可以在预测统计中为我们服务,这正是分位数回归的意义所在——估计中位数(或其他分位数)而不是平均值。 通过选择任何特定的分位数阈值,我们既可以缓和异常值,也可以调整错误的正/负权衡。我们还可以...
  • 摘要 贝叶斯回归分位数在最近的文献中受到广泛关注,本文实现了贝叶斯系数估计和回归分位数(RQ)中的变量选择,带有lasso和自适应...自引入以来,分位数回归一直是理论界非常关注的话题,也在许多研究领域得到了.
  • 最后,将投资组合的损益从小到大排序,得到损益分布,通过给定置信度下的分位数求出VaR。 简单历史模拟法是将目前投资组合的价格按历史时间段的收益率重新抽样,计算组合的损失以及VaR。计算1-day VaR可抽取历史日...
  • 分位数 VAR 模型估计:自回归分布滞后模型脉冲响应函数计算各分位点脉冲图绘制代码实现原理使用pyqt5生成 GUI 界面使用statsmodels进行分位数回归使用pandas将结果保存在 excel 文件Display主窗口界面主要包括:工具...
  • 多元线性模型的分位数回归

    千次阅读 2021-04-11 10:26:42
    分位数回归学习笔记一、为什么要使用分位数回归?二、分位数回归基本模型三、分位数回归估计--单纯形法1.损失函数2.目标函数3.算法推导4.实际案例分析与python代码 一、为什么要使用分位数回归?    ...
  • 分位数回归(Quantile Regression)

    千次阅读 2020-06-29 14:56:31
    数据采用分位数回归 在执行回归分析时,仅对问题进行数值预测还不够,您还需要表达您对该预测的信心。例如,如果您正在查看特定市场中房屋的价格,并且您的模型预测房屋的售价为262,458.45美元,那么您对模型的预测...
  • 分位数回归基准 BenchOpt是一个软件包,用于简化优化算法的比较并使其更加透明和可重现。 此基准专用于L1正则化分位数回归问题: \min_{w} \frac{1}{n} \sum_{i=1}^{n} \text{pinball}(y_i, x_i^\top w) + \frac{\...
  • 线性回归模型属于经典的统计学模型,该模型的应用场景是根据已知的变量(自变量)来预测某个连续的数值变量(因变量)。例如,餐厅根据每天的营业数据(包括菜谱价格、就餐人数、预定人数、特价菜折扣等)预测就餐...
  • 机器学习:Python实现随机森林回归

    千次阅读 2020-12-22 20:49:53
    随机森林(Random forests)或随机决策森林(Random decision forests)是一种用于分类、回归和其他任务的集成学习方法,通过在训练时构建大量决策树并输出作为单个树的类(分类)或平均预测(回归)模式的类(class)来操作。...
  • 使用python进行面板回归,顾名思义就是,使用python这种语言进行面板数据的回归。 一、面板数据解释: 1、面板数据,即Panel Data,也叫“平行数据”,是指在时间序列上取多个截面,在这些截面上同时选取样本观测...
  • 基于R语言的分位数回归(quantile regression)

    万次阅读 多人点赞 2017-12-18 17:45:21
    分位数回归(quantile regression) 这一讲,我们谈谈分位数回归的知识,我想大家传统回归都经常见到。分位数回归可能大家见的少一些,其实这个方法也很早了,大概78年代就有了,但是那个时候这个理论还不完善。到...
  • 本文采用python sklearn库中,作为quantile regression的示例代码。以下为详细解析: import numpy as np import matplotlib.pyplot as plt from sklearn.ensemble import GradientBoostingRegressor %matplotlib ...
  • Python实现随机森林回归

    万次阅读 多人点赞 2020-11-02 13:19:38
    这里,将介绍如何在Python中构建和使用Random Forest回归,而不是仅仅显示代码,同时将尝试了解模型的工作原理。 1.1 随机森林概述 随机森林是一种基于集成学习的监督式机器学习算法。集成学习是一种学习类型,可以...
  • 电气论文实现:深度学习分位数回归实现电力负荷区间预测(有代码、有数据)。

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 4,104
精华内容 1,641
关键字:

分位数回归python

友情链接: HttpAbstractParamBean.rar