Średnia geometryczna: czy jest wbudowana?

106

Próbowałem znaleźć wbudowaną średnią geometryczną, ale nie mogłem.

(Oczywiście wbudowany nie zaoszczędzi mi czasu podczas pracy w powłoce, ani nie podejrzewam żadnej różnicy w dokładności; w przypadku skryptów staram się używać wbudowanych tak często, jak to możliwe, gdzie (skumulowany) wzrost wydajności jest często zauważalny.

Na wypadek, gdyby nie było żadnego (co wątpię), oto moje.

gm_mean = function(a){prod(a)^(1/length(a))}
Doug
źródło
11
Uważaj na liczby ujemne i przepełnienia. prod (a) bardzo szybko spadnie lub przepełni. Próbowałem zmierzyć czas, używając dużej listy i szybko uzyskałem Inf, używając twojej metody vs 1.4 z exp (mean (log (x))); problem zaokrągleń może być dość poważny.
Tristan
po prostu szybko napisałem powyższą funkcję, ponieważ byłem pewien, że 5 minut po wysłaniu tego Q ktoś powie mi, że R jest wbudowany w gm. Więc nie jest wbudowany, więc na pewno warto poświęcić czas na ponowne kodowanie w świetle twoich uwag. +1 ode mnie.
Doug
1
Właśnie oznaczyłem tę geometryczną średnią i wbudowaną , 9 lat później.
smci

Odpowiedzi:

79

Oto wektoryzowana funkcja z tolerancją na zero i NA, służąca do obliczania średniej geometrycznej w R. Pełne meanobliczenia z udziałem length(x)są konieczne w przypadkach, w których xzawiera wartości niedodatnie.

gm_mean = function(x, na.rm=TRUE){
  exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}

Podziękowania dla @ ben-bolker za odnotowanie na.rmprzejścia i @Gregor za upewnienie się, że działa poprawnie.

Myślę, że niektóre komentarze dotyczą fałszywej równoważności NAwartości w danych i zer. W aplikacji, o której myślałem, są takie same, ale oczywiście nie jest to generalnie prawda. Tak więc, jeśli chcesz uwzględnić opcjonalną propagację zer i traktować length(x)inaczej w przypadku NAusuwania, poniższa alternatywa jest nieco dłuższą alternatywą dla powyższej funkcji.

gm_mean = function(x, na.rm=TRUE, zero.propagate = FALSE){
  if(any(x < 0, na.rm = TRUE)){
    return(NaN)
  }
  if(zero.propagate){
    if(any(x == 0, na.rm = TRUE)){
      return(0)
    }
    exp(mean(log(x), na.rm = na.rm))
  } else {
    exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
  }
}

Należy zauważyć, że sprawdza również wszelkie wartości ujemne i zwraca bardziej pouczające i odpowiednie informacje, biorąc pod uwagę, NaNże średnia geometryczna nie jest zdefiniowana dla wartości ujemnych (ale jest dla zer). Dziękuję komentatorom, którzy pozostali przy mojej sprawie w tej sprawie.

Paul McMurdie
źródło
2
czy nie byłoby lepiej przejść na.rmjako argument (tj. pozwolić użytkownikowi zdecydować, czy chce być tolerancyjny na NA, czy nie, w celu zachowania spójności z innymi funkcjami podsumowującymi R)? Boję się automatycznego wykluczania zer - również bym to zrobił.
Ben Bolker
1
Być może masz rację, jeśli chodzi o podanie na.rmopcji. Zaktualizuję odpowiedź. Jeśli chodzi o wykluczanie zer, średnia geometryczna jest nieokreślona dla wartości niedodatnich, w tym zer. Powyższe jest częstą poprawką dla średniej geometrycznej, w której zerom (lub w tym przypadku wszystkim różnym od zera) przypisywana jest wartość fikcyjna 1, co nie ma wpływu na iloczyn (lub równoważnie zero w sumie logarytmicznej).
Paul McMurdie
* Miałem na myśli wspólną poprawkę dla wartości niedodatnich, przy czym zero jest najczęściej używane, gdy używana jest średnia geometryczna.
Paul McMurdie
1
Twoje na.rmprzejście nie działa zgodnie z kodem ... zobacz gm_mean(c(1:3, NA), na.rm = T). Musisz usunąć & !is.na(x)z podzbioru wektorów, a ponieważ pierwszy argument sumto ..., musisz przekazać na.rm = na.rmnazwę, a także musisz wykluczyć 0„i NA” z wektora w lengthwywołaniu.
Gregor Thomas
2
Uwaga: ponieważ xzawiera tylko zero (s), na przykład x <- 0, exp(sum(log(x[x>0]), na.rm = TRUE)/length(x))podaje 1średnią geometryczną, która nie ma sensu.
adatum
88

Nie, ale jest kilka osób, które takie napisały, na przykład tutaj .

Inną możliwością jest użycie tego:

exp(mean(log(x)))
Mark Byers
źródło
Kolejną zaletą używania exp (mean (log (x))) jest to, że możesz pracować z długimi listami dużych liczb, co jest problematyczne, gdy używasz bardziej oczywistego wzoru używającego prod (). Zauważ, że prod (a) ^ (1 / length (a)) i exp (mean (log (a))) dają tę samą odpowiedź.
lukeholman
link został naprawiony
PatrickT
15

Możemy użyć pakietu psych i wywołać funkcję geometryczną .

AliCivil
źródło
1
psych::geometric.mean()
smci
Te funkcje powinny przyjąć serie, a nie ich wzrost, przynajmniej jako opcja, powiedziałbym.
Christoph Hanck,
12

Plik

exp(mean(log(x)))

będzie działać, chyba że jest 0 w x. Jeśli tak, dziennik wygeneruje -Inf (-Infinite), co zawsze daje średnią geometryczną równą 0.

Jednym z rozwiązań jest usunięcie wartości -Inf przed obliczeniem średniej:

geo_mean <- function(data) {
    log_data <- log(data)
    gm <- exp(mean(log_data[is.finite(log_data)]))
    return(gm)
}

Można do tego użyć linijki jednowierszowej, ale oznacza to dwukrotne obliczenie dziennika, co jest nieefektywne.

exp(mean(log(i[is.finite(log(i))])))
Alan James Salmoni
źródło
po co obliczać dziennik dwa razy, skoro możesz: exp (mean (x [x! = 0])))
zzk
1
w obu podejściach średnia jest błędna, ponieważ mianownik średniej sum(x) / length(x)jest błędny, jeśli przefiltrujesz x, a następnie przekażesz go do mean.
Paul McMurdie
Myślę, że filtrowanie jest złym pomysłem, chyba że wyraźnie zamierzasz to zrobić (np. Gdybym pisał funkcję ogólnego przeznaczenia , nie ustawiłbym filtrowania jako domyślnego) - OK, jeśli jest to jednorazowy fragment kodu i masz myślał o tym, co bardzo dokładnie filtrowany zera spośród faktycznie oznacza w kontekście problemu (!)
Ben Bolker
Z definicji średnia geometryczna zbioru liczb zawierających zero powinna wynosić zero! math.stackexchange.com/a/91445/221143
Chris
6

