library(ggplot2)
# length N={10, 100, 1000, 10000}
#######################--------------------define functions-----------------
# define a function of ordinal patterns
ordinal_pattern=function(x,dim){
# Generate ordinal numbers to assign.
ordinal_numbers=seq(0,(dim-1),by=1)
# Compute all possible permutations of the ordinal numbers.
# Maximum size of possible_pattern=dim!
possible_pattern=(combinat::permn(ordinal_numbers))
# Initialize result. Result is the output.
result=0
result[1:length(possible_pattern)]=0
# Loop for computation of ordinal pattern
for(i in 1:(length(x)-(dim-1))){
temp=x[i:(i+(dim-1))]
tempseq=seq(0,dim-1,by=1)
tempdata=data.frame(temp,tempseq)
tempdata=tempdata[order(temp),]
for(j in 1: length(possible_pattern)){
if (all(possible_pattern[[j]]==tempdata$tempseq)){
result[j]=result[j]+1
}
}
}
return(result)
}
# define a function of probability of each permutation
get_prob_dist <- function(perms) {
n <- length(perms)
probs <- NULL
for (i in 1:n) {
probs <- c(probs, (perms[i]/sum(perms)))
}
return(probs)
}
# permutation entropy function
# define a function of permutation entropy
permu_entropy=function(op){
# Compute maximum entropy. maximum entropy = log(dim!)
# or maximum entropy = log(length(ordinal_pattern))
entropy_max=log(length(op))
# Normalized permutation entropy
npe = entropy::entropy(op)
npe_norm=entropy::entropy(op)/entropy_max
res = data.frame(pe = npe, npe = npe_norm)
return(res)
}
#KL distance function
kld_function = function(x, y, n){
w = matrix(NA, nrow = n, ncol = 1)
w[1,] = x[1,]
for (i in 2:length(x)){
w[i,] = w[i-1,] + x[i,]
}
### or using cumsum() function to calculate W
permus_w = ordinal_pattern(w, 3)
permus_prob_w = get_prob_dist(permus_w)
#random walk z = sum of y_i
z = matrix(NA, nrow = n, ncol = 1)
z[1,] = y[1,]
for (i in 2:length(y)){
z[i,] = z[i-1,] + y[i,]
}
permus_z = ordinal_pattern(z, 3)
permus_prob_z = get_prob_dist(permus_z)
Dist = NULL
for (j in 1:6){
Dist[j] = permus_prob_w[j]*log((permus_prob_w[j])/(permus_prob_z[j]))
}
total_dist = sum(Dist)
return(total_dist)
}
##################-------------Generation-----------------
dim = 3
#generating random variables of steps
set.seed(234)
n1=10
x1 = matrix(runif(n1, min = -0.3, max = 0.7))
px1 = x1[x1>0,]
p_x1 = length(px1)/n1
p_x1
## [1] 0.7
y1 = matrix(runif(n1, min = -0.5, max = 0.5))
n2=100
x2 = matrix(runif(n2, min = -0.3, max = 0.7))
px2 = x2[x2>0,]
p_x2 = length(px2)/n2
p_x2
## [1] 0.73
y2 = matrix(runif(n2, min = -0.5, max = 0.5))
n3=1000
x3 = matrix(runif(n3, min = -0.3, max = 0.7))
px3 = x3[x3>0,]
p_x3 = length(px3)/n3
p_x3
## [1] 0.674
y3 = matrix(runif(n3, min = -0.5, max = 0.5))
n4=10000
x4 = matrix(runif(n4, min = -0.3, max = 0.7))
px4 = x4[x4>0,]
p_x4 = length(px4)/n4
p_x4
## [1] 0.7003
y4 = matrix(runif(n4, min = -0.5, max = 0.5))
##compute possible # of ordinal patterns
op_x1 = ordinal_pattern(x1, dim)
op_x2 = ordinal_pattern(x2, dim)
op_x3 = ordinal_pattern(x3, dim)
op_x4 = ordinal_pattern(x4, dim)
op_y1 = ordinal_pattern(y1, dim)
op_y2 = ordinal_pattern(y2, dim)
op_y3 = ordinal_pattern(y3, dim)
op_y4 = ordinal_pattern(y4, dim)
###compute permutation entropy
pe_x1 = permu_entropy(op_x1)
pe_x2 = permu_entropy(op_x2)
pe_x3 = permu_entropy(op_x3)
pe_x4 = permu_entropy(op_x4)
pe_y1 = permu_entropy(op_y1)
pe_y2 = permu_entropy(op_y2)
pe_y3 = permu_entropy(op_y3)
pe_y4 = permu_entropy(op_y4)
resl_x = rbind(pe_x1, pe_x2, pe_x3, pe_x4)
resl_x
## pe npe
## 1 1.320888 0.7372018
## 2 1.780233 0.9935667
## 3 1.790323 0.9991982
## 4 1.791598 0.9999101
resl_y = rbind(pe_y1, pe_y2, pe_y3, pe_y4)
resl_y
## pe npe
## 1 1.494175 0.8339150
## 2 1.754533 0.9792236
## 3 1.786027 0.9968006
## 4 1.791373 0.9997842
##compute KLD
kl1 = kld_function(y1, x1, n1)
kl1
## [1] NaN
kl2 = kld_function(y2, x2, n2)
kl2
## [1] 0.472363
kl3 = kld_function(y3, x3, n3)
kl3
## [1] 0.2495222
kl4 = kld_function(y4, x4, n4)
kl4
## [1] 0.2629817
resl_kl = rbind(kl1, kl2, kl3, kl4)
resl_kl
## [,1]
## kl1 NaN
## kl2 0.4723630
## kl3 0.2495222
## kl4 0.2629817