Как сгладить данные и заставить монотонность

14

У меня есть некоторые данные, которые я хотел бы сгладить так, чтобы сглаженные точки монотонно уменьшались. Мои данные резко уменьшаются и затем начинают плато. Вот пример использования R

df <- data.frame(x=1:10, y=c(100,41,22,10,6,7,2,1,3,1))
ggplot(df, aes(x=x, y=y))+geom_line()

график данных для сглаживания

Какую хорошую технику сглаживания я мог бы использовать? Кроме того, было бы хорошо, если бы я мог заставить первую сглаженную точку быть близкой к моей наблюдаемой точке.

Бен
источник
1
Я заметил, что ваш пример значения являются целыми числами. Ваши реальные ценности имеют значение? Если бы они были, тогда (хотя это и не является гарантией монотонности, для данных, подобных этим, они обычно все равно будут их предоставлять), может пригодиться что-то вроде этого:plot(y~x,data=df); f=fitted( glm( y~ns(x,df=4), data=df,family=quasipoisson)); lines(df$x,f)
Glen_b
Аналогичный вопрос с ответом: stats.stackexchange.com/questions/206073/…
kjetil b halvorsen

Ответы:

18

Вы можете сделать это с помощью наказываться шлицы с ограничениями монотонности через mono.con()и pcls()функций в mgcv пакете. Нам нужно немного поработать, потому что эти функции не так удобны для пользователя, как gam()показано ниже, но шаги показаны ниже, в основном на примере из ?pcls, модифицированного в соответствии с данными, которые вы дали:

df <- data.frame(x=1:10, y=c(100,41,22,10,6,7,2,1,3,1))

## Set up the size of the basis functions/number of knots
k <- 5
## This fits the unconstrained model but gets us smoothness parameters that
## that we will need later
unc <- gam(y ~ s(x, k = k, bs = "cr"), data = df)

## This creates the cubic spline basis functions of `x`
## It returns an object containing the penalty matrix for the spline
## among other things; see ?smooth.construct for description of each
## element in the returned object
sm <- smoothCon(s(x, k = k, bs = "cr"), df, knots = NULL)[[1]]

## This gets the constraint matrix and constraint vector that imposes
## linear constraints to enforce montonicity on a cubic regression spline
## the key thing you need to change is `up`.
## `up = TRUE` == increasing function
## `up = FALSE` == decreasing function (as per your example)
## `xp` is a vector of knot locations that we get back from smoothCon
F <- mono.con(sm$xp, up = FALSE)   # get constraints: up = FALSE == Decreasing constraint!

Теперь нам нужно заполнить объект, который передается, pcls()содержащим детали штрафной модели с ограничениями, которую мы хотим подогнать.

## Fill in G, the object pcsl needs to fit; this is just what `pcls` says it needs:
## X is the model matrix (of the basis functions)
## C is the identifiability constraints - no constraints needed here
##   for the single smooth
## sp are the smoothness parameters from the unconstrained GAM
## p/xp are the knot locations again, but negated for a decreasing function
## y is the response data
## w are weights and this is fancy code for a vector of 1s of length(y)
G <- list(X = sm$X, C = matrix(0,0,0), sp = unc$sp,
          p = -sm$xp, # note the - here! This is for decreasing fits!
      y = df$y,
          w = df$y*0+1)
G$Ain <- F$A    # the monotonicity constraint matrix
G$bin <- F$b    # the monotonicity constraint vector, both from mono.con
G$S <- sm$S     # the penalty matrix for the cubic spline
G$off <- 0      # location of offsets in the penalty matrix

Теперь мы можем наконец сделать примерку

## Do the constrained fit 
p <- pcls(G)  # fit spline (using s.p. from unconstrained fit)

pсодержит вектор коэффициентов для базисных функций, соответствующих сплайну. Чтобы визуализировать подобранный сплайн, мы можем предсказать по модели в 100 точках в диапазоне х. Мы делаем 100 значений, чтобы получить красивую плавную линию на графике.

## predict at 100 locations over range of x - get a smooth line on the plot
newx <- with(df, data.frame(x = seq(min(x), max(x), length = 100)))

