估计概率密度函数
大概查看服从哪些基本的分布
library(fitdistrplus)
library(logspline)
descdist(data, discrete = FALSE)
一、参数法
1、正态分布
检验:
shapiro.test(data)
方法1:
fit.normal <- fitdist(data_tmp$VE,"norm")
summary(fit.normal)
方法2:
fit.normal_b<-Bnormal(data_tmp$VE, priors=list(muMean=13.651040, muSD=3.998584))
summary(fit.normal_b)
2、伽马分布
检验:
library(goft)
gamma_test(data)
方法1:
fit.gamma <- fitdist(data_tmp$VA, distr = "gamma", method = "mle")
summary(fit.gamma)
方法2:
mcmc <- function(y, iter, sigma.rw, lambda.alpha, lambda.beta, nu.alpha, nu.beta) {
# find n from the data (ie the number of observations)
n <- length(y)
# compute the sum of the observations and the sum of the log of the observations
sum.y <- sum(y)
sum.log.y <- sum(log(y))
# first create a table which will store the samples; first column will be alpha, second will be beta
out <- matrix(NA, nrow=iter, ncol=2)
# initial values
alpha.cur <- 1
beta.cur <- 1
out[1,] <- c(alpha.cur, beta.cur)
# mcmc loop starts here
for (i in 2:iter) {
###############
# update alpha (assume beta is fixed)
###############
# propose a new value for alpha
alpha.can <- rnorm(1, alpha.cur, sigma.rw)
# if it is negative reject straight away else compute the M-H ratio
if (alpha.can > 0) {
# evaluate the loglikelihood at the current values of alpha.
loglik.cur <- n*alpha.cur*log(beta.cur) - n*lgamma(alpha.cur) + (alpha.cur-1)*sum.log.y
# compute the log-likelihood at the candidate value of alpha
loglik.can <- n*alpha.can*log(beta.cur) - n*lgamma(alpha.can) + (alpha.can-1)*sum.log.y
# log prior densities
log.prior.alpha.cur <- log.gamma.density(alpha.cur, lambda.alpha, nu.alpha)
log.prior.alpha.can <- log.gamma.density(alpha.can, lambda.alpha, nu.alpha)
logpi.cur <- loglik.cur + log.prior.alpha.cur
logpi.can <- loglik.can + log.prior.alpha.can
# M-H ratio
# draw from a U(0,1)
u <- runif(1)
if (log(u) < logpi.can - logpi.cur) {
alpha.cur <- alpha.can
}
}
###############
# update beta
###############
beta.cur <- rgamma(1, alpha.cur*n + lambda.beta, rate = sum.y + nu.beta)
# store the samples
out[i,] <- c(alpha.cur, beta.cur)
}
return(out)
}
res <- mcmc(y = data_tmp$VA, iter = 10000, sigma.rw = 1, lambda.alpha = 1, lambda.beta = 1, nu.alpha = 1e-4, nu.beta = 1e-4)
res <- res[-c(1:1000),]
mean(res[,1])
mean(res[,2])
sd(res[,1])
sd(res[,2])
3、贝塔分布
方法1:
fit.beta <- fitdist(data_tmp$VA, distr = "beta", method="mle")
summary(fit.beta)
方法2:
library("LearnBayes")
library(Lahman)
findBeta <- function(quantile1,quantile2,quantile3)
{
# find the quantiles specified by quantile1 and quantile2 and quantile3
quantile1_p <- quantile1[[1]]; quantile1_q <- quantile1[[2]]
quantile2_p <- quantile2[[1]]; quantile2_q <- quantile2[[2]]
quantile3_p <- quantile3[[1]]; quantile3_q <- quantile3[[2]]
# find the beta prior using quantile1 and quantile2
priorA <- beta.select(quantile1,quantile2)
priorA_a <- priorA[1]; priorA_b <- priorA[2]
# find the beta prior using quantile1 and quantile3
priorB <- beta.select(quantile1,quantile3)
priorB_a <- priorB[1]; priorB_b <- priorB[2]
# find the best possible beta prior
diff_a <- abs(priorA_a - priorB_a); diff_b <- abs(priorB_b - priorB_b)
step_a <- diff_a / 100; step_b <- diff_b / 100
if (priorA_a < priorB_a) { start_a <- priorA_a; end_a <- priorB_a }
else { start_a <- priorB_a; end_a <- priorA_a }
if (priorA_b < priorB_b) { start_b <- priorA_b; end_b <- priorB_b }
else { start_b <- priorB_b; end_b <- priorA_b }
steps_a <- seq(from=start_a, to=end_a, length.out=1000)
steps_b <- seq(from=start_b, to=end_b, length.out=1000)
max_error <- 10000000000000000000
best_a <- 0; best_b <- 0
for (a in steps_a)
{
for (b in steps_b)
{
# priorC is beta(a,b)
# find the quantile1_q, quantile2_q, quantile3_q quantiles of priorC:
priorC_q1 <- qbeta(c(quantile1_p), a, b)
priorC_q2 <- qbeta(c(quantile2_p), a, b)
priorC_q3 <- qbeta(c(quantile3_p), a, b)
priorC_error <- abs(priorC_q1-quantile1_q) +
abs(priorC_q2-quantile2_q) +
abs(priorC_q3-quantile3_q)
if (priorC_error < max_error)
{
max_error <- priorC_error; best_a <- a; best_b <- b
}
}
}
print(paste("The best beta prior has a=",best_a,"b=",best_b))
}
quantile1 <- list(p=0.5, x=median(data_tmp$VA))
x1<-quantile(data_tmp$VA,probs = 0.8)
x2<-quantile(data_tmp$VA,probs = 0.2)
quantile2 <- list(p=0.8,x=x1)
quantile3 <- list(p=0.2,x=x2)
findBeta(quantile1,quantile2,quantile3)
二、非参数法
1、直方图法
a1<-ggplot(data_tmp, aes(x=VA)) +
geom_histogram(binwidth=0.005,aes(y=..density..),colour="black",fill="white")+
geom_density(alpha=.3,col='red') +
theme_classic()+
annotate('text', x = 0.8, y = 8,
label = "Delta ==0.005 ",parse = TRUE,size=4)
a2<-ggplot(data_tmp, aes(x=VA)) +
geom_histogram(binwidth=0.02,aes(y=..density..),colour="black",fill="white")+
geom_density(alpha=.3,col='red') +
theme_classic()+
annotate('text', x = 0.8, y = 3.8,
label = "Delta ==0.02 ",parse = TRUE,size=4)
a3<-ggplot(data_tmp, aes(x=VA)) +
geom_histogram(binwidth=0.15,aes(y=..density..),colour="black",fill="white")+
geom_density(alpha=.3,col='red') +
theme_classic()+
annotate('text', x = 0.8, y = 3.8,
label = "Delta ==0.15 ",parse = TRUE,size=4)
abc<-grid.arrange(a1,a2,a3,
ncol=1,nrow=3,
as.table=TRUE)
只画下左边半个图
2、Kernel
png(file = "kernel.png",width = 1080, height = 1080, units = "px")
par(mfrow=c(1,2))
plot(density(data_tmp$VA,bw = 0.005, kernel = "gaussian"),col="red",main="Kernel Estimation",xlab="VA")
lines(density(data_tmp$VA,bw = 0.01, kernel = "gaussian"),col="blue")
lines(density(data_tmp$VA,bw = 0.05, kernel = "gaussian"),col="green")
rug(data_tmp$VA)
legend("topright", legend=c("bw=0.005", "bw=0.01","bw=0.05"),
col=c("red", "blue","green"), lty=1:1, cex=0.7)
只画下左边半个图
3、KNN
library( basicTrendline)
dk_a=function(A,x,k){
na=nrow(A)
or=1:na
dis=NULL
for(i in 1:na)
{dis=c(dis,(abs(x-A[i,279])))}
disnew=sort(dis)
dk=disnew[k]
return(dk)
}
knear_a=function(A,x,k){
knear=c()
for(i in 1:nrow(A)){
knear[i]=(k-1)/(2*nrow(A)*dk_a(A,x[i],k))}
knear}
set.seed(123)
data_k<-data_tmp[sample(1:nrow(data_tmp),5000),]
VA<-sort(data_k$VA)
k1<-knear_a(data_k,VA,100)
k2<-knear_a(data_k,VA,500)
k3<-knear_a(data_k,VA,800)
k4<-knear_a(data_k,VA,1000)
par(mfrow=c(2,2))
plot(VA,k1,type="l",xlab="VA",ylab="Density",main='k=100')
plot(VA,k2,type="l",xlab="VA",ylab="Density",main='k=500')
plot(VA,k3,type="l",xlab="VA",ylab="Density",main='k=800')
plot(VA,k4,type="l",xlab="VA",ylab="Density",main='k=1000')