Объясните функцию quantile () в R

Существуют различные способы вычисления квантилей, когда вы задаете ему вектор и не имеете известной функции CDF.

Обдумайте вопрос, что делать, если ваши наблюдения не попадают точно в квантили.

«Типы» только определяют, как это сделать. Итак, методы говорят: «использовать линейную интерполяцию между статистикой k-го порядка и p (k)».

Итак, что такое p (k)? Один парень говорит: «Мне нравится использовать k / n». Другой парень говорит: «Мне нравится использовать (k-1) / (n-1)» и т. д. Каждый из этих методов имеет разные свойства, которые лучше подходят для той или иной проблемы.

\ Alpha и \ beta - это просто способы параметризации функций p. В одном случае это 1 и 1. В другом случае это 3/8 и -1/4. Я не думаю, что буква p когда-либо является константой в документации. Просто они не всегда явно показывают зависимость.

Посмотрите, что происходит с разными типами, когда вы вставляете такие векторы, как 1: 5 и 1: 6.

(также обратите внимание, что даже если ваши наблюдения попадают точно в квантили, некоторые типы все равно будут использовать линейную интерполяцию).

Вы по понятным причинам запутались. Эта документация ужасна. Мне пришлось вернуться к статье, основанной на (Hyndman, R.J .; Fan, Y. (ноябрь 1996). «Sample Quantiles in Statistical Packages». Американский статистик 50 (4): 361–365. DOI: 10.2307 / 2684934), чтобы понять. Начнем с первой проблемы.

where 1 <= i <= 9, (j-m)/n <= p < (j-m+1)/ n, x[j] is the jth order statistic, n is the sample size, and m is a constant determined by the sample quantile type. Here gamma depends on the fractional part of g = np+m-j.

Первая часть взята прямо из статьи, но авторы документации упустили j = int(pn+m). Это означает, что Q[i](p) зависит только от статистики двух порядков, наиболее близкой к тому, что p составляет часть пути через (отсортированные) наблюдения. (Для тех, кто, как и я, не знаком с этим термином, «статистика порядка» серии наблюдений - это отсортированная серия.)

Кроме того, последнее предложение просто неверно. Следует читать

Here gamma depends on the fractional part of np+m, g = np+m-j

Что касается m, то тут все просто. m зависит от того, какой из 9 алгоритмов был выбран. Так же, как Q[i] является функцией квантиля, m следует рассматривать как m[i]. Для алгоритмов 1 и 2 m равен 0, для 3 m равен -1/2, а для остальных это будет в следующей части.

For the continuous sample quantile types (4 through 9), the sample quantiles can be obtained by linear interpolation between the kth order statistic and p(k):

p(k) = (k - alpha) / (n - alpha - beta + 1), where α and β are constants determined by the type. Further, m = alpha + p(1 - alpha - beta), and gamma = g.

Это действительно сбивает с толку. То, что в документации называется p(k), отличается от предыдущей модели p. p(k) - это положение на графике. В статье авторы пишут это как pk, что помогает. Тем более, что в выражении для mp - это исходный p, а m = alpha + p * (1 - alpha - beta). Концептуально для алгоритмов 4–9 точки (pk, x[k]) интерполируются для получения решения (p, Q[i](p)). Каждый алгоритм отличается только алгоритмом для pk.

Что касается последнего бита, R просто указывает, что использует S.

В исходной статье приводится список из 6 «желаемых свойств для функции выборочного квантиля» и указывается, что предпочтение отдается №8, которое удовлетворяет всем с точностью до 1. №5 удовлетворяет всем им, но им это не нравится по другим причинам ( более феноменологичен, чем выведен из принципов). №2 - это то, что такие гики, как я, считают квантилями, и это то, что описано в Википедии.

