Quantilregression als lineares Programmierproblem formulieren?

8

Wie formuliere ich die Quantilregression als lineares Programmierproblem? Wenn ich mir das mittlere Quantilproblem anschaue, weiß ich, dass es das ist

minimize i=1n|β0+Xiβ1Yi|transforms into minimize i=1neis.t.eiβ0+Xiβ1Yiei(β0+Xiβ1Yi)
aber wie transformiere ich die Minimierung für andere Quantile?

machazthegamer
quelle
Einige relevante R-Codes finden Sie hier .
kjetil b halvorsen

Antworten:

14

Sie verwenden den Quantilregressionsschätzer

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

Dabei ist τ(0,1) eine Konstante, nach der das Quantil geschätzt werden muss, und die Funktion ρτ(.) ist definiert als

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

Um den Zweck von betrachten Sie zunächst, dass die Residuen als Argumente verwendet werden, wenn diese als . Die Summe im Minimierungsproblem kann daher umgeschrieben werden alsρτ(.)ϵi=yixiθ

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

so dass positive Residuen, die mit der Beobachtung oberhalb der vorgeschlagenen Quantilregressionslinie assoziiert sind, das Gewicht von während negative Residuen, die mit Beobachtungen unterhalb der vorgeschlagenen Quantilregressionslinie assoziiert sind werden mit gewichtet .yixiθτyixiθ(1τ)

Intuitiv:

Mit positive und negative Residuen mit dem gleichen Gewicht "bestraft" und eine gleiche Anzahl von Beobachtungen befindet sich optimal über und unter der "Linie", so dass die Linie die mittlere Regression ist "Linie".τ=0.5xiβ^

Wenn jedes positive Residuum 9-mal so gewichtet wie das eines negativen Residuums mit einem Gewicht von und somit optimal für jede Beobachtung über der "Linie" ungefähr 9 unter der Linie platziert werden. Daher repräsentiert die "Linie" das 0,9-Quantil. (Für eine genaue Aussage siehe THM. 2.2 und Korollar 2.1 in Koenker (2005) "Quantile Regression")τ=0.91τ=0.1xiβ^

Die beiden Fälle sind in diesen Darstellungen dargestellt. Linkes Feld und rechtes Feld .τ=0.5τ=0.9

Rho-Funktion tau = 0,5 und tau = 0,9

Lineare Programme werden überwiegend mit dem Standardformular analysiert und gelöst

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

Um zu einem linearen Programm in Standardform zu gelangen, besteht das erste Problem darin, dass in einem solchen Programm (1) alle Variablen, über die die Minimierung durchgeführt wird, positiv sein sollten. Um dies zu erreichen, werden Residuen unter Verwendung von Slack-Variablen in positive und negative Teile zerlegt:z

ϵi=uivi

wobei positiver Teil und ist der negative Teil. Die Summe der Residuen, denen durch die Prüffunktion Gewichte zugewiesen wurden, wird dann als gesehenui=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,

Dabei ist und und Vektor alle Koordinaten gleich .u=(u1,...,uN)v=(v1,...,vN)1NN×11

Die Residuen müssen die Bedingungen erfüllen, dieN

yixiθ=ϵi=uivi

Dies führt zur Formulierung als lineares Programm

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

wie in Koenker (2005) "Quantile Regression", Seite 10, Gleichung (1.20) angegeben.

Es fällt jedoch auf, dass immer noch nicht darauf beschränkt ist, positiv zu sein, wie es im linearen Programm auf Standardform (1) erforderlich ist. Daher wird wieder eine Zerlegung in einen positiven und einen negativen Teil verwendetθR

θ=θ+θ

wobei wiederum ein positiver Teil und ein negativer Teil ist. Die Bedingungen können dann geschrieben werden alsθ+=max(0,θ)θ=max(0,θ)N

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

Dabei ist .IN=diag{1N}

