R语言入门——画密度曲线

0引言

给出下面得密度函数，

得出下面得密度图：

任务说明

给出相应得密度函数，还原上述图片，并给出置信区间。

1、案例

1.1代码

1.1.1包的载入以及密度函数的书写

library(ggplot2)  # 载入包
fun <- function(Thetahat, Theta = 0.5, n = 22){
  k = "/"("*"(
             (n-2)*((1-Theta^2)^((n-1)/2)),
             (1-Thetahat^2)^((n-4)/2)
             ),pi
         )
  f <- function(w) "/"(1,(cosh(w)-Theta*Thetahat)^(n-1))
  Int <- integrate(f,0,Inf)
  return(k*Int$val)
}

1.1.2 密度函数向量化

Gfun <- function(f = fun, x = x, Theta = 0.5){
  n = length(x)
  y = c()
  for(i in 1:n){
    y <- c(y,fun(Thetahat = x[i], Theta = Theta))
  }
  return(y)
}

1.2.3 置信区间的选取

SelDDa <- function(x){
  n = length(x);X = cumsum(x)
  d = n - sum(X > cumsum(x)[n]*0.025)
  u = sum(X > cumsum(x)[n]*0.975)
  x[c(1:d,(n-u+1):n)] = 0
  return(x)
}

1.2.4数据的构造

x <- seq(-1,1,0.001)
y_0 <- Gfun(f = fun, x = x, Theta = 0)
y_0.5 <- Gfun(f = fun, x = x, Theta = 0.5)
y_f0.5 <- Gfun(f = fun, x = x, Theta = -0.5)
p = length(x)
Data <- data.frame(x = rep(x,3),
    y = c(y_0,y_0.5,y_f0.5),
    z = factor(rep(c(0,0.5,-0.5),each = p)),
    y0 = c(rep(0,p),SelDDa(y_0.5),rep(0,p)),
    y1 = c(SelDDa(y_0),rep(0,p),rep(0,p)),
    y2 = c(rep(0,p),rep(0,p),SelDDa(y_f0.5)))

1.2.5 画图

P <- qplot(x, y, col = z, facets = z ~ ., data = Data,
xlab = expression(hat(theta)),ylab = expression(f[theta](hat(theta))),
show.legend = F)
P + geom_ribbon(aes(ymin = 0, ymax = y0),
show.legend = F,alpha=0.3,
fill = c("blue")) + geom_ribbon(aes(ymin = 0, ymax = y1),
show.legend = F,alpha=0.3,
fill = c("green")) + geom_ribbon(aes(ymin = 0, ymax = y2),
show.legend = F,alpha=0.3,
fill = c("red")) +
 geom_vline(xintercept=c(0.5), linetype="dotted", size = 1)

1.2图片

2、小结

通过本文你可以学到ggplot的一些画图操作：在图中添加希腊字母、公式、添加垂直虚线、用颜色填充面积、去除默认图例等。
在生成数据时，你可以学会构造画图所需数据类型、通过数据求分位数、使用R语言对一元函数进行求定积分等操作。
大家有更好的方法以及想法，欢迎评论区留言。