关键词:R语言,经验分布函数,生存函数
经验分布函数的数学定义:
R自定义函数--经验分布函数:
scores <- c(63, 72, 74, 79, 82, 82, 87, 89, 90, 90)
# R的自定义函数,函数名为cdf_table
cdf_table <- function(x){
x <- sort(x) #将x按照从小到大的顺序排列,不删掉重复的值。
n <- length(x) #x的长度
tab <- unname(c(table(x))) #x每个数出现的次数,即频数。
pct = tab / n #计算频率
d <- data.frame( #形成一个数据框
x = sort(unique(x)), #删掉x中重复的值,确保x中每个数都是唯一的,然后对这些唯一的数,按照从小到大的顺序进行排列
freq = tab, #变量名为freq,该列中给出x中数值出现的频数(出现的次数)
pct = pct, #变量名为pct,该列给出x中数值对应的频率
cumfreq = cumsum(tab), #变量名为cumfreq,该列中给出频数的累加和,即,该行的第二个值=tab的第一个值,加上 tab的第二个值。
cumpct = cumsum(pct)) #变量名为cumpct,该列中给出频率的累加和。
d #该函数返回d(数据框)
}
knitr::kable(cdf_table(scores)) #生成表格
得到如下表格:
其中
- x表示原始数据(1)先unique,即删掉重复的数值(2)再将剩下的数值按照从小到大的顺序进行排列。
- freq表示x对应的频数;
- pct表示x对应的频率;
- cumfreq表示x对应的频数累加和;
- cumpct表示x对应的频率累加和。
绘制经验分布函数:
d <- data.frame(
scores = c(63, 72, 74, 79, 82, 82, 87, 89, 90, 90)) #原始数据,数据类型为数据框;
dfreq <- cdf_table(d$scores) #返回数据对应的数据框,包含原始数据的取值,每个值对应的频数,频率,频数累加和,频率累加和
p <- ggplot(data = dfreq, mapping = aes(
x = x, y = cumpct )) #基于作图的数据为dfreq(数据框), mapping给出了数据要映射的x轴和y轴数据
#对上图做进一步细化设定
p + stat_ecdf(geom = "step") +
scale_y_continuous(
limits = c(-.01, 1.01), #调整y轴坐标轴的显示范围
expand = c(0, 0), #expand = c(0,0)不放缩,这个命令自己可以尝试修改不同的数字,看下图的效果
name = " 累计百分比") #y轴显示的标签名称
绘制上图对应的生存函数:
surv_table <- function(x){
x <- sort(x)
n <- length(x)
tab <- unname(c(table(x)))
pct = tab / n
d <- data.frame(
x = sort(unique(x)),
freq = tab,
pct = pct,
cumfreq = cumsum(tab),
cumpct = cumsum(pct),
cumsurv = 1-cumsum(pct))
d
}
knitr::kable(surv_table(scores))
scores = c(63, 72, 74, 79, 82, 82, 87, 89, 90, 90)
surv_req <-surv_table(scores)
f1 <- c(-Inf,0,0,0,0,1) # 手动增加-Inf处的取值
f2 <- c(Inf,0,0,0,1,0) #手动增加Inf处的取值
surv_req <- rbind(f1,surv_req,f2)
p <- ggplot(data = surv_req, mapping = aes(
x = x, y = cumsurv ))
p + geom_step() +
scale_y_continuous(
limits = c(-.01, 1.01),
expand = c(0, 0),
name = " 累计百分比")
在绘制生存函数的时候,这里使用了geom_step()函数,需要补充上,阶梯函数在-Inf和Inf处的取值
其他一些可以参考的知识:
关于ggplot函数,可以参考北京大学李东风老师编写的Rbook。下图是李东风老师关于ggplot作图的介绍:
R-base包中的table()函数:
table(x)返回x的每个不同值的频率(出现次数),table函数用于创建列联表,用于统计x中每个不同取值的频率。返回的是一个table类。
unname()函数,用于删除R对象中,变量的名称。
## Stat_ecdf 原函数
stat_ecdf <- function(mapping = NULL,
data = NULL,
geom = "step",
position = "identity",
...,
n = NULL,
pad = TRUE,
na.rm = FALSE,
show.legend = NA,
inherit.aes = TRUE) {
layer(
data = data,
mapping = mapping,
stat = StatEcdf,
geom = geom,
position = position,
show.legend = show.legend,
inherit.aes = inherit.aes,
params = list2(
n = n,
pad = pad,
na.rm = na.rm,
...
)
)
}
#' @rdname ggplot2-ggproto
#' @format NULL
#' @usage NULL
#' @export
StatEcdf <- ggproto("StatEcdf", Stat,
required_aes = c("x|y"),
default_aes = aes(x = after_stat(ecdf), y = after_stat(ecdf)),
setup_params = function(self, data, params) {
params$flipped_aes <- has_flipped_aes(data, params, main_is_orthogonal = FALSE, main_is_continuous = TRUE)
has_x <- !(is.null(data$x) && is.null(params$x))
has_y <- !(is.null(data$y) && is.null(params$y))
if (!has_x && !has_y) {
cli::cli_abort("{.fn {snake_class(self)}} requires an {.field x} or {.field y} aesthetic.")
}
params
},
compute_group = function(data, scales, n = NULL, pad = TRUE, flipped_aes = FALSE) {
data <- flip_data(data, flipped_aes)
# If n is NULL, use raw values; otherwise interpolate
if (is.null(n)) {
x <- unique0(data$x)
} else {
x <- seq(min(data$x), max(data$x), length.out = n)
}
if (pad) {
x <- c(-Inf, x, Inf)
}
data_ecdf <- ecdf(data$x)(x)
df_ecdf <- data_frame0(
x = x,
y = data_ecdf,
ecdf = data_ecdf,
.size = length(x)
)
df_ecdf$flipped_aes <- flipped_aes
flip_data(df_ecdf, flipped_aes)
}
)
推荐常用R语言包及其函数检索网址:R Package Documentation
李东风老师R书,R语言教程 (pku.edu.cn)