Effectuer une régression linéaire, mais forcer la solution à passer par certains points de données particuliers

14

Je sais comment effectuer une régression linéaire sur un ensemble de points. Autrement dit, je sais comment adapter un polynôme de mon choix, à un ensemble de données donné (au sens LSE). Cependant, ce que je ne sais pas, c'est comment forcer ma solution à passer par certains points particuliers de mon choix. J'ai déjà vu cela se faire auparavant, mais je ne me souviens pas de la façon dont la procédure a été appelée, et encore moins de la façon dont elle a été effectuée.

Comme exemple très simple et concret, disons que j'ai 100 points dispersés sur le plan xy, et je choisis de faire passer un polynôme d'ordre quelconque à travers eux. Je sais très bien effectuer cette régression linéaire. Cependant, disons que je veux `` forcer '' ma solution, pour parcourir, disons, trois de mes points de données aux coordonnées , x = 19 et x = 89 (et leurs coordonnées y correspondantes bien sûr).X=3X=19X=89

Comment appelle-t-on cette procédure générale et comment dois-je la connaître?

Éditer:

Je voudrais ajouter que je cherche un moyen concret de le faire. J'ai écrit un programme qui effectue en fait la régression linéaire de deux manières, en inversant la matrice de covariance directement ou par descente de gradient. Ce que je demande, c'est comment, exactement, étape par étape, je modifie ce que j'ai fait, de sorte que je force la solution polynomiale à passer par des points spécifiques?

Merci!

Spacey
la source
Pourquoi appelez-vous cela "linéaire" si vous utilisez un polynôme? Chaque point que vous souhaitez faire passer est une contrainte qui réduira votre degré de liberté. Vous pouvez ensuite utiliser un algorithme d'optimisation contraint.
curious_cat
4
Il est linéaire car vous trouvez des coefficients pour une combinaison linéaire . Par exemple, si vous souhaitez ajuster vos données à un cube, alors vous trouvez les coefficients (les ) de y = c 0 + c 1 x + c 2 x 2 + c 3 x 3 . cy=c0+c1X+c2X2+c3X3
Spacey
1
@Mohammad: Une autre façon d'approximer ce que vous voulez consisterait à utiliser une solution des moindres carrés pondérés et à attribuer des poids très importants aux points sur lesquels vous souhaitez que la ligne de régression passe. Cela devrait forcer la solution à passer très étroitement aux points que vous choisissez.
Jason R
@JasonR Ravi de vous voir ici. Oui WLS est en effet un concurrent intéressant. Je suis allé avec des réponses whubers en raison de la factorisation polynomiale intelligente, et parce qu'elle maintient bien la structure d'erreur.
Spacey

Réponses:

19

Le modèle en question peut être écrit

y=p(X)+(X-X1)(X-X)(β0+β1X++βpXp)+ε

est un polynôme de degré d - 1 passant par des points prédéterminés ( x 1 , y 1 ) , , ( x d , y d ) et ε est aléatoire. (Utilisez le polynôme interpolateur de Lagrange .) Écriture ( x - x 1 ) ( x - x d ) = rp(Xje)=yje-1(X1,y1),,(X,y)ε nous permet de réécrire ce modèle en(X-X1)(X-X)=r(X)

y-p(X)=β0r(X)+β1r(X)X+β2r(X)X2++βpr(X)Xp+ε,

qui est un problème de régression multiple OLS standard avec la même structure d'erreur que l'original où les variables indépendantes sont les quantités r ( x ) x i , i = 0 , 1 , , p . Calculez simplement ces variables et exécutez votre logiciel de régression familier , en veillant à l'empêcher d'inclure un terme constant. Les mises en garde habituelles concernant les régressions sans terme constant s'appliquent; en particulier, le R 2 peut être artificiellement élevé; les interprétations habituelles ne s'appliquent pas.p+1r(X)Xje, je=0,1,,pR2

