我正在嘗試繪制自定義 CDF 函式的 qq 函式，形式為 f(x)= p*e^-x (1-p)*e^-5x，p 介于 0 和 1 之間，并且 x 更大比 0. 資料具有值 0<x<500,000
在這樣做時，我偶然發現了 uniroot 函式，它旨在解決反函式。然而，我遇到了一些錯誤，非常感謝這方面的幫助。
我的目標是按照 R 的基本函式qnorm qexp等撰寫一個函式，但使用qcustom.
# Error 1: (Error in f(x) : argument "p" is missing, with no default)
inverse <- function(f, lower=0.01, upper=500000) {
function(y) {
uniroot(function(x) f(x) - y, lower=lower, upper=upper)[1]
}
}
qcustom <- inverse(function(x, p, a, b) {
a1 <- p*exp(-a*x)
a2 <- (1 - p)*exp(-b*x)
result <- 1 - a1 - a2
return(result)
})
qcustom(0.5) ##when i add in further parameters, e.g p = 0.5 it says it's an unused parameter##
uj5u.com熱心網友回復：
讓我們從頭開始：
custom_cdf <- function(x, p, a, b) {
p*exp(-a*x) (1-p)*exp(-b*x)
}
這有引數c(p, a, b)，讓我們嘗試找到一個對你指出你想要的c(a,b)有意義的配對range(x) = c(0, 500e3)。
plot.new()
curve(custom_cdf(x, p = 0.5, a = 1e-4, b = 1e-5), xlim = c(0, 500e3))
一組合理的值可以是：
a <- 1e-4
b <- 1e-5
讓我們嘗試多個值p：
first <- TRUE
plot.new()
for (p in seq.default(0.1, 1, length.out = 25)) {
curve(custom_cdf(x, p = p, a = 1e-4, b = 1e-5), xlim = c(0, 500e3), add = !first)
first <- FALSE
}
rm('p')
對于 CDF 的倒數，我們可以：
xx <- seq.default(0, 500e3, length.out = 2500)
# xx <- seq.default(0, 500e3, by = 10)
# xx <- seq.default(0, 500e3, by = 1)
yy <- custom_cdf(xx, p = 0.5, a = a, b = b)
plot.new()
plot(yy, xx, main = "inverse of CDF = Quantile Function",type = 'l')
I.e. we can just reverse the plotting. So now, this is a picture of the inverse function.
Let's use uniroot to get one inverse value.
uniroot(
function(x) custom_cdf(x, p = 0.5, a = a, b = b) - 0.2,
interval = c(0, 100e50),
extendInt = "yes"
) -> found_20pct_point
#'
#' Let's plot that point:
#'
points(custom_cdf(found_20pct_point[[1]], p = 0.5, a = a, b = b), found_20pct_point[[1]])
Alright, so now we understand everything and we can do what you tried to do:
custom_quantile <- function(p, a, b) {
particular_cdf <- function(x) custom_cdf(x, p, a, b)
function(prob) {
uniroot(
function(x) particular_cdf(x) - prob,
interval = c(0, 100e50),
extendInt = "yes"
)$root
}
}
custom_quantile_specific <-
custom_quantile(p = 0.5, a = a, b = b)
custom_quantile_specific(0.2)
And this works as it yields something close to the desired value.
But regular quantile functions in R does depend on the parameters of the distribution, so this is good enough:
custom_quantile <- function(prob, p, a, b) {
particular_cdf <- function(x) custom_cdf(x, p, a, b)
uniroot(
function(x) particular_cdf(x) - prob,
interval = c(0, 100e50),
extendInt = "yes"
)$root
}
custom_quantile(0.2, p = 0.5, a = a, b = b)
轉載請註明出處，本文鏈接：https://www.uj5u.com/ruanti/334667.html
標籤：r
上一篇：將行轉移到r中的列