distributions.R (10318B)
1 dweibull2 <- function (x, a1=1, b1=1, c1=1, a2=1, b2=1, c2=1, g=0) { 2 #ifelse(x<g, 3 # a1* abs(c1*b1)*(abs(b1*(x-g))**(c1-1))* exp(-(abs(b1*(x-g))**c1)), 4 # a1* abs(c2*b2)*(abs(b2*(x-g))**(c2-1))* exp(-(abs(b2*(x-g))**c2))) 5 ifelse(x<g, 6 a1*abs(b1*c1)*exp(-(abs(b1*(abs(x-g)-abs((c1-1)/c1)**(1/c1)))**c1)), 7 a2*abs(b2*c2)*exp(-(abs(b2*(abs(x-g)-abs((c2-1)/c2)**(1/c2)))**c2))) 8 } 9 10 dweibull4 <- function (x, a1=1, b1=1, c1=1, a2=1, b2=1, c2=1, g1=0, 11 a3=1, b3=1, c3=1, a4=1, b4=1, c4=1, g2=0) { 12 ifelse(x<0, 13 dweibull2(x,a1,b1,c1,a2,b2,c2,g1), 14 dweibull2(x,a3,b3,c3,a4,b4,c4,g2)) 15 } 16 17 dvonmises <- function (x, x_mean=0, a=1, b=1, c=1) { 18 a * exp(b*cos(c*(x-x_mean))) / (2*pi*besselI(b,0)) 19 } 20 21 dgumbel <- function (x, location=0, scale=1) { 22 z <- (x-location)/scale 23 (1/scale) * exp(-(z + exp(-z))) 24 } 25 26 # https://www.tutorialspoint.com/r/r_mean_median_mode.htm 27 # https://stackoverflow.com/questions/2547402/how-to-find-the-statistical-mode 28 getmode <- function(v) { 29 #uniqv <- unique(v) 30 #uniqv[which.max(tabulate(match(v, uniqv)))] 31 as.numeric(names(sort(-table(v)))[1]) 32 } 33 34 # https://stackoverflow.com/questions/27418461/calculate-the-modes-in-a-multimodal-distribution-in-r 35 find_peaks_indices <- function (x) { 36 modes <- c() 37 for (i in 2:(length(x)-1)){ 38 if ((x[i] > x[i-1]) & (x[i] > x[i+1])) { 39 modes <- c(modes,i) 40 } 41 } 42 modes 43 } 44 45 find_real_peaks_indices <- function (x, max_modes=2) { 46 n <- 1024 47 #x_density <- density(x,n=n,adjust=0.125) 48 #y <- x_density$y 49 y <- x 50 modes <- find_peaks_indices(y) 51 # find the first N highest peaks 52 modes_xy <- data.frame(x=modes,y=y[modes]) 53 modes_xy <- modes_xy[order(modes_xy$y, decreasing=TRUE),] 54 modes_xy <- modes_xy[c(1:max_modes), ] 55 # scale indices to [0,1] 56 (modes_xy$x-1)/(n-1) 57 } 58 59 find_max_x <- function (x, y) { 60 n <- length(y)-1 61 indices <- c(0:n) 62 index <- indices[y == max(y)]/n 63 x[floor(index*(length(x)-1)+1)] 64 } 65 66 normalize_hist <- function(pdf_x) { 67 #data.frame(density=pdf_x$density/max(abs(pdf_x$density)),x=pdf_x$mids) 68 data.frame(density=pdf_x$density,x=pdf_x$mids) 69 } 70 71 fit_velocity_pdf <- function (xx, x, density, sign=-1, plot=TRUE, axis="x", dist="weibull2") { 72 a0 <- 0 73 b0 <- 0 74 c0 <- 1 75 c1 <- 3 76 indices <- find_real_peaks_indices(x, max_modes=1) 77 if (FALSE) { 78 pdf <- normalize_hist(hist(x, plot=FALSE, breaks=1000)) 79 print(indices) 80 indices = floor(indices * length(pdf$x)) 81 x_modes <- pdf$x[indices] 82 print(x_modes) 83 model <- nls(pdf$density ~ dweibull4(pdf$x,a1,b1,c1,a2,b2,c2,x_modes[[1]],a3,b3,c3,a4,b4,c4,x_modes[[2]]), 84 data=pdf, 85 start=list(a1=1,b1=1,c1=1,a2=1,b2=1,c2=1, 86 a3=1,b3=1,c3=1,a4=1,b4=1,c4=1), 87 control=list(warnOnly=TRUE), 88 algorithm="port", 89 lower=c(a0,b0,c0,a0,b0,c0,a0,b0,c0,a0,b0,c0), 90 upper=c(1000,1000,3,1000,2000,3,1000,1000,3,1000,2000,3)) 91 c = summary(model)$coefficients 92 if (plot) { 93 pdf_x_est <- dweibull4(pdf$x,c[[1]],c[[2]],c[[3]],c[[4]],c[[5]],c[[6]],x_modes[[1]], 94 c[[7]],c[[8]],c[[9]],c[[10]],c[[11]],c[[12]],x_modes[[2]]) 95 plot(pdf$x, pdf$density, xlab=paste('Velocity', axis), ylab='PDF') 96 lines(pdf$x, pdf_x_est, col='red') 97 } 98 c 99 } else { 100 #x_mode <- getmode(x) 101 #pdf_x <- normalize_hist(hist(x, plot=FALSE, breaks=100)) 102 pdf_x <- data.frame(x=xx,density=x) 103 c <- c() 104 pdf_x_set <- c() 105 if (dist == "weibull2") { 106 x_mode <- find_max_x(pdf_x$x,pdf_x$density) 107 #indices = floor(indices * length(pdf_x$x)) + 1 108 #x_mode <- pdf_x$x[indices][[1]] 109 c <- tryCatch( 110 { 111 model <- nls(pdf_x$density ~ dweibull2(pdf_x$x,a,b,c,d,e,f,g=x_mode), 112 data=pdf_x, 113 start=list(a=1,b=1,c=1,d=1,e=1,f=1), 114 control=list(warnOnly=TRUE), 115 algorithm="port", 116 lower=c(a0,b0,c0,a0,b0,c0), upper=c(1000,1000,c1,1000,2000,c1)) 117 summary(model)$coefficients 118 }, 119 error=function(cond) { 120 message(cond) 121 c(1,1,1,1,1,1) 122 }) 123 pdf_x_est <- dweibull2(pdf_x$x,c[[1]],c[[2]],c[[3]],c[[4]],c[[5]],c[[6]],x_mode) 124 } else if (dist == "vonmises") { 125 c <- tryCatch({ 126 model <- nls(density ~ dvonmises(2*pi*pdf_x$x/(max(pdf_x$x)),0,a,b,c), 127 data=pdf_x, 128 start=list(a=1,b=1,c=1), 129 control=list(warnOnly=TRUE), 130 algorithm="port", 131 lower=c(0,1,0), upper=c(1000,1000,pi)) 132 c <- summary(model)$coefficients 133 },error=function(cond) { 134 message(cond) 135 c(1,1,1) 136 }) 137 pdf_x_est <- dvonmises(2*pi*pdf_x$x/max(pdf_x$x),0,c[[1]],c[[2]],c[[3]]) 138 } else if (dist == "gumbel") { 139 c <- tryCatch( 140 { 141 model <- nls(pdf_x$density ~ c*dgumbel(pdf_x$x,location=a,scale=b), 142 data=pdf_x, 143 start=list(a=1,b=1,c=1), 144 control=list(warnOnly=TRUE), 145 algorithm="port", 146 lower=c(-1000,0.1,c=1), upper=c(1000,1000,1000)) 147 c = summary(model)$coefficients 148 }, 149 error=function(cond) { 150 message(cond) 151 c <- c(1,1,1,1) 152 }) 153 pdf_x_est <- c[[3]]*dgumbel(pdf_x$x,location=c[[1]],scale=c[[2]]) 154 } else if (dist == "norm") { 155 fit <- fitdist(pdf_x$x, "norm") 156 c <- coef(fit) 157 pdf_x_est <- dnorm(pdf_x$x,mean=c[[1]],sd=c[[2]]) 158 } else if (dist == "rayleigh") { 159 c <- tryCatch( 160 { 161 model <- nls(pdf_x$density ~ b*dweibull(sign*pdf_x$x,scale=a,shape=2), 162 data=pdf_x, 163 start=list(a=1,b=1), 164 control=list(warnOnly=TRUE), 165 algorithm="port", 166 lower=c(0.1,0), upper=c(1000,1000)) 167 c = summary(model)$coefficients 168 }, 169 error=function(cond) { 170 message(cond) 171 c <- c(1,1,1,1) 172 }) 173 pdf_x_est <- c[[2]]*dweibull(sign*pdf_x$x,scale=c[[1]],shape=2) 174 } else { 175 x_mode <- find_max_x(pdf_x$x,pdf_x$density) 176 c <- tryCatch( 177 { 178 model <- nls(pdf_x$density ~ c*dweibull(sign*pdf_x$x,scale=a,shape=b), 179 data=pdf_x, 180 start=list(a=1,b=1,c=1), 181 control=list(warnOnly=TRUE), 182 algorithm="port", 183 lower=c(0.1,0.1,0), 184 upper=c(1000,1000,1000)) 185 c = summary(model)$coefficients 186 }, 187 error=function(cond) { 188 message(cond) 189 c <- c(1,1,1,1) 190 }) 191 pdf_x_est <- c[[3]]*dweibull(sign*pdf_x$x,scale=c[[1]],shape=c[[2]]) 192 } 193 #if (plot) { 194 # plot(sign*pdf_x$x, pdf_x$density, xlab=paste('Velocity', axis), ylab='PDF') 195 # lines(density(sign*x), col='blue') 196 # lines(sign*pdf_x$x, pdf_x_est, col='red') 197 #} 198 #if (sign < 0) { 199 # print(pdf_x$x) 200 # pdf_x$x = -rev(pdf_x$x) 201 # pdf_x$density = rev(pdf_x$density) 202 # pdf_x_est = rev(pdf_x_est) 203 #} 204 list(coefficients=c, 205 actual_density=pdf_x$density, 206 estimated_density=pdf_x_est, 207 x=pdf_x$x, 208 residual=max(abs(pdf_x$density-pdf_x_est))) 209 } 210 } 211 212 fit_velocity_acf <- function (velocity, plot=TRUE, axis="x", timestamp=0, write=TRUE, 213 debug=FALSE) { 214 acf_x <- acf(velocity, type="covariance", plot=FALSE) 215 acf_x <- data.frame(acf=acf_x$acf, lag=acf_x$lag) 216 if (debug) { 217 plot(acf_x$lag, acf_x$acf, xlab='Lag', ylab=paste('DEBUG ACF', axis)) 218 } 219 v_mean <- mean(velocity) 220 c <- tryCatch({ 221 model <- nls(acf ~ a*exp(-abs(b*lag)**c), 222 data=acf_x, 223 start=list(a=1,b=1,c=1), 224 control=list(warnOnly=TRUE), 225 algorithm="port", 226 lower=c(0,0,1e-6), upper=c(1000,1000,10)) 227 c <- summary(model)$coefficients 228 },error=function(cond) { 229 message(cond) 230 c(1,1,1) 231 }) 232 acf_est <- c[[1]]*exp(-abs(c[[2]]*acf_x$lag)**c[[3]]) 233 if (plot) { 234 plot(acf_x$lag, acf_x$acf, xlab='Lag', ylab=paste('ACF', axis)) 235 lines(acf_x$lag, acf_est, col='red') 236 } 237 if (write) { 238 path <- file.path('build', 'acf', axis) 239 make_directory(path) 240 filename <- file.path(path, sprintf("%d", timestamp)) 241 write.table(data.frame(lag=acf_x$lag,actual=acf_x$acf,estimated=acf_est), 242 filename, row.names=FALSE, quote=FALSE) 243 } 244 list(coefficients=c, 245 rmse=rmse(acf_est,acf_x$acf), 246 lag=acf_x$lag, 247 actual=acf_x$acf, 248 estimated=acf_est) 249 } 250 251 fit_direction <- function (theta, density, plot=TRUE) { 252 pdf_x <- data.frame(x=theta,density=density); 253 #print(pdf_x) 254 scale <- max(abs(theta)) 255 c <- tryCatch({ 256 model <- nls(density ~ dvonmises(2*pi*x/scale-m,0,a,b,c), 257 data=pdf_x, 258 start=list(a=1,b=1,c=1,m=0), 259 control=list(warnOnly=TRUE), 260 algorithm="port", 261 lower=c(0,1,0,-10), 262 upper=c(1000,1000,1000,10)) 263 c <- summary(model)$coefficients 264 },error=function(cond) { 265 message(cond) 266 c(1,1,1) 267 }) 268 pdf_x_est <- dvonmises(2*pi*pdf_x$x/scale,c[[4]],c[[1]],c[[2]],c[[3]]) 269 if (plot) { 270 plot(pdf_x$x, pdf_x$density, xlab='Direction', ylab='PDF') 271 lines(pdf_x$x, pdf_x_est, col='red') 272 } 273 list(coefficients=c, 274 actual_density=pdf_x$density, 275 estimated_density=pdf_x_est, 276 x=pdf_x$x, 277 residual=max(abs(pdf_x$density-pdf_x_est))) 278 } 279