Formular regressão quantílica como problema de programação linear?

8

Como formular regressão quantílica como um problema de programação linear? Ao olhar para o problema mediano dos quantis, sei que é

minimize i=1n|β0+Xiβ1Yi|transforms into minimize i=1neis.t.eiβ0+Xiβ1Yiei(β0+Xiβ1Yi)
mas como transformar a minimização para outros quantis?

machazthegamer
fonte
Algum código R relevante pode ser encontrado aqui .
Kjetil b halvorsen

Respostas:

14

Você usa o estimador de regressão quantil

β^(τ):=argminθRKi=1Nρτ(yixiθ).

onde τ(0,1) é constante escolhido de acordo com o qual o quantil precisa ser estimado e a função ρτ(.) é definida como

ρτ(r)=r(τI(r<0)).

Para ver o objetivo de Considere primeiro que ele leva os resíduos como argumentos, quando são definidos como . A soma do problema de minimização pode, portanto, ser reescrita comoρτ(.)ϵi=yixiθ

i=1Nρτ(ϵi)=i=1Nτ|ϵi|I[ϵi0]+(1τ)|ϵi|I[ϵi<0]

de modo que resíduos positivos associados à observação acima da linha de regressão quantílica sugerida recebem o peso de enquanto resíduos negativos associados à observação abaixo da linha de regressão quantílica sugerida são ponderados com .yixiθτyixiθ(1τ)

Intuitivamente:

Com resíduos positivos e negativos são "punidos" com o mesmo peso e um número igual de observação está acima e abaixo da "linha" no ideal, de modo que a linha seja a regressão mediana "linha".τ=0.5xiβ^

Quando cada resíduo positivo é ponderado 9 vezes o de um resíduo negativo com peso e, portanto, é ideal para todas as observações acima da "linha" aproximadamente 9 ser colocado abaixo da linha. Portanto, a "linha" representa o quantil 0,9. (Para uma declaração exata disso, consulte THM. 2.2 e Corollary 2.1 em Koenker (2005) "Regressão quantílica")τ=0.91τ=0.1xiβ^

Os dois casos são ilustrados nessas parcelas. Painel esquerdo e painel direito .τ=0.5τ=0.9

função rho tau = 0,5 e tau = 0,9

Os programas lineares são predominantemente analisados ​​e resolvidos usando o formulário padrão

(1)  minz  cz  subject to Az=b,z0

Para chegar a um programa linear na forma padrão, o primeiro problema é que em um programa (1) todas as variáveis ​​sobre as quais a minimização é realizada devem ser positivas. Para isso, os resíduos são decompostos em partes positivas e negativas usando variáveis ​​de folga:z

ϵi=uivi

onde parte positiva e é a parte negativa. A soma dos pesos residuais atribuídos pela função de verificação é então vista comoui=max(0,ϵi)=|ϵi|I[ϵi0]vi=max(0,ϵi)=|ϵi|I[ϵi<0]

i=1Nρτ(ϵi)=i=1Nτui+(1τ)vi=τ1Nu+(1τ)1Nv,

onde e e é vetor todas as coordenadas igual a .u=(u1,...,uN)v=(v1,...,vN)1NN×11

Os resíduos devem satisfazer as restrições queN

yixiθ=ϵi=uivi

Isso resulta na formulação como um programa linear

minθRK,uR+N,vR+N{τ1Nu+(1τ)1Nv|yi=xiθ+uivi,i=1,...,N},

como indicado na equação de Koenker (2005) "Regressão quantílica", página 10 (1.20).

No entanto, é perceptível que ainda não está restrito a ser positivo, conforme exigido no programa linear no formulário padrão (1). Portanto, novamente se decompõe em parte positiva e negativa.θR

θ=θ+θ

onde novamente é parte positiva e é parte negativa. As restrições podem então ser escritas comoθ+=max(0,θ)θ=max(0,θ)N

y:=[y1yN]=[x1xN](θ+θ)+INuINv,

onde .IN=diag{1N}

Em seguida, defina e a matriz de design armazenando dados em variáveis ​​independentes comob:=yX

X:=[x1xN]

Para reescrever restrição:

b=X(θ+θ)+INuINv=[X,X,IN,IN][θ+θuv]

Definir a matriz(N×2K+2N)

A:=[X,X,IN,IN]
e introduza e como variáveis ​​sobre as quais minimizar, para que façam parte de para obterθ+θz

b=A[θ+θuv]=Az

Como e afetam apenas o problema de minimização através da restrição a da dimensão deve ser introduzida como parte do vetor coeficiente que pode ser definido adequadamente comoθ+θ02K×1c

c=[0τ1N(1τ)1N],

garantindo assim quecz=0(θ+θ)=0+τ1Nu+(1τ)1Nv=i=1Nρτ(ϵi).

Portanto e são então definidos e o programa, conforme indicado em completamente especificado.c,Ab(1)

Provavelmente é melhor digerido usando um exemplo. Para resolver isso em R, use o pacote quantreg de Roger Koenker. Aqui está também uma ilustração de como configurar o programa linear e resolver com um solucionador de programas lineares:

base=read.table("http://freakonometrics.free.fr/rent98_00.txt",header=TRUE)
attach(base)
library(quantreg)
library(lpSolve)
tau <- 0.3


