Metody statystyczne do bardziej wydajnego drukowania danych, gdy istnieją miliony punktów?

Uważam, że generowanie wykresów przez R może zająć dużo czasu, gdy obecne są miliony punktów - nic dziwnego, biorąc pod uwagę, że punkty są drukowane indywidualnie. Ponadto takie wykresy są często zbyt zagracone i gęste, aby były przydatne. Wiele punktów nakłada się i tworzy czarną masę, a wiele czasu spędza na wykreślaniu większej liczby punktów w tę masę.

Czy istnieją jakieś statystyczne alternatywy dla reprezentowania dużych danych w standardowym wykresie rozrzutu? Rozważałem wykres gęstości, ale jakie są inne alternatywy?$n$

Jest to trudne zadanie bez gotowych rozwiązań (dzieje się tak oczywiście, ponieważ wykres gęstości jest tak kuszącym powrotem, jak nikt tak naprawdę nie obchodzi). Więc co możesz zrobić?

Jeśli naprawdę się pokrywają (tzn. Mają dokładnie takie same współrzędne X i Y), a nie używasz alfa, najlepszym pomysłem byłoby po prostu zmniejszenie nakładania się przy użyciu unique(w przypadku alfa można to zsumować w takich grupach).

Jeśli nie, możesz ręcznie zaokrąglić współrzędne do najbliższych pikseli i użyć poprzedniej metody (ale jest to brudne rozwiązanie).

Na koniec możesz wykonać wykres gęstości, aby użyć go do podpróbkowania punktów w najbardziej gęstych obszarach. To z drugiej strony nie stworzy dokładnie tej samej fabuły i może wprowadzić artefakty, jeśli nie zostanie dokładnie dostrojone.

Spójrz na pakiet hexbin , który implementuje papier / metodę Dana Carra. Pdf winieta ma więcej szczegółów, które cytuję poniżej:

1. Przegląd

Sortowanie sześciokątne jest formą dwuwymiarowego histogramu przydatnego do wizualizacji struktury w zestawach danych o dużym n. Podstawowa koncepcja binoksu sześciokątnego jest niezwykle prosta;

1. płaszczyzna xy nad zestawem (zakres (x), zakres (y)) jest mozaikowana regularną siatką sześciokątów.
2. liczba punktów przypadających na każdy sześciokąt jest zliczana i zapisywana w strukturze danych
3. sześciokąty o liczbie> 0 są wykreślane przy użyciu kolorowej rampy lub zmieniając promień sześciokąta proporcjonalnie do liczby. Podstawowy algorytm jest niezwykle szybki i skuteczny do wyświetlania struktury zbiorów danych za pomocą$n \ge 10^6$

Jeśli rozmiar siatki i cięcia w kolorowej rampie są wybierane w sprytny sposób, struktura nieodłącznie związana z danymi powinna pojawić się na spiętych wykresach. Te same zastrzeżenia dotyczą łączenia sześciokątnego, jak w przypadku histogramów i należy zachować ostrożność przy wyborze parametrów łączenia

Muszę przyznać, że nie do końca rozumiem twój ostatni akapit:

„Nie szukam wykresu gęstości (chociaż są one często przydatne), chciałbym uzyskać taki sam wynik jak proste wywołanie wykresu, ale w miarę możliwości znacznie szybciej niż miliony wykresów.”

Nie jest również jasne, jakiego rodzaju wykresu (funkcji) szukasz.

Biorąc pod uwagę, że masz zmienne metryczne, przydatne mogą się okazać sześciokąty lub pola słonecznika. W celu uzyskania dalszych informacji patrz

• Grafika dużych zestawów danych autorstwa Unwin / Theus / Hofmann
• Quick-R na „ Wykresach rozproszenia o wysokiej gęstości ”
• stat_hexbin w ggplot2

Inną bezpośrednią odpowiedzią na to pytanie jest pakiet rgl, który może wykreślić miliony punktów za pomocą OpenGL. Określ także rozmiar punktu (np. 3) i pomniejsz, aby zobaczyć te centra mas jako monolityczne bloki, lub powiększ i zobacz strukturę tego, co kiedyś było monolityczne - rozmiary punktów są stałe, ale odległości między nimi na ekranie zależy od powiększenia. Można także użyć poziomów alfa.

Oto plik, do którego dzwonię bigplotfix.R. Jeśli je podasz, zdefiniuje opakowanie, dla plot.xyktórego „kompresuje” dane wykresu, gdy jest ono bardzo duże. Opakowanie nie robi nic, jeśli dane wejściowe są małe, ale jeśli dane wejściowe są duże, to dzieli je na części i po prostu wykreśla maksymalną i minimalną wartość xiy dla każdej części. Sourcing ponownie bigplotfix.Rwiąże, graphics::plot.xyaby wskazać opakowanie (wielokrotne pozyskiwanie jest OK).