Кстати, в ответ на Дрейв отвечает Mathematica делает вещи значительно иначе. Думаю, я понимаю отображение. Хотя систему Mathematica легче понять, (а) легче выстрелить себе в ногу с бессмысленными параметрами, и (б) она не может выполнить алгоритм №2 R. (Здесь Страница квантилей Mathworld, в котором говорится, что Mathematica не может выполнить №2, но дает более простое обобщение всех других алгоритмов в терминах четырех параметров.)

Я считаю, что справочная документация R ясна после исправлений, отмеченных в комментарии @ RobHyndman, но я нашел это немного подавляющим. Я отправляю этот ответ на случай, если он поможет кому-то быстро перейти к вариантам и их предположениям.

Чтобы разобраться в quantile(x, probs=probs), я хотел проверить исходный код. Это тоже было сложнее, чем я ожидал в R, поэтому я просто взял его из репозиторий github, который выглядел достаточно свежим, чтобы работать с ним. Меня интересовало поведение по умолчанию (тип 7), поэтому я аннотировал некоторые из них, но не делал того же для каждого варианта.

Вы можете увидеть, как метод «типа 7» интерполирует шаг за шагом как в коде, так и я добавил несколько строк, чтобы напечатать некоторые важные значения по ходу.

quantile.default <-function(x, probs = seq(0, 1, 0.25), na.rm = FALSE, names = TRUE
         , type = 7, ...){
    if (is.factor(x)) { #worry about non-numeric data
        if (!is.ordered(x) || ! type %in% c(1L, 3L))
            stop("factors are not allowed")
        lx <- levels(x)
    } else lx <- NULL
    if (na.rm){
        x <- x[!is.na(x)]
    } else if (anyNA(x)){
        stop("missing values and NaN's not allowed if 'na.rm' is FALSE")
        }
    eps <- 100*.Machine$double.eps #this is to deal with rounding things sensibly
    if (any((p.ok <- !is.na(probs)) & (probs < -eps | probs > 1+eps)))
        stop("'probs' outside [0,1]")

    #####################################
    # here is where terms really used in default type==7 situation get defined

    n <- length(x) #how many observations are in sample?

    if (na.p <- any(!p.ok)) { # set aside NA & NaN
        o.pr <- probs
        probs <- probs[p.ok]
        probs <- pmax(0, pmin(1, probs)) # allow for slight overshoot
    }

    np <- length(probs) #how many quantiles are you computing?

    if (n > 0 && np > 0) { #have positive observations and # quantiles to compute
        if (type == 7) { # be completely back-compatible

            index <- 1 + (n - 1) * probs #this gives the order statistic of the quantiles
            lo <- floor(index)  #this is the observed order statistic just below each quantile
            hi <- ceiling(index) #above
            x <- sort(x, partial = unique(c(lo, hi))) #the partial thing is to reduce time to sort, 
            #and it only guarantees that sorting is "right" at these order statistics, important for large vectors 
            #ties are not broken and tied elements just stay in their original order
            qs <- x[lo] #the values associated with the "floor" order statistics
            i <- which(index > lo) #which of the order statistics for the quantiles do not land on an order statistic for an observed value

            #this is the difference between the order statistic and the available ranks, i think
            h <- (index - lo)[i] # > 0  by construction 
            ##      qs[i] <- qs[i] + .minus(x[hi[i]], x[lo[i]]) * (index[i] - lo[i])
            ##      qs[i] <- ifelse(h == 0, qs[i], (1 - h) * qs[i] + h * x[hi[i]])
            qs[i] <- (1 - h) * qs[i] + h * x[hi[i]] # This is the interpolation step: assemble the estimated quantile by removing h*low and adding back in h*high. 
            # h is the arithmetic difference between the desired order statistic amd the available ranks
            #interpolation only occurs if the desired order statistic is not observed, e.g. .5 quantile is the actual observed median if n is odd. 
            # This means having a more extreme 99th observation doesn't matter when computing the .75 quantile


            ###################################
            # print all of these things

            cat("floor pos = ", c(lo))
            cat("\nceiling pos = ", c(hi))
            cat("\nfloor values= ", c(x[lo]))
            cat( "\nwhich floors not targets? ", c(i))
            cat("\ninterpolate between ", c(x[lo[i]]), ";", c(x[hi[i]]))
            cat( "\nadjustment values= ", c(h))
            cat("\nquantile estimates:")

    }else if (type <= 3){## Types 1, 2 and 3 are discontinuous sample qs.
                nppm <- if (type == 3){ n * probs - .5 # n * probs + m; m = -0.5
                } else {n * probs} # m = 0

                j <- floor(nppm)
                h <- switch(type,
                            (nppm > j),     # type 1
                            ((nppm > j) + 1)/2, # type 2
                            (nppm != j) | ((j %% 2L) == 1L)) # type 3

                } else{
                ## Types 4 through 9 are continuous sample qs.
                switch(type - 3,
                       {a <- 0; b <- 1},    # type 4
                       a <- b <- 0.5,   # type 5
                       a <- b <- 0,     # type 6
                       a <- b <- 1,     # type 7 (unused here)
                       a <- b <- 1 / 3, # type 8
                       a <- b <- 3 / 8) # type 9
                ## need to watch for rounding errors here
                fuzz <- 4 * .Machine$double.eps
                nppm <- a + probs * (n + 1 - a - b) # n*probs + m
                j <- floor(nppm + fuzz) # m = a + probs*(1 - a - b)
                h <- nppm - j

                if (any(sml <- abs(h) < fuzz)) h[sml] <- 0

            x <- sort(x, partial =
                          unique(c(1, j[j>0L & j<=n], (j+1)[j>0L & j<n], n))
            )
            x <- c(x[1L], x[1L], x, x[n], x[n])
            ## h can be zero or one (types 1 to 3), and infinities matter
            ####        qs <- (1 - h) * x[j + 2] + h * x[j + 3]
            ## also h*x might be invalid ... e.g. Dates and ordered factors
            qs <- x[j+2L]
            qs[h == 1] <- x[j+3L][h == 1]
            other <- (0 < h) & (h < 1)
            if (any(other)) qs[other] <- ((1-h)*x[j+2L] + h*x[j+3L])[other]

            } 
    } else {
        qs <- rep(NA_real_, np)}

    if (is.character(lx)){
        qs <- factor(qs, levels = seq_along(lx), labels = lx, ordered = TRUE)}
    if (names && np > 0L) {
        names(qs) <- format_perc(probs)
    }
    if (na.p) { # do this more elegantly (?!)
        o.pr[p.ok] <- qs
        names(o.pr) <- rep("", length(o.pr)) # suppress <NA> names
        names(o.pr)[p.ok] <- names(qs)
        o.pr
    } else qs
}

####################

# fake data
x<-c(1,2,2,2,3,3,3,4,4,4,4,4,5,5,5,5,5,5,5,5,5,6,6,7,99)
y<-c(1,2,2,2,3,3,3,4,4,4,4,4,5,5,5,5,5,5,5,5,5,6,6,7,9)
z<-c(1,2,2,2,3,3,3,4,4,4,4,4,5,5,5,5,5,5,5,5,5,6,6,7)

#quantiles "of interest"
probs<-c(0.5, 0.75, 0.95, 0.975)

# a tiny bit of illustrative behavior
quantile.default(x,probs=probs, names=F)
quantile.default(y,probs=probs, names=F) #only difference is .975 quantile since that is driven by highest 2 observations
quantile.default(z,probs=probs, names=F) # This shifts everything b/c now none of the quantiles fall on an observation (and of course the distribution changed...)... but 
#.75 quantile is stil 5.0 b/c the observations just above and below the order statistic for that quantile are still 5. However, it got there for a different reason.

#how does rescaling affect quantile estimates?
sqrt(quantile.default(x^2, probs=probs, names=F))
exp(quantile.default(log(x), probs=probs, names=F))