Execute regressão linear, mas force a solução a passar por alguns pontos de dados específicos

13

Eu sei como executar uma regressão linear em um conjunto de pontos. Ou seja, eu sei como ajustar um polinômio de minha escolha a um determinado conjunto de dados (no sentido LSE). No entanto, o que não sei é como forçar minha solução a passar por alguns pontos específicos de minha escolha. Já vi isso ser feito antes, mas não me lembro como o procedimento foi chamado, muito menos como foi feito.

Como um exemplo muito simples e concreto, digamos que tenho 100 pontos espalhados no plano xy e opto por ajustar um polinômio de qualquer ordem através deles. Eu sei como executar essa regressão linear muito bem. No entanto, digamos que eu quero 'força' a minha solução, que passar por, digamos, três dos meus pontos de dados em coordenadas x , e , (e suas correspondentes coordenadas y claro).x = 19 x = 89x=3x=19x=89

Como é chamado esse procedimento geral, como é feito e há alguma armadilha específica da qual preciso estar ciente?

Editar:

Gostaria de acrescentar que estou procurando uma maneira concreta de fazer isso. Eu escrevi um programa que realmente faz a regressão linear de uma de duas maneiras, invertendo a matriz de covariância diretamente ou através da descida do gradiente. O que estou perguntando é: como, exatamente, passo a passo, modifico o que fiz, de modo a forçar a solução polinomial a passar por pontos específicos?

Obrigado!

