-
Notifications
You must be signed in to change notification settings - Fork 0
/
stat.R
115 lines (103 loc) · 3.16 KB
/
stat.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
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
#Goodness of Fit Measures Used
#KS, MSE of distirbution, Log-likelihood
#
#Ari Purwanto Sarwo Prasojo & Puguh Prasetyoputra (2019)
#Research Center for Population, Indonesian Institute of Sciences
#________________________________________________________________
#Calculate KS distance ----
ks.stat <- function(x,distr,...){
#________________________________________________________________
#ks.stat
#Calculating ks distance for lognormal and weibull distribution
#
#/Usage
#ks.stat(x,distr,...)
#/Arguments
#x : univariate continuous random sample
#distr : choices of distribution ("lognormal" or "weibull")
#... : additional arguments (parameters of distribution)
#
#/Value
#ks : ks distance
#________________________________________________________________
#assign distr pars
pars <- list(...)
for(par in names(pars)){
assign(par, pars[[par]])
}
#sort data
x <- sort(x)
n <- length(x)
#calculating empirical c.d.f
ecdfx <- ecdf(x)
Sx <- ecdfx(x)
Sx1 <- c(0, Sx[-n])
#calculating theoretical c.d.f
if(distr == "lognormal"){
Fx <- plnorm(x,meanlog,sdlog)
}else if(distr == "weibull"){
Fx <- pweibull(x,shape,scale)
}
#max distance between emprical & theoretical c.d.f
ks <- max(pmax(abs(Fx-Sx1), abs(Sx-Fx)))
return(ks)
}
#Calculate MSE of distribution fitting ----
mse.stat <- function(x,distr,...){
#________________________________________________________________
#mse.stat
#Calculating mean square error (mse) of distribution fitting
#
#/Usage
#mse.stat(x,distr,...)
#/Arguments
#x : univariate continuous random sample
#distr : choices of distribution ("lognormal" or "weibull")
#... : additional arguments (parameters of distribution)
#
#/Value
#mse : mse of distribution fitting
#________________________________________________________________
#assign distr pars
pars <- list(...)
for(par in names(pars)){
assign(par, pars[[par]])
}
#sort data
x <- sort(x)
n <- length(x)
#calculating empirical c.d.f
ecdfx <- ecdf(x)
Sx <- ecdfx(x)
#calculating theoretical c.d.f
if(distr == "lognormal"){
Fx <- plnorm(x,meanlog,sdlog)
}else if(distr == "weibull"){
Fx <- pweibull(x,shape,scale)
}
mse <- sum((Fx-Sx)^2)/n
return(mse)
}
#Calculate Log-likelihood ----
loglik <- function(x,distr,theta){
#________________________________________________________________
#loglik
#Calculating log-likelihood of lognormal or weibull distribution
#
#/Usage
#loglik(x,distr,theta)
#/Arguments
#x : univariate continuous random sample
#distr : choices of distribution ("lognormal" or "weibull")
#theta : parameters of distirbution (meanlog,sdlog) or (shape,scale)
#
#/Value
#llh : log-likelihood value
#________________________________________________________________
if(distr == "lognormal"){
llh <- sum(dlnorm(x,meanlog = theta[1],sdlog = theta[2],log = TRUE))
}else if(distr == "weibull"){
llh <- sum(dweibull(x,shape = theta[1],scale = theta[2],log = TRUE))
}
return(llh)
}