Czy prognozy modelu losowego lasu mają przedział prognozy?

Jeśli uruchomię randomForestmodel, mogę następnie przewidywać na podstawie modelu. Czy istnieje sposób na uzyskanie przedziału prognoz dla każdej z prognoz, tak że wiem, jak „pewny” model ma odpowiedź. Jeśli jest to możliwe, czy jest to po prostu oparte na zmienności zmiennej zależnej dla całego modelu, czy będzie miał szersze i węższe przedziały w zależności od konkretnego drzewa decyzyjnego, które zastosowano dla danej prognozy?

Jest to częściowo odpowiedź na @Sashikanth Dareddy (ponieważ nie zmieści się w komentarzu), a częściowo odpowiedź na oryginalny post.

Pamiętaj, co to jest przedział prognozy, to przedział lub zestaw wartości, w których przewidujemy, że przyszłe obserwacje będą leżały. Zasadniczo przedział predykcji składa się z 2 głównych elementów, które określają jego szerokość, elementu reprezentującego niepewność dotyczącą przewidywanej średniej (lub innego parametru), który jest częścią przedziału ufności, oraz elementu reprezentującego zmienność poszczególnych obserwacji wokół tej średniej. Przedział ufności jest dość solidny ze względu na centralne twierdzenie graniczne, aw przypadku losowego lasu pomaga również ładowanie początkowe. Jednak przedział przewidywania jest całkowicie zależny od założeń dotyczących dystrybucji danych, biorąc pod uwagę, że zmienne predykcyjne, CLT i ładowanie początkowe nie mają wpływu na tę część.

Przedział prognozy powinien być szerszy, a odpowiadający mu przedział ufności byłby również szerszy. Inne rzeczy, które mogłyby wpłynąć na szerokość przedziału prognozy, to założenia o równej wariancji lub nie, musi to wynikać z wiedzy badacza, a nie z losowego modelu lasu.

Przedział prognozy nie ma sensu dla wyniku kategorycznego (można zrobić zestaw prognoz zamiast przedziału, ale przez większość czasu prawdopodobnie nie byłby bardzo pouczający).

Możemy zobaczyć niektóre problemy dotyczące przedziałów prognoz, symulując dane, w których znamy dokładną prawdę. Rozważ następujące dane:

set.seed(1)

x1 <- rep(0:1, each=500)
x2 <- rep(0:1, each=250, length=1000)

y <- 10 + 5*x1 + 10*x2 - 3*x1*x2 + rnorm(1000)

Te konkretne dane są zgodne z założeniami regresji liniowej i są dość proste w przypadku losowego dopasowania lasu. Wiemy z „prawdziwego” modelu, że gdy oba predyktory wynoszą 0, to średnia wynosi 10, wiemy również, że poszczególne punkty mają rozkład normalny ze standardowym odchyleniem wynoszącym 1. Oznacza to, że 95% przedział predykcji oparty na doskonałej wiedzy dla punkty te wynosiłyby od 8 do 12 (właściwie od 8,04 do 11,96, ale zaokrąglanie upraszcza). Każdy szacowany interwał przewidywania powinien być szerszy niż ten (brak doskonałej informacji zwiększa szerokość w celu kompensacji) i powinien obejmować ten zakres.

Spójrzmy na odstępy czasu od regresji:

fit1 <- lm(y ~ x1 * x2)

newdat <- expand.grid(x1=0:1, x2=0:1)

(pred.lm.ci <- predict(fit1, newdat, interval='confidence'))
#        fit       lwr      upr
# 1 10.02217  9.893664 10.15067
# 2 14.90927 14.780765 15.03778
# 3 20.02312 19.894613 20.15162
# 4 21.99885 21.870343 22.12735

(pred.lm.pi <- predict(fit1, newdat, interval='prediction'))
#        fit      lwr      upr
# 1 10.02217  7.98626 12.05808
# 2 14.90927 12.87336 16.94518
# 3 20.02312 17.98721 22.05903
# 4 21.99885 19.96294 24.03476

Widzimy, że istnieje pewna niepewność w szacowanych średnich (przedział ufności), co daje nam przedział predykcji, który jest szerszy (ale obejmuje) zakres od 8 do 12.

Spójrzmy teraz na przedział oparty na indywidualnych prognozach poszczególnych drzew (należy się spodziewać, że będą one szersze, ponieważ losowy las nie korzysta z założeń (które, jak wiemy, są prawdziwe dla tych danych), jakie ma regresja liniowa):

