Построение функций CDF и квантилей с учетом PDF

Как бы я построил функции CDF и Quantile в R, если у меня есть PDF. В настоящее время у меня есть следующее (но я думаю, что должен быть лучший способ сделать это):

## Probability Density Function
p <- function(x) {
  result <- (x^2)/9
  result[x < 0 | x > 3] <- 0
  result
}

plot(p, xlim = c(0,3), main="Probability Density Function")

## Cumulative Distribution Function
F <- function(a = 0,b){
  result <- ((b^3)/27) - ((a^3)/27)
  result[a < 0 ] <- 0
  result[b > 3] <- 1
  result
}

plot(F(,x), xlim=c(0,3), main="Cumulative Distribution Function")

## Quantile Function
Finv <- function(p) {
  3*x^(1/3)
}

person Amanda R.    schedule 21.03.2017    source источник
comment
Возможно полезно: stats::integrate   -  person    schedule 21.03.2017
comment
Вы также можете посмотреть на library(Ryacas) или library(Rsympy)   -  person    schedule 21.03.2017
comment
Как-то непонятно, чего вы добиваетесь. Я имею в виду, похоже, вы уже знаете, что для перехода от pdf к cdf вам нужно интегрироваться. Не все функции легко интегрируются. Ваша функция имеет аналитическое решение. Это значение, которое вы хотите вернуть? Или вас просто интересуют численные приближения для любого pdf в целом? R на самом деле не понимает статистику; просто в него встроено множество статистических функций.   -  person MrFlick    schedule 21.03.2017
comment
@MrFlick Мне интересно, есть ли лучший способ найти CDF. Я также ищу функцию квантиля.   -  person Amanda R.    schedule 21.03.2017
comment
@MrFlick @dash2 Что я должен указать для значения по умолчанию b?   -  person Amanda R.    schedule 21.03.2017
comment
Лучше, чем что? Уже знаете CDF? Это в значительной степени лучший сценарий - когда у вас есть хорошая закрытая от для вашего CDF. Base R не может помочь вам с этим. Это не очень хорошо для символической алгебры (но некоторые из пакетов, рекомендованных @dash2, могут использоваться таким образом). Что для вас значит лучше?   -  person MrFlick    schedule 21.03.2017
comment
@MrFlick Какое значение по умолчанию для b я должен установить?   -  person Amanda R.    schedule 21.03.2017
comment
@MrFlick Подходит ли это для квантильной функции?   -  person Amanda R.    schedule 21.03.2017
comment
Для b нет значения по умолчанию. Это зависит от того, для какого региона вы хотите рассчитать площадь. Может быть, переместить b в качестве первого параметра и оставить a необязательным? И что касается функции квантили, мы здесь не для того, чтобы проверять вашу домашнюю работу по статистике. Если у вас есть конкретный вопрос или проблема, связанная с программированием, сообщите об этом. Дайте образец ввода и желаемый результат.   -  person MrFlick    schedule 21.03.2017


Ответы (1)


Как предложил @dash2, CDF потребуется, чтобы вы интегрировали PDF, по сути, вам нужно найти область под кривой.

Вот общее решение, которое должно помочь. В качестве примера я использую распределение Гаусса — вы должны иметь возможность передать ему любую общую функцию.

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

# CDF Function
CDF <- function(FUNC = p, plot = T, area = 0.5, LOWER = -10, UPPER = 10, SIZE = 1000){

    # Create data
    x <- seq(LOWER, UPPER, length.out = SIZE)
    y <- p(x)

    area.vec <- c()
    area.vec[1] <- 0

    for(i in 2:length(x)){
        x.vec <- x[1:i]
        y.vec <- y[1:i]

        area.vec[i] = integrate(p, lower = x[1], upper = x[i])$value
    }

    # Quantile
    quantile = x[which.min(abs(area.vec - area))]

    # Plot if requested
    if(plot == TRUE){

        # PDF
        par(mfrow = c(1, 2))
        plot(x, y, type = "l", main = "PDF", col = "indianred", lwd = 2)
        grid()

        # CDF
        plot(x, area.vec, type = "l", main = "CDF", col = "slateblue",
             xlab = "X", ylab = "CDF", lwd = 2)

        # Quantile 
        mtext(text = paste("Quantile at ", area, "=",
                           round(quantile, 3)), side = 3)
        grid()

        par(mfrow = c(1, 1))
    }
}

# Sample data
# PDF Function - Gaussian distribution
p <- function(x, SD = 1, MU = 0){
    y <- (1/(SD * sqrt(2*pi)) * exp(-0.5 * ((x - MU)/SD) ^ 2))
    return(y)
}

# Call to function
CDF(p, area = 0.5, LOWER = -5, UPPER = 5)

введите здесь описание изображения

person royr2    schedule 21.03.2017