-
Notifications
You must be signed in to change notification settings - Fork 1
/
IQDfunctions.R
115 lines (90 loc) · 3.5 KB
/
IQDfunctions.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
##############################################################
# Calculate integrated quadratic distance (IQD) between the
# empirical cumulative distribution functions of two data
# vectors with the crps_sample function.
##############################################################
IQD_crps <- function(modelVec,obsVec){
IQD <- NA
# Vectors have values (not only NA's)
if(any(!is.na(modelVec)) & any(!is.na(obsVec))){
# Remove NA's
modelVec <- modelVec[!is.na(modelVec)]
obsVec <- obsVec[!is.na(obsVec)]
# IQD of ecdf difference
scoreModel <- 0
scoreObs <- 0
for(i in 1:length(obsVec)){
# method="edf" gives the empirical distribution function
scoreModel <- scoreModel+crps_sample(y=obsVec[i],dat=modelVec,method="edf")
scoreObs <- scoreObs+crps_sample(y=obsVec[i],dat=obsVec,method="edf")
}
IQD <- as.numeric((scoreModel-scoreObs)/length(obsVec))
}
IQD
}
##############################################################
# Calculate integrated quadratic distance (IQD) between the
# empirical cumulative distribution functions of two data
# vectors with the integrate function.
##############################################################
IQD_integrate<- function(modelVec,obsVec){
IQD <- NA
# Vectors have values (not only NA's)
if(any(!is.na(modelVec)) & any(!is.na(obsVec))){
# Remove NA's
modelVec <- modelVec[!is.na(modelVec)]
obsVec <- obsVec[!is.na(obsVec)]
# Limits of integration
lowerLimit <- range(c(modelVec,obsVec))[1]
upperLimit <- range(c(modelVec,obsVec))[2]
# IQD of ecdf difference
IQD <- integrate(Diff_ecdf,lower=lowerLimit,upper=upperLimit,subdivisions=5000,modelData=modelVec,obsData=obsVec)$value
}
IQD
}
#############################################
# Evaluate difference between two empirical
# cumulative distribution functions.
#############################################
Diff_ecdf <- function(x,modelData,obsData){
F <- ecdf(modelData)
G <- ecdf(obsData)
(F(x)-G(x))^2
}
##############################################################
# Calculate weighted integrated quadratic distance (IQD) between
# the empirical cumulative distribution functions of two data
# vectors with the integrate function.
##############################################################
wIQD_integrate<- function(modelVec,obsVec, wFct){
IQD <- NA
# Vectors have values (not only NA's)
if(any(!is.na(modelVec)) & any(!is.na(obsVec))){
# Remove NA's
modelVec <- modelVec[!is.na(modelVec)]
obsVec <- obsVec[!is.na(obsVec)]
# Limits of integration
lowerLimit <- range(c(modelVec,obsVec))[1]
upperLimit <- range(c(modelVec,obsVec))[2]
# IQD of ecdf difference
IQD <- integrate(wDiff_ecdf,lower=lowerLimit,upper=upperLimit,subdivisions=5000,modelData=modelVec,obsData=obsVec,wFct=wFct)$value
}
IQD
}
#############################################
# Evaluate weighted difference between two
# empirical cumulative distribution functions.
#############################################
wDiff_ecdf <- function(x,modelData,obsData,wFct){
F <- ecdf(modelData)
G <- ecdf(obsData)
(F(x)-G(x))^2 * wFct(x)
}
####################################################
# Print IQD values calculated with the crps_sample
# function and the integrate function.
####################################################
PrintIQD <- function(xVec,yVec){
print(paste("IQD_crps = ",round(IQD_crps(xVec,yVec),8),sep=""))
print(paste("IQD_integrate = ",round(IQD_integrate(xVec,yVec),8),sep=""))
}