Spacey
fonte
Por que você o chama de "linear" se está usando um polinômio? Cada ponto que você deseja que ele passe é uma restrição que reduzirá seu grau de liberdade. Você pode usar um algoritmo de otimização restrito.
22253 curious_cat
3
É linear porque você está encontrando coeficientes para uma combinação linear . Por exemplo, se você deseja ajustar seus dados a um cúbico, está encontrando os coeficientes (os c 's) de y=c0 0+c1x+c2x2+c3x3 .
Spacey
1
@ Mohammad: Uma outra maneira de aproximar o que você deseja seria usar uma solução de mínimos quadrados ponderados e atribuir pesos muito grandes aos pontos pelos quais você deseja que a linha de regressão passe. Isso deve forçar a solução a passar muito de perto para os pontos que você escolher.
Jason R
@ JasononR Bom vê-lo aqui. Sim, o WLS é realmente um candidato interessante. Fui com a resposta dos whubers por causa da fatoração polinomial inteligente e porque mantém a estrutura do erro de maneira adequada.
Spacey

Respostas:

18

O modelo em questão pode ser escrito

y=p(x)+(xx1)(xxd)(β0+β1x++βpxp)+ε

onde é um polinômio de grau d - 1 passando por pontos predeterminados ( x 1 , y 1 ) , , ( x d , y d ) e ε é aleatório. (Use o polinômio de interpolação de Lagrange .) Escrita ( x - x 1 ) ( x - x d ) = rp(xi)=yid1(x1,y1),,(xd,yd)ε nos permite reescrever este modelo como(xx1)(xxd)=r(x)

yp(x)=β0r(x)+β1r(x)x+β2r(x)x2++βpr(x)xp+ε,

que é um problema de regressão múltipla OLS padrão com a mesma estrutura de erro que o original , onde as variáveis independentes são os quantidades de r ( x ) x i , i = 0 , 1 , ... , p . Simplesmente calcule essas variáveis ​​e execute seu software de regressão familiar , evitando que ele inclua um termo constante. As advertências usuais sobre regressões sem termo constante se aplicam; em particular, o R 2 pode ser artificialmente elevado; as interpretações usuais não se aplicam.p+1r(x)xi, i=0,1,,pR2

(Na verdade, a regressão através da origem é um caso especial desta construção, em que , ( x 1 , y 1 ) = ( 0 , 0 ) , e p ( x ) = 0 , de modo que o modelo é y = β 0 x + + β p x p + 1 + ε . )d=1(x1,y1)=(0,0)p(x)=0y=β0x++βpxp+1+ε.


Aqui está um exemplo trabalhado (em R)

# Generate some data that *do* pass through three points (up to random error).
x <- 1:24
f <- function(x) ( (x-2)*(x-12) + (x-2)*(x-23) + (x-12)*(x-23) )  / 100
y0 <-(x-2) * (x-12) * (x-23) * (1 + x - (x/24)^2) / 10^4  + f(x)
set.seed(17)
eps <- rnorm(length(y0), mean=0, 1/2)
y <- y0 + eps
data <- data.frame(x,y)

# Plot the data and the three special points.
plot(data)
points(cbind(c(2,12,23), f(c(2,12,23))), pch=19, col="Red", cex=1.5)

# For comparison, conduct unconstrained polynomial regression
data$x2 <- x^2
data$x3 <- x^3
data$x4 <- x^4

fit0 <- lm(y ~ x + x2 + x3 + x4, data=data)
lines(predict(fit0), lty=2, lwd=2)

# Conduct the constrained regressions
data$y1 <- y - f(x)
data$r <- (x-2)*(x-12)*(x-23)
data$z0 <- data$r
data$z1 <- data$r * x
data$z2 <- data$r * x^2

fit <- lm(y1 ~ z0 + z1 + z2 - 1, data=data)
lines(predict(fit) + f(x), col="Red", lwd=2)

Plot

Os três pontos fixos são mostrados em vermelho sólido - eles não fazem parte dos dados. O ajuste mínimo de quadrados do polinômio de quarta ordem sem restrições é mostrado com uma linha pontilhada preta (possui cinco parâmetros); o ajuste restrito (da ordem cinco, mas com apenas três parâmetros livres) é mostrado com a linha vermelha.

Inspecionar a saída dos mínimos quadrados ( summary(fit0)e summary(fit)) pode ser instrutivo - deixo isso para o leitor interessado.

whuber
fonte
whuber, isso é interessante ... Eu mentiria se dissesse que ainda o entendi completamente, mas estou digerindo enquanto falamos. Se entendi corretamente, aqui estou basicamente resolvendo os s como de costume, mas eles estão sendo multiplicados por r ( x ) x i , em vez de apenas x i s como antes, sim? Se isso estiver correto, como você está calculando r ( x ) exatamente? Obrigado. βr(x)xixir(x)
Spacey
Eu adicionei um exemplo trabalhado, Mohammad.
whuber
Oh, perfeito. Eu estudarei isso Usando seu exemplo, ainda seria possível forçar o poli a passar por pontos que fazem parte dos dados, certo?
Spacey
Absolutamente isso pode ser feito: mas seja duplamente cauteloso ao interpretar os valores-p ou quaisquer outras estatísticas, porque agora suas restrições são baseadas nos próprios dados.
whuber
Sua postagem me deixou acordada ontem à noite. Eu me ensinei o LIP. (LIP é interessante. É como uma decomposição de Fourier, mas com polis).
Spacey
9

(xi,yi)xEu de todo x-value e yEu de todo y-valor. Agora o ponto está na origem do plano de coordenadas. Você simplesmente ajusta uma linha de regressão enquanto suprime a interceptação (forçando a interceptação a ser (0,0). Como essa é uma transformação linear, você pode facilmente transformar tudo depois posteriormente, se desejar.

Se você deseja forçar uma linha a passar por dois pontos em um plano XY, isso também é bastante fácil. Quaisquer dois pontos podem ser ajustados com uma linha. Você pode usar a fórmula da inclinação dos pontos para calcular sua inclinação e, em seguida, usar um dos pontos, a inclinação e a equação de uma linha para encontrar a interceptação.

Observe que talvez não seja possível ajustar uma linha reta através de três pontos em um plano de coordenadas. No entanto, podemos garantir que eles se encaixam perfeitamente com uma parábola (ou seja, usando os doisX e X2) Também há álgebra para isso, mas à medida que avançamos, pode ser mais fácil ajustar um modelo ao software, incluindo apenas esses três (mais) pontos no conjunto de dados. Da mesma forma, você pode obter a linha reta que melhor se aproxima desses três pontos, ajustando um modelo que tenha acesso apenas a esses três pontos.


Sinto-me compelido a mencionar, a essa altura, que isso pode não ser uma boa coisa a se fazer (a menos que sua teoria forneça razões muito sólidas para fazê-lo). Você também pode analisar a regressão bayesiana , onde pode permitir que seu modelo encontre a melhor combinação de informações em seus dados e algumas informações anteriores (que você pode usar para influenciar fortemente sua interceptação para zero, por exemplo, sem forçando).

gung - Reinstate Monica
fonte
1
Gung, thanks for your answer. I have modified my question a little. I did not know about Bayesian Regression, but will take a look at it. I am afraid I do not completely understand how, concretely, from an algorithmic view point, the one point and two point case you mentioned. Specifically, for the one point, I understand removing and re-adding xi and yi to each point before and after a block, but I do not understand how to do that block exactly. For the 2-point case, I am afraid I do not understand at all what to do there. Thanks.
Spacey
2
Although throwing in three more points and weighting them (a la Glen_b's answer) could create such a fit, interpreting any of the statistical output would be problematic: some adjustments would be needed.
whuber
6

To add a little extra information to @gung's excellent coverage of the linear case, in the higher order polynomial case there are several ways you could do it either exactly or approximately (but pretty much as accurately as you need).

Primeiro, observe que os graus de liberdade do polinômio (ou mesmo de qualquer função ajustada) devem ser pelo menos tão grandes quanto o número de pontos "conhecidos". Se os graus de liberdade forem iguais, você não precisará dos dados, pois a curva está completamente determinada. Se houver mais pontos "conhecidos", você não poderá resolvê-lo (a menos que todos estejam exatamente no mesmo polinômio do grau especificado, caso em que qualquer subconjunto de tamanho adequado será suficiente). Daqui em diante, falarei sobre quando o polinômio tem mais df do que os pontos conhecidos (como um cúbico - com 4df - e três pontos conhecidos, para que o cúbico não seja sobredeterminado por pontos conhecidos nem completamente determinado por eles) .

1) "a curva deve passar por esse ponto" é uma restrição linear nos parâmetros, resultando em estimativa restrita ou em mínimos quadrados restritos (embora ambos os termos possam incluir outras coisas além de restrições lineares, como restrições de positividade). Você pode incorporar restrições lineares por

  (a) reformular a parametrização para incluir implicitamente cada restrição, resultando em um modelo de ordem inferior.

  (b) usando ferramentas padrão que podem incorporar restrições lineares nos parâmetros de ajuste de mínimos quadrados. (geralmente por meio de algo como a fórmula fornecida no link acima)

2) Another way is via weighted regression. If you give the known points sufficiently large weight, you can get essentially the same fit as in (1). This is often readily implemented, can be substantially quicker than reparameterizing, and can be done in packages that don't offer constrained fitting.

All of @gung's caveats apply

Glen_b -Reinstate Monica
fonte
Glen_b, eu não havia considerado a regressão ponderada. Pode ser o caminho a seguir. Eu coloquei na minha lista de tarefas. Acredito que posso me ensinar isso sem incidentes. Em relação a (1), você poderia gentilmente expandir esse aspecto da reparmaterização? Além disso, como você chama isso que estou tentando fazer, onde forço o polinômio a passar por certos pontos? Parte do problema é que eu não sei o que procurar no Google. Se eu souber como isso foi chamado, talvez eu possa aumentar o que você está dizendo com material on-line. Obrigado.
Spacey
Veja minhas edições acima, que incluem alguns termos de pesquisa e um link com mais alguns detalhes.
Glen_b -Reinstala Monica
2
A regressão ponderada com +1 é uma boa ideia. Alguns ajustes das estatísticas de saída, como estimativas de erro do RMS, podem ser necessários.
whuber
@whuber +1, de fato, se as estatísticas (como s2, F, R2... erros padrão etc.) destinam-se a se relacionar apenas aos pontos "não conhecidos" (o que provavelmente é o que se deseja), além das estimativas de parâmetros e dos valores ajustados, as estatísticas brutas que saem quase todas estarão erradas . Eu originalmente digitei uma frase relacionada a isso, mas parece que a excluí antes de postar; é importante mencionar isso.
Glen_b -Reinstala Monica,
Thanks for your answer Glen_b, although I have accepted @whuber 's, I still learned a lot from yours.
Spacey