2016-04-21 3 views
5

Voici le code que je couraissmooth.spline(): modèle ajusté ne correspond pas à degré spécifié par l'utilisateur de la liberté

fun <- function(x) {1 + 3*sin(4*pi*x-pi)} 
set.seed(1) 
num.samples <- 1000 
x <- runif(num.samples) 
y <- fun(x) + rnorm(num.samples) * 1.5 
fit <- smooth.spline(x, y, all.knots=TRUE, df=3) 

Malgré df=3, quand je vérifié le modèle ajusté, la sortie était

Call: 
smooth.spline(x = x, y = y, df = 3, all.knots = TRUE) 
Smoothing Parameter spar= 1.499954 lambda= 0.002508571 (26 iterations) 
Equivalent Degrees of Freedom (Df): 9.86422 

Quelqu'un pourrait-il m'aider? Merci!

+0

Avez-vous envisagé la possibilité que les degrés de liberté que vous fournissez soient un objectif que l'algorithme tente d'optimiser (entre autres critères) et qui soit aussi proche que possible de l'algorithme? – joran

Répondre

4

Notez qu'à partir de R-3.4.0 (2017-04-21), smooth.spline peut accepter la spécification directe de λ par un argument nouvellement ajouté lambda. Mais il sera encore converti en l'interne spar pendant l'estimation. Donc, la réponse suivante n'est pas affectée.


paramètre Lissage λ/spar se trouve dans le centre de contrôle de lissé

Lissé est contrôlée par le paramètre de lissage λ. smooth.spline() utilise un paramètre de lissage interne spar plutôt que λ:

spar = s0 + 0.0601 * log(λ) 

Une telle transformation est nécessaire logarithme pour faire la minimisation sans contrainte, comme GCV/CV. L'utilisateur peut spécifier spar pour spécifier indirectement λ. Lorsque spar se développe linéairement, λ va croître de façon exponentielle. Ainsi, il est rarement nécessaire d'utiliser une grande valeur spar.

Le degré de liberté df, est également définie en termes de λ:

edf

X est la matrice de modèle avec base B-spline et S est la matrice de pénalité.

Vous pouvez avoir un contrôle sur leurs relations avec votre ensemble de données:

spar <- seq(1, 2.5, by = 0.1) 
a <- sapply(spar, function (spar_i) unlist(smooth.spline(x, y, all.knots=TRUE, spar = spar_i)[c("df","lambda")])) 

Let croquis df ~ spar, λ ~ spar et log(λ) ~ spar:

par(mfrow = c(1,3)) 
plot(spar, a[1, ], type = "b", main = "df ~ spar", 
    xlab = "spar", ylab = "df") 
plot(spar, a[2, ], type = "b", main = "lambda ~ spar", 
    xlab = "spar", ylab = "lambda") 
plot(spar, log(a[2,]), type = "b", main = "log(lambda) ~ spar", 
    xlab = "spar", ylab = "log(lambda)") 

plot

Notez la croissance radicale de λ avec spar, la relation linéaire entre log(λ) et spar, et la relation relativement lisse entre df et spar.


smooth.spline() itérations montage de pneu pour spar

Si nous spécifions manuellement la valeur de spar, comme ce que nous avons fait dans le sapply(), pas d'itérations d'ajustement est fait pour la sélection spar; sinon, smooth.spline() doit parcourir un certain nombre de valeurs spar. Si nous

  • spécifions cv = TRUE/FALSE, les itérations d'ajustement visent à minimiser le score CV/GCV;
  • spécifier df = mydf, ajustements itérations vise à minimiser (df(spar) - mydf)^2.

La minimisation du GCV est facile à suivre. Nous ne nous soucions pas du score GCV, mais attention à la spar correspondante. Au contraire, en minimisant (df(spar) - mydf)^2, on se préoccupe souvent de la valeur df à la fin de l'itération plutôt que spar! Mais en gardant à l'esprit qu'il s'agit d'un problème de minimisation, nous ne sommes jamais garanti que le df final correspond à notre valeur cible mydf.


Pourquoi vous mettre df = 3, mais obtenir df = 9.864?

La fin de l'itération, peut soit frapper implique un minimum, ou d'atteindre la limite de recherche, ou d'atteindre nombre maximum d'itérations.

Nous sommes loin de la limite d'itérations maximale (par défaut 500); Pourtant, nous ne frappons pas le minimum. Eh bien, nous pourrions atteindre la limite. Ne pas se concentrer sur df, pensez à spar.

smooth.spline(x, y, all.knots=TRUE, df=3)$spar # 1.4999 

Selon ?smooth.spline, par défaut, smooth.spline() recherches spar entre [-1.5, 1.5]. C'est-à-dire, lorsque vous mettez df = 3, la minimisation se termine à la limite de la recherche, plutôt que de frapper df = 3. Jetez un oeil à notre graphique de la relation entre df et spar, encore une fois. De la figure, il semble que nous ayons besoin d'une valeur de spar près de 2 pour aboutir à df = 3.

Let utiliser l'argument control.spar:

fit <- smooth.spline(x, y, all.knots=TRUE, df=3, control.spar = list(high = 2.5)) 
# Smoothing Parameter spar= 1.859066 lambda= 0.9855336 (14 iterations) 
# Equivalent Degrees of Freedom (Df): 3.000305 

Maintenant, vous voyez, vous vous retrouvez avec df = 3. Et nous avons besoin d'un spar = 1.86.


Une meilleure suggestion: Ne pas utiliser all.knots = TRUE

Regardez, vous avez 1000 données. Avec all.knots = TRUE, vous utiliserez 1000 paramètres.Souhaitant finir avec df = 3 implique que 997 sur 1000 paramètres sont supprimés. Imaginez la taille d'un λ d'où spar dont vous avez besoin!

Essayez plutôt d'utiliser la spline de régression pénalisée. 200 paramètres à supprimer 3 est sans aucun doute beaucoup plus facile:

fit <- smooth.spline(x, y, nknots = 200, df=3) ## using 200 knots 
# Smoothing Parameter spar= 1.317883 lambda= 0.9853648 (16 iterations) 
# Equivalent Degrees of Freedom (Df): 3.000386 

Maintenant, vous vous retrouvez avec df = 3 sans spar contrôle.