Należy pamiętać, że plot.xyjest to funkcja „koń pociągowy” dla standardowych metod kreślenia takich jak plot(), lines()i points(). Dzięki temu możesz nadal używać tych funkcji w kodzie bez modyfikacji, a duże wykresy będą automatycznie kompresowane.

To jest przykładowy wynik. Zasadniczo jest to plot(runif(1e5))z punktami i liniami oraz z wprowadzoną tutaj kompresją i bez niej. Wykres „skompresowanych punktów” nie trafia w środkowy obszar ze względu na charakter kompresji, ale wykres „linii skompresowanych” wygląda znacznie bliżej nieskompresowanego oryginału. Czasy dotyczą png()urządzenia; z jakiegoś powodu punkty są znacznie szybsze w pngurządzeniu niż w X11urządzeniu, ale przyspieszenia w X11są porównywalne ( X11(type="cairo")było wolniejsze niż X11(type="Xlib")w moich eksperymentach).

Powodem, dla którego to napisałem, jest to, że zmęczyło plot()mnie przypadkowe uruchamianie dużego zbioru danych (np. Pliku WAV). W takich przypadkach musiałbym wybierać między odczekaniem kilku minut na zakończenie wykresu a zakończeniem sesji R sygnałem (w ten sposób tracąc moją historię poleceń i zmienne). Teraz, jeśli pamiętam, aby załadować ten plik przed każdą sesją, w takich przypadkach mogę uzyskać użyteczny wykres. Mały komunikat ostrzegawczy wskazuje, kiedy dane wydruku zostały „skompresowane”.

# bigplotfix.R
# 28 Nov 2016

# This file defines a wrapper for plot.xy which checks if the input
# data is longer than a certain maximum limit. If it is, it is
# downsampled before plotting. For 3 million input points, I got
# speed-ups of 10-100x. Note that if you want the output to look the
# same as the "uncompressed" version, you should be drawing lines,
# because the compression involves taking maximum and minimum values
# of blocks of points (try running test_bigplotfix() for a visual
# explanation). Also, no sorting is done on the input points, so
# things could get weird if they are out of order.
test_bigplotfix = function() {
  oldpar=par();
  par(mfrow=c(2,2))
  n=1e5;
  r=runif(n)
  bigplotfix_verbose<<-T
  mytitle=function(t,m) { title(main=sprintf("%s; elapsed=%0.4f s",m,t["elapsed"])) }
  mytime=function(m,e) { t=system.time(e); mytitle(t,m); }

  oldbigplotfix_maxlen = bigplotfix_maxlen
  bigplotfix_maxlen <<- 1e3;

  mytime("Compressed, points",plot(r));
  mytime("Compressed, lines",plot(r,type="l"));
  bigplotfix_maxlen <<- n
  mytime("Uncompressed, points",plot(r));
  mytime("Uncompressed, lines",plot(r,type="l"));
  par(oldpar);
  bigplotfix_maxlen <<- oldbigplotfix_maxlen
  bigplotfix_verbose <<- F
}

bigplotfix_verbose=F

downsample_xy = function(xy, n, xlog=F) {
  msg=if(bigplotfix_verbose) { message } else { function(...) { NULL } }
  msg("Finding range");
  r=range(xy$x);
  msg("Finding breaks");
  if(xlog) {
    breaks=exp(seq(from=log(r[1]),to=log(r[2]),length.out=n))
  } else {
    breaks=seq(from=r[1],to=r[2],length.out=n)
  }
  msg("Calling findInterval");
  ## cuts=cut(xy$x,breaks);
  # findInterval is much faster than cuts!
  cuts = findInterval(xy$x,breaks);
  if(0) {
    msg("In aggregate 1");
    dmax = aggregate(list(x=xy$x, y=xy$y), by=list(cuts=cuts), max)
    dmax$cuts = NULL;
    msg("In aggregate 2");
    dmin = aggregate(list(x=xy$x, y=xy$y), by=list(cuts=cuts), min)
    dmin$cuts = NULL;
  } else { # use data.table for MUCH faster aggregates
    # (see http://stackoverflow.com/questions/7722493/how-does-one-aggregate-and-summarize-data-quickly)
    suppressMessages(library(data.table))
    msg("In data.table");
    dt = data.table(x=xy$x,y=xy$y,cuts=cuts)
    msg("In data.table aggregate 1");
    dmax = dt[,list(x=max(x),y=max(y)),keyby="cuts"]
    dmax$cuts=NULL;
    msg("In data.table aggregate 2");
    dmin = dt[,list(x=min(x),y=min(y)),keyby="cuts"]
    dmin$cuts=NULL;
    #  ans = data_t[,list(A = sum(count), B = mean(count)), by = 'PID,Time,Site']
  }
  msg("In rep, rbind");
  # interleave rows (copied from a SO answer)
  s <- rep(1:n, each = 2) + (0:1) * n
  xy = rbind(dmin,dmax)[s,];
  xy
}