Używam dokładnie tego, co mówi Mark. W ten sposób, nawet z tapply, możesz korzystać z wbudowanej meanfunkcji, bez konieczności definiowania swojej! Na przykład, aby obliczyć średnie geometryczne wartości $ danych dla poszczególnych grup:

exp(tapply(log(data$value), data$group, mean))
TMS
źródło
3

Ta wersja zawiera więcej opcji niż inne odpowiedzi.

  • Pozwala użytkownikowi rozróżnić wyniki, które nie są (rzeczywistymi) liczbami, i te, które nie są dostępne. Jeśli obecne są liczby ujemne, odpowiedź nie będzie liczbą rzeczywistą, więc NaNjest zwracana. Jeśli są to wszystkie NAwartości, funkcja zwróci NA_real_zamiast tego odzwierciedlenie, że rzeczywista wartość jest dosłownie niedostępna. Jest to subtelna różnica, ale taka, która może dać (nieco) solidniejsze wyniki.

  • Pierwszy opcjonalny parametr zero.rmma na celu zezwolenie użytkownikowi na to, aby zera wpływały na dane wyjściowe bez zerowania go. Jeśli zero.rmjest ustawiona na FALSEi etajest ustawiona na NA_real_(wartość domyślna), zera zmniejszają wynik do jedynki. Nie mam na to żadnego teoretycznego uzasadnienia - po prostu wydaje się, że bardziej sensowne jest nie ignorowanie zer, ale „zrobienie czegoś”, co nie obejmuje automatycznego zerowania wyniku.

  • etato sposób postępowania z zerami zainspirowany następującą dyskusją: https://support.bioconductor.org/p/64014/

