Geometric Mean: apakah ada built-in?

106

Saya mencoba menemukan built-in untuk mean geometris tetapi tidak bisa.

(Jelas built-in tidak akan menyelamatkan saya kapan saja saat bekerja di shell, saya juga tidak curiga ada perbedaan dalam akurasi; untuk skrip saya mencoba menggunakan built-in sesering mungkin, di mana (kumulatif) peningkatan kinerja sering kali terlihat.

Jika tidak ada satu (yang saya ragu adalah kasusnya), ini milik saya.

gm_mean = function(a){prod(a)^(1/length(a))}
doug
sumber
11
Hati-hati dengan angka negatif dan luapan. prod (a) akan turun atau meluap dengan sangat cepat. Saya mencoba mengatur waktu ini menggunakan daftar besar dan dengan cepat mendapatkan Inf menggunakan metode Anda vs 1.4 dengan exp (mean (log (x))); masalah pembulatan bisa sangat parah.
Tristan
saya baru saja menulis fungsi di atas dengan cepat karena saya yakin bahwa 5 menit setelah memposting Q ini, seseorang akan memberi tahu saya R built-in untuk gm. Jadi tidak ada built-in jadi ada baiknya meluangkan waktu untuk membuat kode ulang sehubungan dengan komentar Anda. +1 dari saya.
doug
1
Saya baru saja menandai ini geometris-mean dan built-in , 9 tahun kemudian.
smci

Jawaban:

77

Berikut adalah fungsi vectorized, zero- dan NA-tolerant untuk menghitung rata-rata geometris di R. meanPerhitungan verbose yang melibatkan length(x)diperlukan untuk kasus-kasus di mana xmengandung nilai-nilai non-positif.

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

Terima kasih kepada @ ben-bolker karena telah mencatat na.rmpass-through dan @Gregor untuk memastikannya berfungsi dengan benar.

Saya pikir beberapa komentar terkait dengan kesetaraan NAnilai palsu dalam data dan nol. Dalam penerapan yang saya pikirkan, mereka sama, tetapi tentu saja ini tidak benar secara umum. Jadi, jika Anda ingin menyertakan penyebaran opsional dari nol, dan memperlakukan secara length(x)berbeda dalam kasus NApenghapusan, berikut ini adalah alternatif yang sedikit lebih panjang untuk fungsi di atas.

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

Perhatikan bahwa ini juga memeriksa nilai negatif apa pun, dan mengembalikan nilai yang lebih informatif dan tepat NaNsehubungan dengan rata-rata geometris tidak ditentukan untuk nilai negatif (tetapi untuk nol). Terima kasih kepada pemberi komentar yang tetap menangani kasus saya tentang ini.

Paul McMurdie
sumber
2
bukankah lebih baik untuk diteruskan na.rmsebagai argumen (yaitu biarkan pengguna memutuskan apakah mereka ingin toleran NA atau tidak, untuk konsistensi dengan fungsi ringkasan R lainnya)? Saya gugup tentang secara otomatis mengecualikan nol - saya akan menjadikannya sebagai pilihan juga.
Ben Bolker
1
Mungkin Anda benar tentang lulus na.rmsebagai opsi. Saya akan memperbarui jawaban saya. Adapun untuk mengecualikan nol, rata-rata geometris tidak ditentukan untuk nilai non-positif, termasuk nol. Di atas adalah fiksasi umum untuk rata-rata geometrik, di mana nol (atau dalam hal ini semua bukan nol) diberi nilai dummy 1, yang tidak berpengaruh pada hasil kali (atau ekuivalen, nol dalam jumlah logaritmik).
Paul McMurdie
* Maksud saya perbaikan umum untuk nilai non-positif, nol menjadi yang paling umum saat rata-rata geometris digunakan.
Paul McMurdie
1
Anda na.rmpass-through tidak bekerja sebagai kode ... lihat gm_mean(c(1:3, NA), na.rm = T). Anda perlu menghapus & !is.na(x)dari subset vektor, dan karena argumen pertama sumadalah ..., Anda harus memberikan na.rm = na.rmnama, dan Anda juga perlu mengecualikan 0's dan NA' dari vektor saat lengthpanggilan.
Gregor Thomas
2
Hati-hati: karena xhanya mengandung nol, seperti x <- 0, exp(sum(log(x[x>0]), na.rm = TRUE)/length(x))memberikan 1mean geometrik, yang tidak masuk akal.
adatum
88

Tidak, tapi ada beberapa orang yang pernah menulisnya, seperti di sini .

Kemungkinan lain adalah menggunakan ini:

exp(mean(log(x)))
Mark Byers
sumber
Keuntungan lain menggunakan exp (mean (log (x))) adalah Anda dapat bekerja dengan daftar panjang bilangan besar, yang bermasalah saat menggunakan rumus yang lebih jelas menggunakan prod (). Perhatikan bahwa prod (a) ^ (1 / length (a)) dan exp (mean (log (a))) memberikan jawaban yang sama.
lukeholman
tautan telah diperbaiki
PatrickT
15

Kita bisa menggunakan paket psych dan memanggil fungsi geometric.mean .

AliCivil
sumber
1
psych::geometric.mean()
smci
Fungsi-fungsi ini harus mengambil rangkaian dan bukan pertumbuhannya, setidaknya sebagai pilihan, menurut saya.
Christoph Hanck
12

Itu

exp(mean(log(x)))

akan bekerja kecuali ada 0 di x. Jika demikian, log akan menghasilkan -Inf (-Infinite) yang selalu menghasilkan mean geometrik 0.

Salah satu solusinya adalah menghapus nilai -Inf sebelum menghitung mean:

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

Anda dapat menggunakan satu baris untuk melakukan ini tetapi itu berarti menghitung log dua kali yang tidak efisien.

exp(mean(log(i[is.finite(log(i))])))
Alan James Salmoni
sumber
mengapa menghitung log dua kali ketika Anda bisa melakukan: exp (mean (x [x! = 0]))
zzk
kedua pendekatan mendapatkan mean salah, karena penyebut untuk mean, sum(x) / length(x)salah jika Anda memfilter x dan meneruskannya mean.
Paul McMurdie
Saya pikir pemfilteran adalah ide yang buruk kecuali Anda secara eksplisit bermaksud melakukannya (misalnya jika saya menulis fungsi tujuan umum, saya tidak akan menjadikan pemfilteran sebagai default) - OK jika ini adalah potongan kode satu kali dan Anda Pikirkan dengan sangat hati-hati tentang apa sebenarnya arti memfilter nol dalam konteks masalah Anda (!)
Ben Bolker
Menurut definisi, rata-rata geometris dari sekumpulan angka yang mengandung nol harus nol! math.stackexchange.com/a/91445/221143
Chris
6

Saya menggunakan persis apa yang dikatakan Mark. Dengan cara ini, bahkan dengan tapply, Anda dapat menggunakan meanfungsi bawaan , tidak perlu menentukan milik Anda! Misalnya, untuk menghitung rata-rata geometris per grup dari data $ value:

exp(tapply(log(data$value), data$group, mean))
TMS
sumber
3

Versi ini memberikan lebih banyak pilihan daripada jawaban lainnya.

  • Ini memungkinkan pengguna untuk membedakan antara hasil yang bukan bilangan (nyata) dan yang tidak tersedia. Jika ada angka negatif, jawabannya bukan bilangan real, jadi NaNdikembalikan. Jika itu semua NAnilai maka fungsinya akan kembali NA_real_untuk mencerminkan bahwa nilai sebenarnya secara harfiah tidak tersedia. Ini adalah perbedaan yang halus, tetapi mungkin menghasilkan (sedikit) hasil yang lebih kuat.

  • Parameter opsional pertama zero.rmdimaksudkan agar pengguna memiliki nol yang memengaruhi keluaran tanpa menjadikannya nol. Jika zero.rmdisetel ke FALSEdan etadisetel ke NA_real_(nilai defaultnya), nol memiliki efek menyusutkan hasilnya ke satu. Saya tidak memiliki pembenaran teoretis untuk ini - sepertinya lebih masuk akal untuk tidak mengabaikan angka nol tetapi untuk "melakukan sesuatu" yang tidak melibatkan otomatis membuat hasil menjadi nol.

  • etaadalah cara menangani angka nol yang terinspirasi dari diskusi berikut: 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
sumber
1
Dapatkah Anda menambahkan beberapa detail yang menjelaskan bagaimana hal ini berbeda dari / meningkatkan solusi yang ada? (Saya pribadi tidak ingin menambahkan ketergantungan yang berat seperti dplyruntuk utilitas seperti itu kecuali jika diperlukan ...)
Ben Bolker
Saya setuju, case_whens itu agak konyol, jadi saya menghapusnya dan ketergantungan mendukung ifs. Saya juga memberikan beberapa elaborasi.
Chris Coffee
1
Saya pergi dengan ide terakhir Anda dan mengubah default nan.rmmenjadi TRUEuntuk menyelaraskan ketiga parameter ".rm``.
Chris Coffee
1
Satu nitpick gaya lainnya. ifelsedirancang untuk vektorisasi. Dengan satu kondisi untuk diperiksa, akan lebih idiomatis untuk digunakanvalue.count <- if(zero.rm) sum(x[!is.na(x)] > 0) else sum(!is.na(x))
Gregor Thomas
Ini juga terlihat lebih bagus dari ifelse. Berubah. Terima kasih!
Chris Coffee
3

Jika ada nilai yang hilang dalam data Anda, ini bukan kasus yang jarang terjadi. Anda perlu menambahkan satu argumen lagi.

Anda dapat mencoba kode berikut:

exp(mean(log(i[ is.finite(log(i)) ]), na.rm = TRUE))
Tian Yi
sumber
1
exp(mean(log(x1))) == prod(x1)^(1/length(x1))
pengguna12882764
sumber