Definieren Sie als nächstes und die Entwurfsmatrix der Daten für unabhängige Variablen als gespeichert werdenb:=yX

X:=[x1xN]

So schreiben Sie die Einschränkung neu:

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

Definieren Sie die Matrix(N×2K+2N)

A:=[X,X,IN,IN]
und und als Variablen ein, über die minimiert werden soll, damit sie Teil von , um zu erhaltenθ+θz

b=A[θ+θuv]=Az

Da und das Minimierungsproblem nur durch die Bedingung beeinflussen, muss a der Dimension als Teil des Koeffizientenvektors , der entsprechend definiert werden kann alsθ+θ02K×1c

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

Dadurch wird sichergestellt, dasscz=0(θ+θ)=0+τ1Nu+(1τ)1Nv=i=1Nρτ(ϵi).

Daher werden dann und definiert und das in angegebene Programm vollständig spezifiziert.c,Ab(1)

Dies wird wahrscheinlich am besten anhand eines Beispiels verdaut. Um dies in R zu lösen, verwenden Sie das Paket quantreg von Roger Koenker. Hier sehen Sie auch, wie Sie das lineare Programm einrichten und mit einem Löser für lineare Programme lösen:

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 für den Präsidenten
quelle
Haben
Ich werde erweitern, habe Geduld, Änderungen kommen :)
Jesper für Präsident
Stellen Sie jetzt eine Folgefrage, wenn etwas ausgearbeitet werden muss.
Jesper für Präsident
4

Ich habe den Code von Jesper Hybel in Python mit cvxopt umgeschrieben. Ich poste es hier, falls jemand anderes dies auch in Python benötigt.

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)
```
Kumpel Uzsoki
quelle
Vielen Dank, Ihre Übersetzung nach Python hat mir sehr geholfen. Gund herscheinen nicht im ursprünglichen R-Code oder in Jespers Bericht. Sind diese Artefakte, wie CVXOPT die Formulierung von Problemen erfordert, oder sind sie in einem LP-Solver enthalten? In meinem Fall bin ich auf ein Hindernis gestoßen, das versucht hat, mit N = 50.000 weiterzulaufen. Gwird in diesem Fall eine riesige quadratische Matrix. Ich denke darüber nach, einen verteilten LP-Solver wie Spark zu verwenden, aber vielleicht ist die Verwendung von LP für die Quantilregression auf dieser Datenskala einfach nicht nachvollziehbar.
peace_within_reach
Follow-up zu meinem früheren Kommentar: Ich konnte die quantreg rqRoutine mit 15 Millionen Datenzeilen ausführen . Ich war beeindruckt, dass jede auf linearer Programmierung basierende Methode so viele Daten verarbeiten kann. In meinem Fall (Schätzung sehr hoher Quantile) müssen jedoch noch mehr Daten verwendet werden. Ich habe rqDrosseln gefunden, wenn ich 20 Millionen oder mehr Reihen benutze. Der Fehler liegt darin long vectors are not supported in .C, dass die Vektoren zu lang werden. Für jeden in der gleichen Situation ist die beste Software, die ich für die Quantilregression auf Big Data gefunden habe, Microsofts LightGBM (Gradient Boosting)
peace_within_reach
1
Ja, G und h gelten nur für die CVXOPT-Formulierung. Sie finden sie auch, wenn Sie die CVXOPT-Dokumentation lesen. Der Ansatz in CVXOPT besteht darin, dass Gleichheitsbeschränkungen in einer Matrix (A) und Ungleichheitsbeschränkungen in einer anderen (G) gespeichert werden. Gleiches gilt für die rechte Matrix (h).
Kumpel Uzsoki
2
Wenn Sie Leistung benötigen, ist es wichtig, solver = 'glpk' zu verwenden - die Geschwindigkeit ändert sich stark. Es kann auch eine gute Idee sein, mit anderen Lösern zu experimentieren, um festzustellen, ob sie zu schnelleren Ergebnissen führen.
Kumpel Uzsoki