geomean <- function(x,
                    zero.rm = TRUE,
                    na.rm = TRUE,
                    nan.rm = TRUE,
                    eta = NA_real_) {
    nan.count <- sum(is.nan(x))
     na.count <- sum(is.na(x))
  value.count <- if(zero.rm) sum(x[!is.na(x)] > 0) else sum(!is.na(x))

  #Handle cases when there are negative values, all values are missing, or
  #missing values are not tolerated.
  if ((nan.count > 0 & !nan.rm) | any(x < 0, na.rm = TRUE)) {
    return(NaN)
  }
  if ((na.count > 0 & !na.rm) | value.count == 0) {
    return(NA_real_)
  }

  #Handle cases when non-missing values are either all positive or all zero.
  #In these cases the eta parameter is irrelevant and therefore ignored.
  if (all(x > 0, na.rm = TRUE)) {
    return(exp(mean(log(x), na.rm = TRUE)))
  }
  if (all(x == 0, na.rm = TRUE)) {
    return(0)
  }

  #All remaining cases are cases when there are a mix of positive and zero
  #values.
  #By default, we do not use an artificial constant or propagate zeros.
  if (is.na(eta)) {
    return(exp(sum(log(x[x > 0]), na.rm = TRUE) / value.count))
  }
  if (eta > 0) {
    return(exp(mean(log(x + eta), na.rm = TRUE)) - eta)
  }
  return(0) #only propagate zeroes when eta is set to 0 (or less than 0)
}
Chris Coffee
źródło
2
Czy możesz dodać kilka szczegółów wyjaśniających, czym różni się to od istniejących rozwiązań / poprawia je? (Osobiście nie chciałbym dodawać ciężkiej zależności, jak w dplyrprzypadku takiego narzędzia, chyba że to konieczne ...)
Ben Bolker
Zgadzam się, te case_whenbyły trochę głupie, więc usunąłem je i zależność na korzyść ifs. Dostarczyłem też pewnych wyjaśnień.
Chris Coffee
1
Poszedłem z twoim drugim pomysłem i zmieniłem domyślne ustawienie nan.rmna, TRUEaby wyrównać wszystkie trzy parametry `` .rm ''.
Chris Coffee
1
Jeszcze jeden stylistyczny chwytak. ifelsejest przeznaczony do wektoryzacji. Z jednym warunkiem do sprawdzenia, byłoby bardziej idiomatyczne w użyciuvalue.count <- if(zero.rm) sum(x[!is.na(x)] > 0) else sum(!is.na(x))
Gregor Thomas
Wygląda ładniej niż ifelseteż. Zmieniony. Dzięki!
Chris Coffee
3

Jeśli w danych brakuje wartości, nie jest to rzadki przypadek. musisz dodać jeszcze jeden argument.

Możesz spróbować następującego kodu:

exp(mean(log(i[ is.finite(log(i)) ]), na.rm = TRUE))
Tian Yi
źródło
1
exp(mean(log(x1))) == prod(x1)^(1/length(x1))
user12882764
źródło