# Problem (1) only one covariate
X <- cbind(1,base$area)
K <- ncol(X)
N <- nrow(X)

A <- cbind(X,-X,diag(N),-diag(N))
c <- c(rep(0,2*ncol(X)),tau*rep(1,N),(1-tau)*rep(1,N))
b <- base$rent_euro
const_type <- rep("=",N)

linprog <- lp("min",c,A,const_type,b)
beta <- linprog$sol[1:K] -  linprog$sol[(1:K+K)]
beta
rq(rent_euro~area, tau=tau, data=base)


# Problem (2) with 2 covariates
X <- cbind(1,base$area,base$yearc)
K <- ncol(X)
N <- nrow(X)

A <- cbind(X,-X,diag(N),-diag(N))
c <- c(rep(0,2*ncol(X)),tau*rep(1,N),(1-tau)*rep(1,N))
b <- base$rent_euro
const_type <- rep("=",N)

linprog <- lp("min",c,A,const_type,b)
beta <- linprog$sol[1:K] -  linprog$sol[(1:K+K)]
beta
rq(rent_euro~ area + yearc, tau=tau, data=base)
Jesper para Presidente
fonte
você se importa de explicar essa minimização em palavras ... Eu não consigo entender direito #
machazthegamer
Vou expandir a ter paciência, as edições estão chegando :) #
247 Jesper for President
Faça a pergunta de acompanhamento agora se algo precisar de elaboração.
Jesper for President
4

Reescrevi o código de Jesper Hybel em Python usando cvxopt. Estou postando aqui, caso alguém também precise disso em Python.

import pandas as pd
import io
import requests
import numpy as np
url="http://freakonometrics.free.fr/rent98_00.txt"
s=requests.get(url).content
base=pd.read_csv(io.StringIO(s.decode('utf-8')), sep='\t')


tau = 0.3

from cvxopt import matrix, solvers

X = pd.DataFrame(columns=[0,1])
X[1] = base["area"] #data points for independent variable area
X[2] = base["yearc"] #data points for independent variable year
X[0] = 1 #intercept

K = X.shape[1]
N = X.shape[0]

# equality constraints - left hand side

A1 = X.to_numpy() # intercepts & data points - positive weights
A2 = X.to_numpy() * - 1 # intercept & data points - negative weights
A3 = np.identity(N) # error - positive
A4 = np.identity(N)*-1 # error - negative

A = np.concatenate((A1,A2,A3,A4 ), axis= 1) #all the equality constraints 

# equality constraints - right hand side
b = base["rent_euro"].to_numpy()

#goal function - intercept & data points have 0 weights
#positive error has tau weight, negative error has 1-tau weight
c = np.concatenate((np.repeat(0,2*K), tau*np.repeat(1,N), (1-tau)*np.repeat(1,N) ))

#converting from numpy types to cvxopt matrix

Am = matrix(A)
bm = matrix(b)
cm = matrix(c)

# all variables must be greater than zero
# adding inequality constraints - left hand side
n = Am.size[1]
G = matrix(0.0, (n,n))
G[::n+1] = -1.0

# adding inequality constraints - right hand side (all zeros)
h = matrix(0.0, (n,1))

#solving the model
sol = solvers.lp(cm,G,h,Am,bm, solver='glpk')

x = sol['x']

#both negative and positive components get values above zero, this gets fixed here
beta = x[0:K] - x[K :2*K]

print(beta)
```
Mate Uzsoki
fonte
Obrigado, sua tradução para Python me ajudou muito. Ge hnão apareça no código R original ou no artigo de Jesper. Esses artefatos de como o CVXOPT exigem a formulação de problemas ou estão implícitos em qualquer solucionador de LP? No meu caso, encontrei um obstáculo tentando continuar com N = 50.000. Gtorna-se uma enorme matriz quadrada neste caso. Estou pensando em usar um solucionador de LP distribuído como o Spark, mas talvez o uso de LP para regressão quantil nessa escala de dados não seja tratável.
peace_within_reach
Seguimento do meu comentário anterior: Consegui executar a quantreg rqrotina com 15 milhões de linhas de dados. Fiquei impressionado que qualquer método baseado em programação linear pudesse lidar com tantos dados. No entanto, no meu caso (estimando quantis muito altos), é necessário usar ainda mais dados do que isso. Encontrei rqestrangulamentos quando uso 20 milhões ou mais de linhas. O erro é long vectors are not supported in .Cque ocorre porque os vetores se tornam muito longos. Para qualquer pessoa na mesma situação, o melhor software que eu encontrei para regressão quantílica sobre big data é LightGBM da Microsoft (gradiente aumentando)
peace_within_reach
11
Sim, G e h são apenas para a formulação CVXOPT. Você também os encontrará se estiver lendo a documentação do CVXOPT. A abordagem no CVXOPT é que restrições de igualdade são armazenadas em uma matriz (A) e restrições de desigualdade em outra (G). O mesmo vale para a matriz do lado direito (h).
Mate Uzsoki 28/09/19
2
Se você precisar de desempenho, é importante usar o solver = 'glpk' - isso muda muito a velocidade. Também pode ser uma boa ideia experimentar outros solucionadores para ver se eles oferecem resultados mais rápidos.
Mate Uzsoki 28/09/19