library(graphics);
# make sure we don't create infinite recursion if someone sources
# this file twice
if(!exists("old_plot.xy")) {
  old_plot.xy = graphics::plot.xy
}

bigplotfix_maxlen = 1e4

# formals copied from graphics::plot.xy
my_plot.xy = function(xy, type, pch = par("pch"), lty = par("lty"),
  col = par("col"), bg = NA, cex = 1, lwd = par("lwd"),
  ...) {

  if(bigplotfix_verbose) {
    message("In bigplotfix's plot.xy\n");
  }

  mycall=match.call();
  len=length(xy$x)
  if(len>bigplotfix_maxlen) {
    warning("bigplotfix.R (plot.xy): too many points (",len,"), compressing to ",bigplotfix_maxlen,"\n");
    xy = downsample_xy(xy, bigplotfix_maxlen, xlog=par("xlog"));
    mycall$xy=xy
  }
  mycall[[1]]=as.symbol("old_plot.xy");

  eval(mycall,envir=parent.frame());
}

# new binding solution adapted from Henrik Bengtsson
# https://stat.ethz.ch/pipermail/r-help/2008-August/171217.html
rebindPackageVar = function(pkg, name, new) {
  # assignInNamespace() no longer works here, thanks nannies
  ns=asNamespace(pkg)
  unlockBinding(name,ns)
  assign(name,new,envir=asNamespace(pkg),inherits=F)
  assign(name,new,envir=globalenv())
  lockBinding(name,ns)
}
rebindPackageVar("graphics", "plot.xy", my_plot.xy);

Może zostanę odrzucony za moją metodę, złe wspomnienia jednego z moich profesorów badawczych krzyczących na ludzi za wyrzucanie dobrych danych przez tłumaczenie ich na kategorie (oczywiście, zgadzam się teraz kilka dni lol), nie wiem. W każdym razie, jeśli mówisz o wykresie rozrzutu, to miałem te same problemy. Teraz, gdy mam dane liczbowe, nie ma dla mnie sensu klasyfikowanie ich do analizy. Ale wizualizacja to inna historia. To, co według mnie działa najlepiej, to najpierw (1) rozbić zmienną niezależną na grupy za pomocą funkcji cięcia. Możesz grać z liczbą grup, a następnie (2) po prostu wykreślić DV względem wyciętej wersji IV. R wygeneruje wykresy pudełkowe zamiast tego obrzydliwego wykresu punktowego. Polecam usunięcie wartości odstających z wykresu (użyj opcji konspekt = FAŁSZ w poleceniu plot). Ponownie, NIGDY nie zmarnowałbym idealnie dobrych danych liczbowych, kategoryzując, a następnie analizując. Robi to zbyt wiele problemów. Chociaż wiem, że to drażliwy temat debaty. Ale robienie tego specjalnie w celu przynajmniej wizualnego zrozumienia danych, nie wyrządzając przy tym wiele szkody. Wykreśliłem dane tak duże, jak 10M i nadal udało mi się to zrozumieć dzięki tej metodzie. Mam nadzieję, że to pomaga! Z poważaniem! widzieliśmy z tego. Wykreśliłem dane tak duże, jak 10M i nadal udało mi się to zrozumieć dzięki tej metodzie. Mam nadzieję, że to pomaga! Z poważaniem! widzieliśmy z tego. Wykreśliłem dane tak duże, jak 10M i nadal udało mi się to zrozumieć dzięki tej metodzie. Mam nadzieję, że to pomaga! Z poważaniem!

W przypadku dużych serii czasowych pokochałem smoothScatter (nie mniej niż część podstawy R). Często muszę podawać dodatkowe dane, a zachowanie podstawowego interfejsu API fabuły jest naprawdę pomocne, na przykład:

set.seed(1)
ra <- rnorm(n = 100000, sd = 1, mean = 0)
smoothScatter(ra)
abline(v=25000, col=2)
text(25000, 0, "Event 1", col=2)

Co daje (jeśli wybaczysz projekt):

Jest zawsze dostępny i działa dobrze z ogromnymi zestawami danych, więc miło jest przynajmniej spojrzeć na to, co masz.