Jak wygładzić dane i wymusić monotoniczność

14

Mam pewne dane, które chciałbym wygładzić, aby wygładzone punkty monotonicznie zmniejszały się. Moje dane gwałtownie spadają, a następnie zaczynają się wyrównywać. Oto przykład z użyciem 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()

wykres danych do wygładzenia

Jakiej dobrej techniki wygładzania mógłbym użyć? Byłoby również miło, gdybym mógł zmusić pierwszy wygładzony punkt do zbliżenia się do mojego obserwowanego punktu.

Ben
źródło
1
Zauważyłem, że twoje przykładowe wartości są liczbami całkowitymi. Czy twoje prawdziwe wartości się liczą? Gdyby były, to (chociaż nie jest to gwarancja monotoniczności, w przypadku danych takich jak i tak zwykle da je), coś takiego może być przydatne:plot(y~x,data=df); f=fitted( glm( y~ns(x,df=4), data=df,family=quasipoisson)); lines(df$x,f)
Glen_b
Podobne pytanie Q z odpowiedzią: stats.stackexchange.com/questions/206073/...
kjetil b halvorsen 31.01.19

Odpowiedzi:

18

Można to zrobić za pomocą splajnowanych splajnów z ograniczeniami monotoniczności za pośrednictwem funkcji mono.con()i pcls()w pakiecie mgcv . Jest trochę do zrobienia, ponieważ te funkcje nie są tak przyjazne dla użytkownika gam(), ale kroki pokazano poniżej, oparte głównie na przykładzie z ?pcls, zmodyfikowane, aby pasowały do ​​podanych przez Ciebie danych przykładowych:

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!

Teraz musimy wypełnić obiekt, który przechodzi do pcls()zawierającego szczegółowe informacje o ograniczonym modelu ograniczonym, który chcemy dopasować

## 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

Teraz możemy w końcu dokonać dopasowania

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

pzawiera wektor współczynników dla funkcji bazowych odpowiadających splajnowi. Aby zwizualizować dopasowany splajn, możemy przewidzieć na podstawie modelu w 100 lokalizacjach w zakresie x. Robimy 100 wartości, aby uzyskać ładną, gładką linię na wykresie.

## 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)))

Do generowania przewidywanych wartości używamy Predict.matrix()macierzy, która po pomnożeniu przez współczynniki pdaje przewidywane wartości z dopasowanego modelu:

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")

Daje to:

wprowadź opis zdjęcia tutaj

Od ciebie zależy, czy dane zostaną uporządkowane w celu drukowania w ggplot ...

Możesz wymusić ściślejsze dopasowanie (aby częściowo odpowiedzieć na pytanie o płynniejsze dopasowanie pierwszego punktu danych) poprzez zwiększenie wymiaru funkcji podstawowej x. Na przykład ustawienie krówności 8( k <- 8) i ponowne uruchomienie kodu powyżej otrzymujemy

wprowadź opis zdjęcia tutaj

Nie można kznacznie zwiększyć tych danych i trzeba uważać, aby nie dopasować; wszystko, co pcls()robisz, to rozwiązanie problemu najmniejszych kwadratów podlegających karze, biorąc pod uwagę ograniczenia i dostarczone funkcje podstawowe, nie wykonuje dla ciebie wyboru gładkości - nie, że wiem o ...)

Jeśli chcesz interpolacji, zobacz podstawową funkcję R, ?splinefunktóra ma splajny Hermite i splajny sześcienne z ograniczeniami monotoniczności. W takim przypadku nie można tego jednak użyć, ponieważ dane nie są ściśle monotoniczne.

Przywróć Monikę - G. Simpson
źródło
Dzięki. Jestem pewien, że twoje rozwiązanie jest odpowiednie, ale jest tak złożone i zaciemnione, że po prostu nie mogę go użyć. splinefunbyła moja pierwsza myśl, jak również (jestem interpolację), ale spline(x=df$x, y=df$y, n=nrow(df), method="monoH.FC")i spline(x=df$x, y=df$y, n=nrow(df), method="hyman")oba błędy raise
Ben
1
Jeśli tylko spróbujesz, jestem pewien, że możesz go użyć; Nie mam pojęcia, co się tu dzieje, ale wymyśliłem to i wskazałem miejsca, w których trzeba coś zmienić. Zakładając, że wiesz trochę R oczywiście . Większość szczegółów ma charakter implementacyjny, który można zignorować, jeśli wszystko, co chcesz, pasuje do monotonicznie ograniczonego splajnu. Czy chcesz, żebym opatrzył kod trochę adnotacjami, aby lepiej podkreślić, co robi każdy krok? Odwołanie ?mono.conzawiera dalsze szczegóły dotyczące metody.
Przywróć Monikę - G. Simpson
Co do tego, dlaczego splinefunpodnosi błąd; Właśnie sobie uświadomiłem, że możesz dopasować monotoniczny splajn, który interpoluje dane, które same w sobie nie są monotoniczne. Obserwacja w x = 6jest większa yniż obserwacje w x = 5. Musisz tylko zignorować tę część odpowiedzi :-)
Przywróć Monikę - G. Simpson
Rozumiem. I nie ma potrzeby - jestem dość doświadczonym użytkownikiem R. Po prostu lubię rozumieć matematykę tego, czego używam, i wydaje się, że pod tym rozwiązaniem dzieje się całkiem sporo. Jeszcze raz dziękuję za pomoc.
Ben
Dodałem kilka notatek, aby wyjaśnić, co każda rzecz jest lub robi; należy przede wszystkim zauważyć, że ograniczenia monotoniczności są nakładane przez określony zestaw ograniczeń nierówności, które mono.conpowracają do splajnu sześciennego. ?pclszawiera przykłady cienkich wypustów płytowych i modeli addytywnych, które są mniej przyjazne dla użytkownika niż powyższe, ale które mogą ujawnić nieco więcej matematyki, jeśli znasz matematykę dla tych rodzajów splajnu (sam nie jestem taki zaznajomiony).
Przywróć Monikę - G. Simpson
13

Niedawny pakiet oszustwa autorstwa Natalii Pya i oparty na pracy Pya & Wood (2015) „Modele z ograniczeniami kształtu” autorstwa Pya & Wood (2015) mogą znacznie ułatwić proces wspomniany w doskonałej odpowiedzi Gavina.

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

Istnieje wiele funkcji bs, których możesz użyć - powyżej użyłem mpd do „monotonicznego zmniejszania splajnu P”, ale ma także funkcje, które wymuszają wypukłość lub wklęsłość osobno lub równolegle z ograniczeniami monotonicznymi.

srepho
źródło