library(randomForest)
fit2 <- randomForest(y ~ x1 + x2, ntree=1001)

pred.rf <- predict(fit2, newdat, predict.all=TRUE)

pred.rf.int <- apply(pred.rf$individual, 1, function(x) {
  c(mean(x) + c(-1, 1) * sd(x), 
  quantile(x, c(0.025, 0.975)))
})

t(pred.rf.int)
#                           2.5%    97.5%
# 1  9.785533 13.88629  9.920507 15.28662
# 2 13.017484 17.22297 12.330821 18.65796
# 3 16.764298 21.40525 14.749296 21.09071
# 4 19.494116 22.33632 18.245580 22.09904

Przedziały są szersze niż przedziały przewidywania regresji, ale nie obejmują całego zakresu. Obejmują one prawdziwe wartości, a zatem mogą być uzasadnione jako przedziały ufności, ale przewidują tylko, gdzie jest średnia (wartość przewidywana), a nie dodany element dla rozkładu wokół tej średniej. W pierwszym przypadku, gdy x1 i x2 są równe 0, przedziały nie spadają poniżej 9,7, jest to bardzo różna od prawdziwego przedziału przewidywania, który spada do 8. Jeśli wygenerujemy nowe punkty danych, będzie ich kilka (znacznie więcej niż 5%), które są w przedziałach prawdziwych i regresyjnych, ale nie mieszczą się w przypadkowych przedziałach leśnych.

Aby wygenerować przedział predykcji, musisz przyjąć pewne mocne założenia dotyczące rozmieszczenia poszczególnych punktów wokół przewidywanych średnich, a następnie możesz wziąć prognozy z poszczególnych drzew (element przedziału ufności bootstrapped), a następnie wygenerować losową wartość z założonej dystrybucja z tym centrum. Kwantyle tych wygenerowanych elementów mogą tworzyć przedział predykcji (ale nadal bym go przetestował, być może trzeba powtórzyć proces jeszcze kilka razy i połączyć).

Oto przykład zrobienia tego przez dodanie normalnych (ponieważ wiemy, że w oryginalnych danych zastosowano normalne) odchylenia do prognoz ze standardowym odchyleniem opartym na szacunkowym MSE z tego drzewa:

pred.rf.int2 <- sapply(1:4, function(i) {
  tmp <- pred.rf$individual[i, ] + rnorm(1001, 0, sqrt(fit2$mse))
  quantile(tmp, c(0.025, 0.975))
})
t(pred.rf.int2)
#           2.5%    97.5%
# [1,]  7.351609 17.31065
# [2,] 10.386273 20.23700
# [3,] 13.004428 23.55154
# [4,] 16.344504 24.35970

Te przedziały zawierają te oparte na doskonałej wiedzy, więc wyglądaj rozsądnie. Będą one jednak w dużym stopniu zależeć od przyjętych założeń (założenia są tutaj ważne, ponieważ wykorzystaliśmy wiedzę o tym, jak dane zostały zasymulowane, mogą nie być tak prawidłowe w rzeczywistych przypadkach danych). Nadal powtarzałbym symulacje kilka razy dla danych, które wyglądają bardziej jak twoje rzeczywiste dane (ale symulowane, abyś poznał prawdę) kilka razy, zanim w pełni zaufam tej metodzie.

Zdaję sobie sprawę, że to stary post, ale przeprowadziłem kilka symulacji na ten temat i pomyślałem, że podzielę się swoimi odkryciami.

$[ \mu + \sigma, \mu - \sigma]$$[\mu+1.96*\sigma,\mu-1.96*\sigma]$

Dokonując tej zmiany w kodzie @GregSnow, otrzymujemy następujące wyniki

set.seed(1)
x1 <- rep( 0:1, each=500 )
x2 <- rep( 0:1, each=250, length=1000 )
y <- 10 + 5*x1 + 10*x2 - 3*x1*x2 + rnorm(1000)

library(randomForest)
fit2 <- randomForest(y~x1+x2)
pred.rf <- predict(fit2, newdat, predict.all=TRUE)
pred.rf.int <- t(apply( pred.rf$individual, 1, function(x){ 
  c( mean(x) + c(-1.96,1.96)*sd(x), quantile(x, c(0.025,0.975)) )}))

pred.rf.int
                          2.5%    97.5%
