Расчет робастных процентилей в R

Question

У меня есть набор данных с параметром (тэта) и плотностью параметра (гхат). Я пытаюсь вычислить 10-й и 90-й процентили.

Проблема в том, что мои 90-ые процентили меньше среднего. Что я делаю неправильно? Есть ли лучший способ получить процентили с учетом имеющихся у меня данных?

df <- structure(list(theta = c(0.01, 0.02, 0.03, 0.04, 0.05, 0.06,
                         0.07, 0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17,
                         0.18, 0.19, 0.2, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28,
                         0.29, 0.3, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39,
                         0.4, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5,
                         0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6, 0.61,
                         0.62, 0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.7, 0.71, 0.72,
                         0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.8, 0.81, 0.82, 0.83,
                         0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94,
                         0.95, 0.96, 0.97, 0.98, 0.99), ghat = c(8.0497665821822e-10,
                                                                 1.26230992655857e-06, 8.26438628871083e-05, 0.00145824530878958,
                                                                 0.0124065393081063, 0.0658007731995211, 0.249038975036493, 0.727932651403715,
                                                                 1.72673415548062, 3.43406909516598, 5.85358700869572, 8.68314227655718,
                                                                 11.3284276546002, 13.094173436367, 13.4757519628881, 12.38746868323,
                                                                 10.1900555710108, 7.5111077282492, 4.9909183994075, 3.01639625008116,
                                                                 1.67228183293732, 0.85731428205356, 0.409565196187043, 0.183684664661787,
                                                                 0.0778924799288354, 0.031448806214798, 0.0121711651792063, 0.00454509949213374,
                                                                 0.00164831719139466, 0.000584215715883644, 0.000203629809008254,
                                                                 7.02266909309361e-05, 2.41086667080809e-05, 8.28732776812513e-06,
                                                                 2.86092948666336e-06, 9.90307440023348e-07, 3.43025256131075e-07,
                                                                 1.18650514458439e-07, 4.0894730452215e-08, 1.40139045176172e-08,
                                                                 4.76385534376833e-09, 1.60270281245794e-09, 5.32357991008131e-10,
                                                                 1.74159006347556e-10, 5.597419206501e-11, 1.76283038485334e-11,
                                                                 5.42584004312429e-12, 1.62772510314568e-12, 4.74620239243908e-13,
                                                                 1.3412948371487e-13, 3.66440610180999e-14, 9.66671282421902e-15,
                                                                 2.46018105914716e-15, 6.03455452969417e-16, 1.42510511150036e-16,
                                                                 3.23637808265161e-17, 7.05860980817371e-18, 1.47641726404577e-18,
                                                                 2.95700659895129e-19, 5.66118457054648e-20, 1.03409880677609e-20,
                                                                 1.79856703040568e-21, 2.97184153530068e-22, 4.65359973268091e-23,
                                                                 6.88718470220651e-24, 9.6048794951896e-25, 1.25807242160492e-25,
                                                                 1.54158396409635e-26, 1.7591743130845e-27, 1.86026900777325e-28,
                                                                 1.81301776633084e-29, 1.61875794813682e-30, 1.31532955754031e-31,
                                                                 9.65524286954306e-33, 6.35046336266651e-34, 3.70829126953137e-35,
                                                                 1.90275953781146e-36, 8.47940966417498e-38, 3.23858070076557e-39,
                                                                 1.04413668032609e-40, 2.79241430850164e-42, 6.07033401824646e-44,
                                                                 1.04749686096541e-45, 1.39538655335294e-47, 1.38828675388999e-49,
                                                                 9.91227489676626e-52, 4.83700986305667e-54, 1.51838997861272e-56,
                                                                 2.84052131966389e-59, 2.87016663788723e-62, 1.37611003830011e-65,
                                                                 2.62683062482798e-69, 1.56095716947688e-73, 2.01382058534212e-78,
                                                                 3.22597830727024e-84, 2.51631532528489e-91, 1.65259162698972e-100,
                                                                 1.77718300016021e-113, 1.08471354600224e-135)), row.names = c(NA,
                                                                                                                               -99L), class = c("tbl_df", "tbl", "data.frame"))

plot(df)

df$theta[which(df$ghat == max(df$ghat))]      # mle
#> [1] 0.15
df$theta[which(cumsum(df$ghat) >= .1)[1]]     # 10th percentile
#> [1] 0.07
df$theta[which(cumsum(df$ghat) >= .9)[1]]     # 90th percentile
#> [1] 0.08

^{Создано в 2020-02-12 с помощью представляет пакет (v0.3.0)}

Answer 1

Прежде всего, если вы сложите свой столбец df$ghat, вы увидите, что он суммирует до 100. Таким образом, ваш 10-й квантиль должен быть >= 10, а ваш 90-й >= 90.

. После этого кажется, что вы в основном оцениваете квантили для сгруппированных данных. Таким образом, способ, которым вы рассчитываете это немного отличается от обычного способа. Вы можете увидеть формулу в математика это весело в легко усваиваемом формате в разделе «Оценка медианы из сгруппированных данных»: https://www.mathsisfun.com/data/frequency-grouped-mean-median-mode.html.

Я мог бы объяснить это своими словами о том, что они говорят на веб-сайте, но там так ясно, что я не хочу разрушать это на самом деле. Чтобы вычислить 10-й и 90-й квантили, вы можете преобразовать их формулу в:

10-й квантиль

90-й квантиль

где ( полужирное"преобразование" переменных в их сообщение в ваш случай):

• L - нижняя граница класса группы, содержащей квантиль
• n - сумма плотности (в данном случае 100)
• B - совокупная плотность групп перед квантильной группой
• G - плотность группы
• w - ширина группы

Обычно в целочисленных группах они размещают нижнюю границу и верхнюю граница группы между элементами, что также выполнимо в вашем случае. В случае группы с theta = 0.10, то есть после theta = 0.09 и до theta = 0.11, соответствующая нижняя граница будет:

ширина w, упомянутая в формулу также можно рассчитать, исходя из границ группы. Учитывая вышеупомянутый пример, нижние границы тета 0.10 и 0.11 соответственно равны L(0.10) = 0.095 и L(0.11) = 0.105, поэтому ширина первой группы будет

Answer 2

10% квантиль

df %>% 
  mutate(cdf=cumsum(ghat)/sum(ghat)) %>%
  filter(cdf <= 10) %>%
  summarise(max(theta))

90% квантиль

df %>% 
  mutate(cdf=cumsum(ghat)/sum(ghat)) %>%
  filter(cdf >= 0.9) %>%
  summarise(min(theta))

ожидаемое значение

df %>% 
  summarise(weighted.mean(theta,ghat))