Для генерации прогнозируемых значений мы используем Predict.matrix(), который генерирует матрицу, которая при умножении на коэффициенты pдает прогнозируемые значения из подобранной модели:

fv <- Predict.matrix(sm, newx) %*% p
newx <- transform(newx, yhat = fv[,1])

plot(y ~ x, data = df, pch = 16)
lines(yhat ~ x, data = newx, col = "red")

Это производит:

введите описание изображения здесь

Я оставлю это на ваше усмотрение, чтобы получить данные в удобной форме для построения графиков с помощью ggplot ...

Вы можете принудительно подгонять (частично ответить на ваш вопрос о сглаживании подгонки к первой точке данных), увеличив размер базовой функции x. Например, установив kравным 8( k <- 8) и перезапустив код выше, мы получим

введите описание изображения здесь

Вы не можете продвинуться kнамного выше для этих данных, и вы должны быть осторожны с перебором; все, что pcls()мы делаем, это решаем проблему штрафованных наименьших квадратов, учитывая ограничения и предоставляемые базовые функции, он не выполняет выбор гладкости для вас - не то, что я знаю ...)

Если вы хотите интерполяцию, то посмотрите базовую функцию R, ?splinefunкоторая имеет сплайны Эрмита и кубические сплайны с ограничениями монотонности. В этом случае вы не можете использовать это, однако, поскольку данные не являются строго монотонными.

Восстановить Монику - Дж. Симпсон
источник
Благодарю. Я уверен, что ваше решение подходящее, но оно настолько сложное и запутанное, что я просто не могу его использовать. splinefunбыла моя первоначальная мысль (я интерполирую), но spline(x=df$x, y=df$y, n=nrow(df), method="monoH.FC")и spline(x=df$x, y=df$y, n=nrow(df), method="hyman")обе ошибки вызывают
Бен
1
Если вы просто попробуете, я уверен, что вы можете использовать это; Я мало представляю, что здесь происходит под капотом, но я решил, и я указал места, где вам нужно что-то изменить. Предполагая, что вы знаете немного R, конечно . Большая часть деталей - это реализация, которую вы можете игнорировать, если все, что вы хотите сделать, соответствует монотонно ограниченному сплайну. Хотели бы вы, чтобы я немного больше комментировал код, чтобы больше рассказать о том, что делает каждый шаг? Ссылка ?mono.conсодержит более подробную информацию о методе.
Восстановить Монику - Г. Симпсон
Что касается того, почему splinefunвозникает ошибка; Я только что понял, что вы можете использовать монотонный сплайн, который интерполирует данные, которые сами по себе не являются монотонными. Наблюдение в x = 6больше y, чем в x = 5. Вы просто должны проигнорировать эту часть ответа :-)
Восстановить Монику - Дж. Симпсон
Понял. И не нужно - я довольно опытный пользователь R. Мне просто нравится понимать математику того, что я использую, и у этого решения, похоже, много чего происходит под капотом. Еще раз спасибо за помощь.
Бен
Я добавил несколько заметок, чтобы объяснить, что каждая вещь делает или делает; главное, на что следует обратить внимание: ограничения монотонности накладываются определенным набором ограничений неравенства, которые mono.conвозвращаются для кубического сплайна. ?pclsесть примеры сплайнов на тонких пластинах и аддитивных моделей, которые менее удобны для пользователя, чем приведенные выше, но которые могут раскрыть некоторые математические знания, если вы знакомы с математикой для этих типов сплайнов (я сам не настолько знаком).
Восстановить Монику - Г. Симпсон
13

Недавний мошеннический пакет Натальи Пи, основанный на статье «Модели с аддитивными формами» от Pya & Wood (2015), может значительно облегчить процесс, упомянутый в превосходном ответе Гевина.

library(scam)
con <- scam(y ~ s(x, k = k, bs = "mpd"), data = df)
plot(con)

Есть ряд функций bs, которые вы можете использовать - в приведенном выше примере я использовал mpd для «монотонного убывающего P-сплайна», но он также имеет функции, которые обеспечивают выпуклость или вогнутость либо отдельно, либо вместе с монотонными ограничениями.

srepho
источник