1  7.826896 16.05521  9.915482 15.31431
2 11.010662 19.35793 12.298995 18.64296
3 14.296697 23.61657 14.749248 21.11239
4 18.000229 23.73539 18.237448 22.10331

Teraz, porównując je z przedziałami generowanymi przez dodanie normalnego odchylenia do prognoz ze standardowym odchyleniem, ponieważ MSE, takie jak @GregSnow, sugeruje, że otrzymamy,

pred.rf.int2 <- sapply(1:4, function(i) {
   tmp <- pred.rf$individual[i,] + rnorm(1000, 0, sqrt(fit2$mse))
   quantile(tmp, c(0.025, 0.975))
   })
t(pred.rf.int2)
          2.5%    97.5%
[1,]  7.486895 17.21144
[2,] 10.551811 20.50633
[3,] 12.959318 23.46027
[4,] 16.444967 24.57601

Przerwa między tymi dwoma podejściami jest teraz bardzo bliska. Wykreślenie przedziału prognoz dla trzech podejść w odniesieniu do rozkładu błędów w tym przypadku wygląda następująco

• Czarne linie = przedziały prognoz z regresji liniowej,
• Czerwone linie = losowe odstępy między lasami obliczane na podstawie indywidualnych prognoz,
• Niebieskie linie = losowe odstępy między lasami obliczane przez dodanie normalnego odchylenia do prognoz

Teraz ponownie uruchommy symulację, ale tym razem zwiększając wariancję terminu błędu. Jeśli nasze obliczenia przedziałów prognoz są dobre, powinniśmy skończyć z szerszymi przedziałami niż powyżej.

set.seed(1)
x1 <- rep( 0:1, each=500 )
x2 <- rep( 0:1, each=250, length=1000 )
y <- 10 + 5*x1 + 10*x2 - 3*x1*x2 + rnorm(1000,mean=0,sd=5)

fit1 <- lm(y~x1+x2)
newdat <- expand.grid(x1=0:1,x2=0:1)
predict(fit1,newdata=newdat,interval = "prediction")
      fit       lwr      upr
1 10.75006  0.503170 20.99695
2 13.90714  3.660248 24.15403
3 19.47638  9.229490 29.72327
4 22.63346 12.386568 32.88035

set.seed(1)
fit2 <- randomForest(y~x1+x2,localImp=T)
pred.rf.int <- t(apply( pred.rf$individual, 1, function(x){ 
  c( mean(x) + c(-1.96,1.96)*sd(x), quantile(x, c(0.025,0.975)) )}))
pred.rf.int
                          2.5%    97.5%
1  7.889934 15.53642  9.564565 15.47893
2 10.616744 18.78837 11.965325 18.51922
3 15.024598 23.67563 14.724964 21.43195
4 17.967246 23.88760 17.858866 22.54337

pred.rf.int2 <- sapply(1:4, function(i) {
   tmp <- pred.rf$individual[i,] + rnorm(1000, 0, sqrt(fit2$mse))
   quantile(tmp, c(0.025, 0.975))
   })
t(pred.rf.int2)
         2.5%    97.5%
[1,] 1.291450 22.89231
[2,] 4.193414 25.93963
[3,] 7.428309 30.07291
[4,] 9.938158 31.63777

Wyjaśnia to teraz, że obliczanie przedziałów prognozowania za pomocą drugiego podejścia jest znacznie dokładniejsze i daje wyniki całkiem zbliżone do przedziału prognozowania z regresji liniowej.

$\mu_i$$MSE_i$$N(\mu_i,RMSE_i)$$N(\sum \mu_i/n, \sum RMSE_i/n)$

mean.rf <- pred.rf$aggregate
sd.rf <- mean(sqrt(fit2$mse))
pred.rf.int3 <- cbind(mean.rf - 1.96*sd.rf, mean.rf + 1.96*sd.rf)
pred.rf.int3
1  1.332711 22.09364
2  4.322090 25.08302
3  8.969650 29.73058
4 10.546957 31.30789

Są one bardzo dobrze dopasowane do przedziałów modelu liniowego, a także sugerowanego podejścia @GregSnow. Należy jednak pamiętać, że podstawowym założeniem wszystkich omawianych metod jest to, że błędy mają rozkład normalny.