(En fait, la régression par l'origine est un cas particulier de cette construction où , ( x 1 , y 1 ) = ( 0 , 0 ) , et p ( x ) = 0 , de sorte que le modèle est y = β 0 x + + β p x p + 1 + ε . )=1(X1,y1)=(0,0)p(X)=0y=β0X++βpXp+1+ε.


Voici un exemple travaillé (en 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)

Terrain

Les trois points fixes sont affichés en rouge continu - ils ne font pas partie des données. L'ajustement des moindres carrés polynomiaux du quatrième ordre non contraint est représenté par une ligne pointillée noire (il a cinq paramètres); l'ajustement contraint (de l'ordre de cinq, mais avec seulement trois paramètres libres) est indiqué par la ligne rouge.

Inspecter la sortie des moindres carrés ( summary(fit0)et summary(fit)) peut être instructif - je laisse cela au lecteur intéressé.

whuber
la source
βr(X)XjeXjer(X)
J'ai ajouté un exemple concret, Mohammad.
whuber
Oh, parfait. Je vais l'étudier. En utilisant votre exemple, il serait toujours possible de forcer le poly à passer par des points qui font partie des données, non?
Spacey
Absolument, cela peut être fait: mais soyez doublement prudent quant à l'interprétation des valeurs de p ou de toute autre statistique, car maintenant vos contraintes sont basées sur les données elles-mêmes.
whuber
Votre message m'a réveillé hier soir. Je me suis enseigné le LIP. (LIP est intéressant. C'est comme une décomposition de Fourier mais avec des polys).
Spacey
9

(Xje,yje)XjeXyjey

Si vous voulez forcer une ligne à passer par deux points dans un plan XY, c'est assez facile à faire aussi. Deux points quelconques peuvent correspondre à une ligne. Vous pouvez utiliser la formule point-pente pour calculer votre pente, puis utiliser l'un des points, la pente et l' équation d'une ligne pour trouver l'ordonnée à l'origine.

XX2


Je me sens obligé de mentionner à ce stade, cependant, que ce n'est peut-être pas une bonne chose à faire (à moins que votre théorie ne fournisse des raisons très solides pour le faire). Vous pouvez également examiner la régression bayésienne , où vous pouvez permettre à votre modèle de trouver la meilleure combinaison des informations dans vos données et certaines informations antérieures (que vous pourriez utiliser pour biaiser fortement votre interception vers zéro, par exemple, sans tout à fait le forçant).

gung - Réintégrer Monica
la source
1
Xjeyje
2
Bien que le fait d' ajouter trois points supplémentaires et de les pondérer ( réponse de la Glen_b) puisse créer un tel ajustement, l'interprétation de n'importe quel résultat statistique serait problématique: certains ajustements seraient nécessaires.
whuber
6

Pour ajouter un peu d'informations supplémentaires à l'excellente couverture de @ gung du cas linéaire, dans le cas polynomial d'ordre supérieur, il existe plusieurs façons de le faire exactement ou approximativement (mais à peu près aussi précisément que nécessaire).

Notons tout d'abord que les degrés de liberté du polynôme (ou bien de toute fonction ajustée) doivent être au moins aussi importants que le nombre de points "connus". Si les degrés de liberté sont égaux, vous n'avez pas du tout besoin des données, car la courbe est complètement déterminée. S'il y a plus de points «connus», vous ne pouvez pas le résoudre (à moins qu'ils ne reposent tous sur exactement le même polynôme du degré spécifié, auquel cas tout sous-ensemble de taille appropriée suffira). À partir de là, je vais juste parler du moment où le polynôme a plus de df que les points connus (comme un cube - avec 4df - et trois points connus, de sorte que le cube n'est ni surdéterminé par les points connus ni complètement déterminé par eux) .

1) "la courbe doit passer par ce point" est une contrainte linéaire sur les paramètres, entraînant une estimation contrainte ou des moindres carrés contraints (bien que les deux termes puissent inclure autre chose que des contraintes linéaires, telles que des contraintes de positivité). Vous pouvez incorporer des contraintes linéaires en

  (a) refonte de la paramétrisation pour inclure implicitement chaque contrainte résultant en un modèle d'ordre inférieur.

  (b) utiliser des outils standard qui peuvent incorporer des contraintes linéaires sur les paramètres d'un ajustement par moindres carrés. (généralement via quelque chose comme la formule donnée au lien ci-dessus)

2) Une autre méthode consiste à effectuer une régression pondérée. Si vous donnez aux points connus un poids suffisamment important, vous pouvez obtenir essentiellement le même ajustement qu'en (1). Ceci est souvent facilement implémenté, peut être beaucoup plus rapide que le reparamétrage et peut être effectué dans des packages qui n'offrent pas d'ajustement contraint.

Toutes les mises en garde de @ gung s'appliquent

Glen_b -Reinstate Monica
la source
Glen_b, je n'avais pas envisagé de régression pondérée. Ce pourrait être la façon de procéder. Je l'ai mis dans ma liste de tâches. Je crois que je peux m'enseigner cela sans incident. En ce qui concerne (1), pourriez-vous bien vouloir développer cet aspect de la re-paramatérisation? Aussi, comment appelez-vous cela que j'essaie de faire, où je force le polynôme à passer par certains points? Une partie du problème est que je ne sais pas quoi chercher sur Google. Si je sais comment cela s'appelait, je pourrais peut-être étoffer ce que vous dites avec du matériel en ligne. Merci.
Spacey
Voir mes modifications ci-dessus, qui incluent des termes de recherche et un lien avec quelques détails supplémentaires.
Glen_b -Reinstate Monica
2
+1 La régression pondérée est une bonne idée. Un ajustement des statistiques de sortie, telles que des estimations de l'erreur RMS, peut être nécessaire.
whuber
s2FR2
Merci pour votre réponse Glen_b, bien que j'aie accepté @whuber, j'ai quand même beaucoup appris de la vôtre.
Spacey