-
Notifications
You must be signed in to change notification settings - Fork 0
/
genFunctions.R
95 lines (80 loc) · 3.34 KB
/
genFunctions.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
genFits <- function(x, type = "norm", method="mge", gof="CvM"){
ini.par = mledist(x[x >= 0], "weibull")
fit <- switch(type,
"norm" = fitdist(x, "norm", method = method, gof = gof, start = list(mean=mean(x), sd=sd(x))),
"logis" = fitdist(x, "logis", method = method, gof = gof),
"extreme" = fextreme(x, 'extreme', method=method, gof=gof, start=list(mean=mean(x), sd=sd(x))),
"gumbel" = fgumbel(x, 'gumbel', method=method, gof=gof, start=list(mean=mode(x), sd=0.78*sd(x))),
"weibull" = fitdist(x, "weibull", method = method, gof = gof, start=list(shape=ini.par$estimate[[1]], scale=ini.par$estimate[[2]])),
return(fit)
)}
dgumbel <- function(x, mean, sd) 1/sd*exp((mean-x)/sd)*exp(-exp((mean-x)/sd)) #f
pgumbel <- function(q, mean, sd) exp(-exp((mean-q)/sd)) #F
qgumbel <- function(p, mean, sd) mean-sd*log(-log(p))
rgumble <- function(n) log(-log(runif(n)))
fgumbel <- function(x, type = "gumbel" , method="mge", gof="CvM", start=list(mean=mode(x), sd=0.78*sd(x))) {
fitdist(x, 'gumbel', method=method, gof=gof, start=list(mean=mode(x), sd=0.78*sd(x)))
}
dextreme <- function(x, mean, sd) 1/sd*exp((-mean+x)/sd)*exp(-exp((-mean+x)/sd))
pextreme <- function(q, mean, sd) 1-exp(-exp((-mean+q)/sd))
qextreme <- function(p, mean, sd) mean+sd*log(-log(1-p))
rextreme <- function(n) log(-log(1-runif(n)))
fextreme <- function(x, type = 'extreme', method="mge", gof="CvM", start=list(mean=mean(x), sd=sd(x))) {
fitdist(x, 'extreme', method=method, gof=gof, start=list(mean=mean(x), sd=sd(x)))
}
# getParams <- function(fit, bootmethod = "nonparam", niter = 5000){
# b4 <-bootdist(fit, bootmethod = bootmethod, niter=niter)
# sumb4 <- summary(b4)
# par41 <- sumb4$fitpart$estimate[1]
# par42 <- sumb4$fitpart$estimate[2]
# return(list(par1 = par41, par2 = par42))
# }
getParams <- function(fit.nor){
sumb4 <- summary(fit.nor)
par41 <- sumb4$estimate[1]
par42 <- sumb4$estimate[2]
return(list(par1 = par41, par2 = par42))
}
mode <- function(v) {
uniqv <- unique(v)
uniqv[which.max(tabulate(match(v, uniqv)))]
}
SSD <- function(x, p.sim) {
inv.ecdf <- function(x){return(ecdf(x)(x))}
p.data <- inv.ecdf(x)
ssd <- sum((p.sim - p.data)^2)
return(ssd)
}
MSE <- function(x, p.sim) {
inv.ecdf <- function(x){return(ecdf(x)(x))}
p.data <- inv.ecdf(x)
ssd <- sum((p.sim - p.data)^2)
mse <- ssd/(length(x)-2)
return(mse)
}
hc <- function(fit, bootfit, p) {
hc <- quantile(fit, probs = p)
hcq <- quantile(bootfit, probs = p)
return(list(hc=hc, hcq=hcq))
}
HC <- function(model.name, scale, hc.list, p) {
if(scale == 2) {
HCp <- as.data.frame(sapply(hc.list, function(i) c(10^i$hc$quantile$p, 10^i$hcq$quantCI$p)))
}else{
HCp <- as.data.frame(sapply(hc.list, function(i) c(i$hc$quantile$p, i$hcq$quantCI$p)))
}
colnames(HCp) <- model.name
rownames(HCp) <- list(paste0("HC", p*100), paste0("HC", p*100, "L"), paste0("HC", p*100, "H"))
HCp <- format(HCp, digits=3, nsmall = 3, justify="left")
return(HCp)
}
fx <- function(scale, x, y) ifelse(scale == 2, return(x), return(y))
CI.boot <- function(scale, b1, CI.level) {
qboot <- quantile(b1, probs = seq(0.001,0.999,0.001), CI.type = "two.sided", CI.level = CI.level)
cis <- qboot$quantCI
cis <- fx(scale, 10^cis, cis)
Pr <- qboot$probs
rownames(cis) <- c('lwr' ,'upr')
CI.Pr <- data.frame(cbind(Pr, t(cis)))
return(CI.Pr)
}