Jeśli używasz R, możesz łatwo stworzyć przedziały predykcji dla prognoz losowej regresji lasów: po prostu skorzystaj z pakietu quantregForest(dostępnego w CRAN ) i przeczytaj artykuł N. Meinshausena na temat tego, jak kwantyle warunkowe można wywnioskować z lasów regresji kwantowej i jak one może być użyty do zbudowania przedziałów prognoz. Bardzo pouczające, nawet jeśli nie pracujesz z R!

Łatwo to rozwiązać za pomocą randomForest.

Najpierw pozwól mi poradzić sobie z zadaniem regresji (zakładając, że twój las ma 1000 drzew). W predictfunkcji masz opcję zwracania wyników z poszczególnych drzew. Oznacza to, że otrzymasz 1000 kolumn wyjściowych. Możemy wziąć średnią z 1000 kolumn dla każdego wiersza - jest to normalna moc wyjściowa, którą RF wytworzyłby w jakikolwiek sposób. Teraz, aby uzyskać interwał przewidywania, powiedzmy +/- 2 std. odchylenia, co musisz zrobić, to dla każdego wiersza od 1000 wartości obliczyć +/- 2 standardowe. odchylenia i uczyńcie je górnymi i dolnymi granicami prognozy.

Po drugie, w przypadku klasyfikacji pamiętaj, że każde drzewo generuje 1 lub 0 (domyślnie), a suma wszystkich 1000 drzew podzielonych przez 1000 daje prawdopodobieństwo klasy (w przypadku klasyfikacji binarnej). Aby ustawić przedział predykcji na prawdopodobieństwo, musisz zmodyfikować min. opcja nodesize (dokładna nazwa tej opcji znajduje się w dokumentacji randomForest) po ustawieniu jej wartości >> 1 poszczególne drzewa będą generować liczby od 1 do 0. Teraz możesz powtórzyć ten sam proces, jak opisano powyżej dla zadanie regresji.

Mam nadzieję, że to ma sens.

Próbowałem kilka opcji (to wszystko WIP):

1. W rzeczywistości zmienną zależną uczyniłem problem klasyfikacji z wynikami jako zakresami, zamiast pojedynczej wartości. Wyniki, które uzyskałem, były słabe w porównaniu do zwykłej wartości. Zrezygnowałem z tego podejścia.

2. Następnie przekonwertowałem go na wiele problemów z klasyfikacją, z których każdy miał dolną granicę dla zakresu (wynikiem tego, czy model przekroczyłby dolną granicę, czy nie), a następnie uruchomiłem wszystkie modele (~ 20), a następnie połączono wynik, aby uzyskać ostateczną odpowiedź jako zakres. Działa to lepiej niż 1 powyżej, ale nie jest tak dobre, jak tego potrzebuję. Nadal pracuję nad ulepszeniem tego podejścia.

Użyłem OOB i pomijając szacunki, aby zdecydować, jak dobre / złe są moje modele.

Problem konstruowania przedziałów prognoz dla losowych prognoz lasów został omówiony w następującym artykule:

Zhang, Haozhe, Joshua Zimmerman, Dan Nettleton i Daniel J. Nordman. „Losowe przedziały prognozowania lasu”. The American Statistician, 2019.

Pakiet R „rfinterval” to jego implementacja dostępna w CRAN.

Instalacja

Aby zainstalować pakiet R rfinterval :

#install.packages("devtools")
#devtools::install_github(repo="haozhestat/rfinterval")
install.packages("rfinterval")
library(rfinterval)
?rfinterval

Stosowanie

Szybki start:

train_data <- sim_data(n = 1000, p = 10)
test_data <- sim_data(n = 1000, p = 10)

output <- rfinterval(y~., train_data = train_data, test_data = test_data,
                     method = c("oob", "split-conformal", "quantreg"),
                     symmetry = TRUE,alpha = 0.1)

### print the marginal coverage of OOB prediction interval
mean(output$oob_interval$lo < test_data$y & output$oob_interval$up > test_data$y)

### print the marginal coverage of Split-conformal prediction interval
mean(output$sc_interval$lo < test_data$y & output$sc_interval$up > test_data$y)

### print the marginal coverage of Quantile regression forest prediction interval
mean(output$quantreg_interval$lo < test_data$y & output$quantreg_interval$up > test_data$y)

Przykład danych:

oob_interval <- rfinterval(pm2.5 ~ .,
                            train_data = BeijingPM25[1:1000, ],
                            test_data = BeijingPM25[1001:2000, ],
                            method = "oob",
                            symmetry = TRUE,
                            alpha = 0.1)
str(oob_interval)