anal.R (22974B)
1 library(fitdistrplus) 2 library(plotrix) 3 library(sn) 4 library(doParallel) 5 library(foreach) 6 library(gplots) # hist2d 7 registerDoParallel(cores=detectCores()) 8 9 source("R/common.R") 10 source("R/distributions.R") 11 12 # https://stackoverflow.com/questions/4954507/calculate-the-area-under-a-curve 13 area_under_curve <- function (x,y) { 14 sum(diff(x) * (head(y,-1)+tail(y,-1)))/2 15 } 16 17 fit_velocity_pdf_bilateral <- function (velocity, axis="x", timestamp=0, dist="weibull") { 18 velocity_hist <- normalize_hist(hist(velocity, plot=FALSE, breaks=100)) 19 x = velocity_hist$x; 20 v = velocity_hist$density; 21 d = density(velocity,adjust=0.5) 22 dist_left <- fit_velocity_pdf(x[x<=0], v[x<=0], density=d$y[d$x<=0], 23 sign=-1, plot=FALSE, axis=axis, dist=dist) 24 dist_right <- fit_velocity_pdf(x[x>=0], v[x>=0], density=d$y[d$x>=0], 25 sign=1, plot=FALSE, axis=axis, dist=dist) 26 #x <- c(dist_left$x, dist_right$x) 27 actual_density <- c(dist_left$actual_density, dist_right$actual_density) 28 estimated_density <- c(dist_left$estimated_density, dist_right$estimated_density) 29 plot(x, actual_density, xlab=paste('Velocity', axis), ylab='PDF') 30 #hist(velocity, xlab=paste('Velocity', axis), ylab='PDF', freq=FALSE, breaks=100) 31 lines(d, col='green', lwd=2) 32 #lines(x, estimated_density, col='red', lwd=2) 33 lines(dist_left$x, dist_left$estimated_density, col='red', lwd=2) 34 lines(dist_right$x, dist_right$estimated_density, col='blue', lwd=2) 35 #printf("RMSE %f\n", rmse(estimated_density, actual_density)) 36 path <- file.path('build', 'velocity', axis, dist) 37 make_directory(path) 38 filename <- file.path(path, sprintf("%d", timestamp)) 39 write.table(data.frame(x=x,actual=actual_density,estimated=estimated_density), 40 filename, row.names=FALSE, quote=FALSE) 41 write("\n",file=filename,append=TRUE) 42 write.table(data.frame(x=d$x,density=d$y), filename, 43 row.names=FALSE, quote=FALSE, append=TRUE) 44 area_left <- area_under_curve(dist_left$x, dist_left$actual_density) 45 area_right <- area_under_curve(dist_right$x, dist_right$actual_density) 46 list(coefficients=dist_left$coefficients, 47 rmse=rmse(estimated_density, actual_density), 48 area_min=min(area_left,area_right), 49 area_max=max(area_left,area_right)) 50 } 51 52 fit_velocity_acf_bilateral <- function (velocity, axis="x", timestamp=0, plot=TRUE) { 53 v = velocity 54 model_1 <- fit_velocity_acf(v[v<=0], axis=axis, timestamp=timestamp, plot=FALSE, write=FALSE); 55 model_2 <- fit_velocity_acf(v[v>=0], axis=axis, timestamp=timestamp, plot=FALSE, write=FALSE); 56 #x <- c(dist_left$x, dist_right$x) 57 lag <- c(-model_1$lag, model_2$lag) 58 estimated <- c(model_1$estimated, model_2$estimated) 59 actual <- c(model_1$actual, model_2$actual) 60 plot(lag, actual, xlab='Lag', ylab='ACF') 61 lines(-model_1$lag, model_1$estimated, col='red', lwd=2) 62 lines(model_2$lag, model_2$estimated, col='blue', lwd=2) 63 path <- file.path('build', 'acf', axis) 64 make_directory(path) 65 filename <- file.path(path, sprintf("%d", timestamp)) 66 #write.table(data.frame(lag=lag,actual=actual,estimated=estimated), 67 # filename, row.names=FALSE, quote=FALSE) 68 write.table(data.frame(lag=model_1$lag,actual=model_1$actual,estimated=model_1$estimated), 69 filename, row.names=FALSE, quote=FALSE) 70 write("\n",file=filename,append=TRUE) 71 write.table(data.frame(lag=model_2$lag,actual=model_2$actual,estimated=model_2$estimated), 72 filename, row.names=FALSE, quote=FALSE, append=TRUE) 73 mean_1 <- abs(mean(v[v<=0])) 74 mean_2 <- abs(mean(v[v>=0])) 75 coefs <- model_1$coefficients 76 if (mean_1 < mean_2) { 77 coefs <- model_2$coefficients 78 } 79 list(coefficients=coefs, 80 rmse=rmse(estimated,actual), 81 lag=lag, 82 actual=actual, 83 estimated=estimated) 84 } 85 86 fit_direction_bilateral <- function (direction, plot=TRUE, timestamp=0, dist="vonmises") { 87 hist <- normalize_hist(hist(direction, plot=FALSE, breaks=1000)) 88 x = hist$x; 89 v = hist$density; 90 d = density(direction,adjust=0.5) 91 dist_1 <- fit_direction(x[-pi<=x&x<=-pi/2], v[-pi<=x&x<=-pi/2], plot=FALSE) 92 dist_2 <- fit_direction(x[-pi/2<=x&x<=0], v[-pi/2<=x&x<=0], plot=FALSE) 93 dist_3 <- fit_direction(x[x>=0&x<=pi/2], v[x>=0&x<=pi/2], plot=FALSE) 94 dist_4 <- fit_direction(x[x>=pi/2&x<=pi], v[x>=pi/2&x<=pi], plot=FALSE) 95 #dist_1 <- fit_direction(x[x>=0], v[x>=0], sign=1, plot=FALSE) 96 #x <- c(dist_1$x, dist_3$x) 97 actual_density <- c(dist_1$actual_density, 98 dist_2$actual_density, 99 dist_3$actual_density, 100 dist_4$actual_density) 101 estimated_density <- c(dist_1$estimated_density, 102 dist_2$estimated_density, 103 dist_3$estimated_density, 104 dist_4$estimated_density) 105 if (plot) { 106 plot(x, actual_density, xlab='Direction', ylab='PDF') 107 #hist(direction, xlab=paste('Velocity', axis), ylab='PDF', freq=FALSE, breaks=100) 108 lines(d, col='green', lwd=2) 109 #lines(x, estimated_density, col='red', lwd=2) 110 lines(dist_1$x, dist_1$estimated_density, col='red', lwd=2) 111 lines(dist_3$x, dist_3$estimated_density, col='blue', lwd=2) 112 } 113 #printf("RMSE %f\n", rmse(estimated_density, actual_density)) 114 path <- file.path('build', 'direction', dist) 115 make_directory(path) 116 filename <- file.path(path, sprintf("%d", timestamp)) 117 write.table(data.frame(x=x,actual=actual_density,estimated=estimated_density), 118 filename, row.names=FALSE, quote=FALSE) 119 write("\n",file=filename,append=TRUE) 120 write.table(data.frame(x=d$x,density=d$y), filename, 121 row.names=FALSE, quote=FALSE, append=TRUE) 122 list(coefficients=dist_1$coefficients, 123 rmse=rmse(estimated_density, actual_density)) 124 } 125 126 time_delta = 60*60*2 127 #time_start = time_string_to_timestamp("2021-03-06T00:00") 128 #time_end = time_string_to_timestamp("2021-03-06T23:00") 129 130 # all of the data 131 time_start = time_string_to_timestamp("2021-02-21T00:00") 132 time_end = time_string_to_timestamp("2021-03-31T23:00") 133 134 # timestamp for ACF validation 135 #time_start = 1616936400 136 #time_end = time_start 137 138 print(time_start) 139 print(time_end) 140 141 142 timestamps <- seq(time_start, time_end, time_delta) 143 #pb <- txtProgressBar(min = 0, max = length(timestamps), style = 3) 144 145 #i <- 0 146 #for (timestamp in timestamps) { 147 cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling='remove') %dopar% { 148 cumulative_data <- data.frame(speed=numeric(0),direction=numeric(0), 149 abs_x_mode=numeric(0),x_sd=numeric(0), 150 y_sd=numeric(0),z_sd=numeric(0), 151 x_mean=numeric(0),y_mean=numeric(0),z_mean=numeric(0), 152 x_var=numeric(0),y_var=numeric(0),z_var=numeric(0), 153 v_mean=numeric(0),d_mean=numeric(0), 154 x_a=numeric(0),x_b=numeric(0),x_c=numeric(0), 155 x_d=numeric(0),x_e=numeric(0),x_f=numeric(0), 156 x_mode=numeric(0),y_mode=numeric(0),z_mode=numeric(0), 157 y_a=numeric(0),y_b=numeric(0),y_c=numeric(0), 158 y_d=numeric(0),y_e=numeric(0),y_f=numeric(0), 159 timestamp=numeric(0), 160 x_rmse=numeric(0), y_rmse=numeric(0), z_rmse=numeric(0), 161 acf_x_rmse=numeric(0), acf_y_rmse=numeric(0), acf_z_rmse=numeric(0), 162 acf_x_a=numeric(0),acf_x_b=numeric(0),acf_x_c=numeric(0), 163 acf_y_a=numeric(0),acf_y_b=numeric(0),acf_y_c=numeric(0), 164 acf_z_a=numeric(0),acf_z_b=numeric(0),acf_z_c=numeric(0), 165 acf_v_a=numeric(0),acf_v_b=numeric(0),acf_v_c=numeric(0), 166 acf_d_a=numeric(0),acf_v_b=numeric(0),acf_v_c=numeric(0), 167 x_area_min=numeric(0),x_area_max=numeric(0), 168 y_area_min=numeric(0),y_area_max=numeric(0), 169 z_area_min=numeric(0),z_area_max=numeric(0)) 170 velocity <- select_samples(timestamp, time_delta) 171 if (nrow(velocity) == 0) { return(NULL) } 172 filename <- sprintf("build/rplot-%d.pdf", timestamp) 173 print(filename) 174 pdf(filename) 175 x_mode <- getmode(velocity$x) 176 y_mode <- getmode(velocity$y) 177 z_mode <- getmode(velocity$z) 178 new_row = nrow(cumulative_data)+1 179 cumulative_data[new_row, "timestamp"] = timestamp 180 cumulative_data[new_row, "x_mode"] = x_mode 181 cumulative_data[new_row, "y_mode"] = y_mode 182 cumulative_data[new_row, "z_mode"] = z_mode 183 cumulative_data[new_row, "speed"] = mean(velocity$speed) 184 cumulative_data[new_row, "direction"] = direction(x_mode, y_mode) 185 186 max_abs_mean <- function (x) { 187 #max(abs(mean(x[x<=0])), abs(mean(x[x>=0]))) 188 mean_1 <- abs(mean(x[x<=0])) 189 mean_2 <- abs(mean(x[x>=0])) 190 if (mean_1 > mean_2) { 191 return(mean(x[x<=0])) 192 } else { 193 return(mean(x[x>=0])) 194 } 195 } 196 197 max_abs_var <- function (x) { 198 mean_1 <- abs(mean(x[x<=0])) 199 mean_2 <- abs(mean(x[x>=0])) 200 if (mean_1 > mean_2) { 201 return(var(x[x<=0])) 202 } else { 203 return(var(x[x>=0])) 204 } 205 } 206 207 cumulative_data[new_row, "abs_x_mode"] = abs(x_mode) 208 cumulative_data[new_row, "x_sd"] = sd(velocity$x) 209 cumulative_data[new_row, "y_sd"] = sd(velocity$y) 210 cumulative_data[new_row, "z_sd"] = sd(velocity$z) 211 cumulative_data[new_row, "x_mean"] = max_abs_mean(velocity$x) 212 cumulative_data[new_row, "y_mean"] = max_abs_mean(velocity$y) 213 cumulative_data[new_row, "z_mean"] = max_abs_mean(velocity$z) 214 cumulative_data[new_row, "x_var"] = max_abs_var(velocity$x) 215 cumulative_data[new_row, "y_var"] = max_abs_var(velocity$y) 216 cumulative_data[new_row, "z_var"] = max_abs_var(velocity$z) 217 cumulative_data[new_row, "v_mean"] = mean(velocity$speed) 218 cumulative_data[new_row, "d_mean"] = mean(velocity$direction) 219 #cumulative_data[new_row, "direction"] = weighted.mean(velocity$direction, speed(velocity$x,velocity$y)) 220 #fit <- fitdist(velocity$x, "sn", start=list(omega=1,alpha=-1)) 221 #fit <- fitdist(velocity$x, "weibull") 222 #fit <- fitdist(velocity$x, "expnorm", start=list(a=1,b=1,c=1)) 223 #fit <- fitdist(velocity$x, "norm") 224 #fit <- fitdist(velocity$x, "gumbel", start=list(a=1, b=1)) 225 #plot(fit) 226 # Velocity X 227 #for (dist in c("weibull2")) { 228 # fit_velocity_pdf_bilateral(velocity$x, axis="x", timestamp=timestamp, dist="weibull2") 229 # fit_velocity_pdf_bilateral(velocity$y, axis="y", timestamp=timestamp, dist="weibull2") 230 # fit_velocity_pdf_bilateral(velocity$z, axis="z", timestamp=timestamp, dist="weibull2") 231 #} 232 x_model <- fit_velocity_pdf_bilateral(velocity$x, axis="x", timestamp=timestamp, dist="weibull") 233 cumulative_data[new_row, "x_a"] = x_model$coefficients[[1]] 234 cumulative_data[new_row, "x_b"] = x_model$coefficients[[2]] 235 cumulative_data[new_row, "x_rmse"] = x_model$rmse 236 cumulative_data[new_row, "x_area_min"] = x_model$area_min 237 cumulative_data[new_row, "x_area_max"] = x_model$area_max 238 #cumulative_data[new_row, "x_c"] = x_coefs[[3]] 239 #cumulative_data[new_row, "x_d"] = x_coefs[[4]] 240 #cumulative_data[new_row, "x_e"] = x_coefs[[5]] 241 #cumulative_data[new_row, "x_f"] = x_coefs[[6]] 242 # Velocity Y 243 y_model <- fit_velocity_pdf_bilateral(velocity$y, axis="y", timestamp=timestamp, dist="weibull") 244 cumulative_data[new_row, "y_a"] = y_model$coefficients[[1]] 245 cumulative_data[new_row, "y_b"] = y_model$coefficients[[2]] 246 cumulative_data[new_row, "y_rmse"] = y_model$rmse 247 cumulative_data[new_row, "y_area_min"] = y_model$area_min 248 cumulative_data[new_row, "y_area_max"] = y_model$area_max 249 #cumulative_data[new_row, "y_c"] = y_coefs[[3]] 250 #cumulative_data[new_row, "y_d"] = y_coefs[[4]] 251 #cumulative_data[new_row, "y_e"] = y_coefs[[5]] 252 #cumulative_data[new_row, "y_f"] = y_coefs[[6]] 253 # Velocity Z 254 z_model <- fit_velocity_pdf_bilateral(velocity$z, axis="z", timestamp=timestamp, dist="weibull") 255 cumulative_data[new_row, "z_rmse"] = z_model$rmse 256 cumulative_data[new_row, "z_area_min"] = z_model$area_min 257 cumulative_data[new_row, "z_area_max"] = z_model$area_max 258 # ACF X 259 x_model <- fit_velocity_acf_bilateral(velocity$x, plot=TRUE, axis="x", timestamp=timestamp) 260 cumulative_data[new_row, "acf_x_a"] = x_model$coefficients[[1]] 261 cumulative_data[new_row, "acf_x_b"] = x_model$coefficients[[2]] 262 cumulative_data[new_row, "acf_x_c"] = x_model$coefficients[[3]] 263 cumulative_data[new_row, "acf_x_rmse"] = x_model$rmse 264 y_model <- fit_velocity_acf_bilateral(velocity$y, plot=TRUE, axis="y", timestamp=timestamp) 265 cumulative_data[new_row, "acf_y_a"] = y_model$coefficients[[1]] 266 cumulative_data[new_row, "acf_y_b"] = y_model$coefficients[[2]] 267 cumulative_data[new_row, "acf_y_c"] = y_model$coefficients[[3]] 268 cumulative_data[new_row, "acf_y_rmse"] = y_model$rmse 269 z_model <- fit_velocity_acf_bilateral(velocity$z, plot=TRUE, axis="z", timestamp=timestamp) 270 cumulative_data[new_row, "acf_z_a"] = z_model$coefficients[[1]] 271 cumulative_data[new_row, "acf_z_b"] = z_model$coefficients[[2]] 272 cumulative_data[new_row, "acf_z_c"] = z_model$coefficients[[3]] 273 cumulative_data[new_row, "acf_z_rmse"] = z_model$rmse 274 v_model <- fit_velocity_acf(velocity$speed, plot=TRUE, axis="speed", timestamp=timestamp) 275 cumulative_data[new_row, "acf_v_a"] = v_model$coefficients[[1]] 276 cumulative_data[new_row, "acf_v_b"] = v_model$coefficients[[2]] 277 cumulative_data[new_row, "acf_v_c"] = v_model$coefficients[[3]] 278 cumulative_data[new_row, "acf_v_rmse"] = v_model$rmse 279 d_model <- fit_velocity_acf(velocity$direction, plot=TRUE, axis="direction", timestamp=timestamp) 280 cumulative_data[new_row, "acf_d_a"] = d_model$coefficients[[1]] 281 cumulative_data[new_row, "acf_d_b"] = d_model$coefficients[[2]] 282 cumulative_data[new_row, "acf_d_c"] = d_model$coefficients[[3]] 283 cumulative_data[new_row, "acf_d_rmse"] = d_model$rmse 284 285 direction_model <- fit_direction_bilateral(velocity$direction, timestamp=timestamp) 286 cumulative_data[new_row, "direction_rmse"] = direction_model$rmse 287 288 old_mar = par("mar") 289 direction_hist <- hist(velocity$direction, plot=FALSE, breaks=100) 290 direction_hist <- data.frame(direction=direction_hist$density, angle=direction_hist$mids) 291 #path <- file.path('build', 'direction') 292 #make_directory(path) 293 #filename <- file.path(path, sprintf("%d-hist-xy", timestamp)) 294 ## add corner points with mean speed 295 #x_max = max(abs(velocity$x)) 296 #y_max = max(abs(velocity$y)) 297 #s = mean(velocity$speed) 298 #r = sqrt(x_max**2 + y_max**2) 299 #write.table(data.frame(x=c(velocity$x/r,r,-r,-r,r), 300 # y=c(velocity$y/r,r,r,-r,-r), 301 # r=c(velocity$speed,s,s,s,s)), 302 # filename, row.names=FALSE, quote=FALSE) 303 #hist_xy <- hist2d(velocity$x, velocity$y, same.scale=TRUE, show=TRUE) 304 #grid_xy <- expand.grid(x=hist_xy$x,y=hist_xy$y) 305 #x_max = max(abs(grid_xy$x)) 306 #y_max = max(abs(grid_xy$y)) 307 #write.table(data.frame(x=c(grid_xy$x/r,r,-r,-r,r), 308 # y=c(grid_xy$y/r,r,r,-r,-r), 309 # r=c(as.vector(hist_xy$counts),0,0,0,0)), 310 # filename, row.names=FALSE, quote=FALSE) 311 polar.plot(direction_hist$direction,rad2deg(direction_hist$angle), 312 start=45,clockwise=TRUE,rp.type="polygon") 313 par(mar=old_mar) 314 #i <- i + 1 315 #setTxtProgressBar(pb, i) 316 dev.off() 317 cumulative_data 318 } 319 #close(pb) 320 321 #rmse_table <- cumulative_data[c("timestamp","x_rmse","y_rmse","z_rmse","speed", 322 # "direction_rmse","direction", 323 # "acf_x_rmse","acf_y_rmse","acf_z_rmse")] 324 write.table(cumulative_data, "build/daily-rmse", row.names=FALSE, quote=FALSE) 325 print("Mean speed") 326 print(mean(cumulative_data$speed)) 327 #print(cumulative_data) 328 print(summary(cumulative_data)) 329 model <- nls(x_sd ~ a*speed+b, 330 data=cumulative_data, 331 control=list(warnOnly=TRUE), 332 start=list(a=1,b=0)) 333 coefs <- summary(model)$coefficients 334 x_sd_est <- sapply(cumulative_data$speed, function (x) { coefs[[1]]*x + coefs[[2]] }) 335 plot(cumulative_data$speed, cumulative_data$x_sd) 336 lines(cumulative_data$speed, x_sd_est, col='red') 337 338 #plot(cumulative_data$acf_x_a, cumulative_data$x_sd) 339 #plot(log(cumulative_data$acf_x_b), cumulative_data$x_sd) 340 #plot(cumulative_data$acf_x_c, cumulative_data$x_sd) 341 342 # turbulence coefficient 343 #plot(hist(cumulative_data$x_area_min)) 344 #plot(hist(cumulative_data$x_area_max)) 345 turbulence_x = cumulative_data$x_area_min / cumulative_data$x_area_max 346 turbulence_y = cumulative_data$y_area_min / cumulative_data$y_area_max 347 turbulence_z = cumulative_data$z_area_min / cumulative_data$z_area_max 348 plot(hist(turbulence_x, breaks=50)) 349 plot(abs(cumulative_data$x_mean),turbulence_x) 350 reg <- lm(turbulence_x~cumulative_data$x_mean) 351 lines(cumulative_data$x_mean, reg$fitted.values, col='red') 352 print(reg) 353 print(mean(turbulence_x)) 354 #plot(loess(turbulence_x ~ abs(x_mean), cumulative_data)) 355 scatter.smooth(abs(cumulative_data$x_mean), turbulence_x) 356 path <- file.path('build', 'turbulence') 357 write.table(data.frame(velocity_x=abs(cumulative_data$x_mean), 358 velocity_y=abs(cumulative_data$y_mean), 359 velocity_z=abs(cumulative_data$z_mean), 360 turbulence_x=turbulence_x, 361 turbulence_y=turbulence_y, 362 turbulence_z=turbulence_z), 363 path, row.names=FALSE, quote=FALSE) 364 #lo <- loess.smooth(cumulative_data$x_mean, cumulative_data$turbulence_x, span=0.5) 365 #lines(lo$x, lo$y) 366 367 # remove outliers 368 cumulative_data <- cumulative_data[cumulative_data$acf_x_b<2 & cumulative_data$acf_y_b<2 & cumulative_data$acf_z_b<2,] 369 370 plot(cumulative_data$x_mean, (cumulative_data$acf_x_a)) 371 reg <- lm((cumulative_data$acf_x_a)~cumulative_data$x_mean) 372 lines(cumulative_data$x_mean, reg$fitted.values, col='red') 373 print(reg) 374 375 plot(cumulative_data$x_mean, (cumulative_data$acf_x_b)) 376 reg <- lm((cumulative_data$acf_x_b)~cumulative_data$x_mean) 377 lines(cumulative_data$x_mean, reg$fitted.values, col='red') 378 print(reg) 379 380 plot(cumulative_data$x_mean, cumulative_data$acf_x_c) 381 reg <- lm(cumulative_data$acf_x_c~cumulative_data$x_mean) 382 lines(cumulative_data$x_mean, reg$fitted.values, col='red') 383 print(reg) 384 385 # ### 386 387 #plot(cumulative_data$x_mean, (cumulative_data$acf_y_a)) 388 #reg <- lm((cumulative_data$acf_y_a)~cumulative_data$x_mean) 389 #lines(cumulative_data$x_mean, reg$fitted.values, col='red') 390 #print(reg) 391 # 392 #plot(cumulative_data$x_mean, (cumulative_data$acf_y_b)) 393 #reg <- lm((cumulative_data$acf_y_b)~cumulative_data$x_mean) 394 #lines(cumulative_data$x_mean, reg$fitted.values, col='red') 395 #print(reg) 396 # 397 #plot(cumulative_data$x_mean, cumulative_data$acf_y_c) 398 #reg <- lm(cumulative_data$acf_y_c~cumulative_data$x_mean) 399 #lines(cumulative_data$x_mean, reg$fitted.values, col='red') 400 #print(reg) 401 402 403 plot(cumulative_data$y_mean, (cumulative_data$acf_y_a)) 404 reg <- lm((cumulative_data$acf_y_a)~cumulative_data$y_mean) 405 lines(cumulative_data$y_mean, reg$fitted.values, col='red') 406 print(reg) 407 408 plot(cumulative_data$y_mean, (cumulative_data$acf_y_b)) 409 reg <- lm((cumulative_data$acf_y_b)~cumulative_data$y_mean) 410 lines(cumulative_data$y_mean, reg$fitted.values, col='red') 411 print(reg) 412 413 plot(cumulative_data$y_mean, cumulative_data$acf_y_c) 414 reg <- lm(cumulative_data$acf_y_c~cumulative_data$y_mean) 415 lines(cumulative_data$y_mean, reg$fitted.values, col='red') 416 print(reg) 417 418 419 plot(cumulative_data$z_mean, (cumulative_data$acf_z_a)) 420 reg <- lm((cumulative_data$acf_z_a)~cumulative_data$z_mean) 421 lines(cumulative_data$z_mean, reg$fitted.values, col='red') 422 print(reg) 423 424 plot(cumulative_data$z_mean, (cumulative_data$acf_z_b)) 425 reg <- lm((cumulative_data$acf_z_b)~cumulative_data$z_mean) 426 lines(cumulative_data$z_mean, reg$fitted.values, col='red') 427 print(reg) 428 429 plot(cumulative_data$z_mean, cumulative_data$acf_z_c) 430 reg <- lm(cumulative_data$acf_z_c~cumulative_data$z_mean) 431 lines(cumulative_data$z_mean, reg$fitted.values, col='red') 432 print(reg) 433 434 #plot(cumulative_data$x_mean, (cumulative_data$acf_x_b)) 435 quit() 436 #plot(cumulative_data$acf_v_a, cumulative_data$v_mean) 437 #plot(log(cumulative_data$acf_v_b), cumulative_data$v_mean) 438 #plot(cumulative_data$acf_v_c, cumulative_data$v_mean) 439 #plot(cumulative_data$acf_d_a, cumulative_data$d_mean) 440 #plot(cumulative_data$acf_d_b, cumulative_data$d_mean) 441 #plot(cumulative_data$acf_d_c, cumulative_data$d_mean) 442 443 plot(cumulative_data$abs_x_mode, cumulative_data$acf_x_b) 444 #plot(cumulative_data$acf_x_b, cumulative_data$y_mean/tan(cumulative_data$direction)) 445 #plot(cumulative_data$speed, cumulative_data$acf_x_b) 446 #plot(cumulative_data$direction, cumulative_data$acf_x_b) 447 model <- nls(acf_x_b ~ a1*abs_x_mode+b1, 448 data=cumulative_data, 449 control=list(warnOnly=TRUE), 450 start=list(a1=1,b1=0), 451 algorithm="port", 452 lower=c(-1000,-1000), upper=c(1000,1000)) 453 coefs <- summary(model)$coefficients 454 print(coefs) 455 b_est <- sapply(cumulative_data$abs_x_mode, function (x) { coefs[[1]]*x + coefs[[2]] }) 456 lines(cumulative_data$abs_x_mode, b_est, col='red') 457 458 459 cumulative_data$abs_x_mode_acf_x_c <- cumulative_data$abs_x_mode**-cumulative_data$acf_x_c 460 plot(cumulative_data$abs_x_mode_acf_x_c, exp(-cumulative_data$acf_x_b**(cumulative_data$acf_x_c))) 461 model <- nls(exp(-acf_x_b**acf_x_c) ~ abs(a1*abs_x_mode_acf_x_c+b1), 462 data=cumulative_data, 463 control=list(warnOnly=TRUE), 464 start=list(a1=1,b1=0), 465 algorithm="port", 466 lower=c(-1000,-1000), upper=c(1000,1000)) 467 coefs <- summary(model)$coefficients 468 print(coefs) 469 #bc_est <- sapply(cumulative_data$abs_x_mode, function (x) { (coefs[[1]]*x + coefs[[2]]) }) 470 cumulative_data <- cumulative_data[order(cumulative_data$abs_x_mode_acf_x_c), ] 471 bc_est <- with (cumulative_data, { 472 (coefs[[1]]*abs_x_mode_acf_x_c + coefs[[2]]) 473 }) 474 lines(cumulative_data$abs_x_mode_acf_x_c, bc_est, col='red') 475 476 477 478 479 #plot(cumulative_data$x_mode, cumulative_data$c) 480 #plot(cumulative_data$x_mode, cumulative_data$acf_x_b**cumulative_data$c) 481 482 #plot(sign(cumulative_data$x_mode)*(cumulative_data$x_mode)**2, cumulative_data$c) 483 #plot(sign(cumulative_data$x_mode)*sqrt(abs(cumulative_data$x_mode)), cumulative_data$c) 484 #plot(log(cumulative_data$speed), cumulative_data$c) 485 #plot(cumulative_data$direction, cumulative_data$x_b) 486 #plot(cumulative_data$direction, cumulative_data$y_b) 487 #plot(cumulative_data$direction, cumulative_data$acf_x_b) 488 #plot(cumulative_data$direction, cumulative_data$acf_y_b) 489 plot(cumulative_data$abs_x_mode, cumulative_data$x_c) 490 plot(cumulative_data$abs_x_mode, cumulative_data$acf_x_c) 491 #cumulative_data$x_b <- ifelse(cumulative_data$x_b > 10, 1, cumulative_data$x_b) 492 #print(summary(cumulative_data$x_b)) 493 with (cumulative_data, { 494 plot((abs_x_mode**-acf_x_c), exp(-x_b**(x_c))) 495 plot((abs_x_mode**-acf_x_c), exp(-acf_x_b**(acf_x_c))) 496 plot(direction, direction(acf_x_b**acf_x_c, acf_y_b**acf_x_c)) 497 plot(direction, direction(x_b**acf_x_c, y_b**acf_x_c)) 498 }) 499