commit 2c275e446ed303e28087642fe1b54716bd0a250a
parent 1895895637899bca6beecc257289ce4460b28c84
Author: Ivan Gankevich <i.gankevich@spbu.ru>
Date: Mon, 19 Apr 2021 13:50:18 +0300
Compute ACF coefficient regressions on mean wind speed.
Diffstat:
| R/anal.R | | | 137 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------- |
1 file changed, 123 insertions(+), 14 deletions(-)
diff --git a/R/anal.R b/R/anal.R
@@ -325,7 +325,7 @@ fit_velocity_acf_bilateral <- function (velocity, axis="x", timestamp=0, plot=TR
estimated <- c(model_1$estimated, model_2$estimated)
actual <- c(model_1$actual, model_2$actual)
plot(lag, actual, xlab='Lag', ylab='ACF')
- lines(-model_1$lag, model_2$estimated, col='red', lwd=2)
+ lines(-model_1$lag, model_1$estimated, col='red', lwd=2)
lines(model_2$lag, model_2$estimated, col='blue', lwd=2)
path <- file.path('build', 'acf', axis)
make_directory(path)
@@ -337,7 +337,13 @@ fit_velocity_acf_bilateral <- function (velocity, axis="x", timestamp=0, plot=TR
write("\n",file=filename,append=TRUE)
write.table(data.frame(lag=model_2$lag,actual=model_2$actual,estimated=model_2$estimated),
filename, row.names=FALSE, quote=FALSE, append=TRUE)
- list(coefficients=model_1$coefficients,
+ mean_1 <- abs(mean(v[v<=0]))
+ mean_2 <- abs(mean(v[v>=0]))
+ coefs <- model_1$coefficients
+ if (mean_1 < mean_2) {
+ coefs <- model_2$coefficients
+ }
+ list(coefficients=coefs,
rmse=rmse(estimated,actual),
lag=lag,
actual=actual,
@@ -413,11 +419,13 @@ fit_direction_bilateral <- function (direction, plot=TRUE, timestamp=0, dist="vo
rmse=rmse(estimated_density, actual_density))
}
+time_delta = 60*60*2
#time_start = time_string_to_timestamp("2021-03-06T00:00")
#time_end = time_string_to_timestamp("2021-03-06T23:00")
-time_start = time_string_to_timestamp("2021-02-21T00:00")
-time_end = time_string_to_timestamp("2021-03-31T23:00")
-time_delta = 60*60*2
+#time_start = time_string_to_timestamp("2021-02-21T00:00")
+#time_end = time_string_to_timestamp("2021-03-31T23:00")
+time_start = 1616936400
+time_end = time_start
print(time_start)
print(time_end)
@@ -430,12 +438,14 @@ timestamps <- seq(time_start, time_end, time_delta)
#for (timestamp in timestamps) {
cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling='remove') %dopar% {
cumulative_data <- data.frame(speed=numeric(0),direction=numeric(0),
- x_mode=numeric(0),abs_x_mode=numeric(0),x_sd=numeric(0),
- y_mean=numeric(0),v_mean=numeric(0),d_mean=numeric(0),
+ abs_x_mode=numeric(0),x_sd=numeric(0),
+ y_sd=numeric(0),z_sd=numeric(0),
+ x_mean=numeric(0),y_mean=numeric(0),z_mean=numeric(0),
+ x_var=numeric(0),y_var=numeric(0),z_var=numeric(0),
+ v_mean=numeric(0),d_mean=numeric(0),
x_a=numeric(0),x_b=numeric(0),x_c=numeric(0),
x_d=numeric(0),x_e=numeric(0),x_f=numeric(0),
- # Y
- y_mode=numeric(0),
+ x_mode=numeric(0),y_mode=numeric(0),z_mode=numeric(0),
y_a=numeric(0),y_b=numeric(0),y_c=numeric(0),
y_d=numeric(0),y_e=numeric(0),y_f=numeric(0),
timestamp=numeric(0),
@@ -453,10 +463,12 @@ cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling
pdf(filename)
x_mode <- getmode(velocity$x)
y_mode <- getmode(velocity$y)
+ z_mode <- getmode(velocity$z)
new_row = nrow(cumulative_data)+1
cumulative_data[new_row, "timestamp"] = timestamp
cumulative_data[new_row, "x_mode"] = x_mode
cumulative_data[new_row, "y_mode"] = y_mode
+ cumulative_data[new_row, "z_mode"] = z_mode
cumulative_data[new_row, "speed"] = mean(velocity$speed)
cumulative_data[new_row, "direction"] = direction(x_mode, y_mode)
#cumulative_data[new_row, "direction"] = weighted.mean(velocity$direction, speed(velocity$x,velocity$y))
@@ -549,9 +561,37 @@ cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling
start=45,clockwise=TRUE,rp.type="polygon")
par(mar=old_mar)
+ max_abs_mean <- function (x) {
+ #max(abs(mean(x[x<=0])), abs(mean(x[x>=0])))
+ mean_1 <- abs(mean(x[x<=0]))
+ mean_2 <- abs(mean(x[x>=0]))
+ if (mean_1 > mean_2) {
+ return(mean(x[x<=0]))
+ } else {
+ return(mean(x[x>=0]))
+ }
+ }
+
+ max_abs_var <- function (x) {
+ mean_1 <- abs(mean(x[x<=0]))
+ mean_2 <- abs(mean(x[x>=0]))
+ if (mean_1 > mean_2) {
+ return(var(x[x<=0]))
+ } else {
+ return(var(x[x>=0]))
+ }
+ }
+
cumulative_data[new_row, "abs_x_mode"] = abs(x_mode)
cumulative_data[new_row, "x_sd"] = sd(velocity$x)
- cumulative_data[new_row, "y_mean"] = mean(velocity$y)
+ cumulative_data[new_row, "y_sd"] = sd(velocity$y)
+ cumulative_data[new_row, "z_sd"] = sd(velocity$z)
+ cumulative_data[new_row, "x_mean"] = max_abs_mean(velocity$x)
+ cumulative_data[new_row, "y_mean"] = max_abs_mean(velocity$y)
+ cumulative_data[new_row, "z_mean"] = max_abs_mean(velocity$z)
+ cumulative_data[new_row, "x_var"] = max_abs_var(velocity$x)
+ cumulative_data[new_row, "y_var"] = max_abs_var(velocity$y)
+ cumulative_data[new_row, "z_var"] = max_abs_var(velocity$z)
cumulative_data[new_row, "v_mean"] = mean(velocity$speed)
cumulative_data[new_row, "d_mean"] = mean(velocity$direction)
#i <- i + 1
@@ -572,14 +612,83 @@ model <- nls(x_sd ~ a*speed+b,
control=list(warnOnly=TRUE),
start=list(a=1,b=0))
coefs <- summary(model)$coefficients
-#print(coefs)
x_sd_est <- sapply(cumulative_data$speed, function (x) { coefs[[1]]*x + coefs[[2]] })
plot(cumulative_data$speed, cumulative_data$x_sd)
lines(cumulative_data$speed, x_sd_est, col='red')
-plot(cumulative_data$acf_x_a, cumulative_data$x_sd)
-plot(log(cumulative_data$acf_x_b), cumulative_data$x_sd)
-plot(cumulative_data$acf_x_c, cumulative_data$x_sd)
+#plot(cumulative_data$acf_x_a, cumulative_data$x_sd)
+#plot(log(cumulative_data$acf_x_b), cumulative_data$x_sd)
+#plot(cumulative_data$acf_x_c, cumulative_data$x_sd)
+
+# remove outliers
+cumulative_data <- cumulative_data[cumulative_data$acf_x_b<2 & cumulative_data$acf_y_b<2 & cumulative_data$acf_z_b<2,]
+
+plot(cumulative_data$x_mean, (cumulative_data$acf_x_a))
+reg <- lm((cumulative_data$acf_x_a)~cumulative_data$x_mean)
+lines(cumulative_data$x_mean, reg$fitted.values, col='red')
+print(reg)
+
+plot(cumulative_data$x_mean, (cumulative_data$acf_x_b))
+reg <- lm((cumulative_data$acf_x_b)~cumulative_data$x_mean)
+lines(cumulative_data$x_mean, reg$fitted.values, col='red')
+print(reg)
+
+plot(cumulative_data$x_mean, cumulative_data$acf_x_c)
+reg <- lm(cumulative_data$acf_x_c~cumulative_data$x_mean)
+lines(cumulative_data$x_mean, reg$fitted.values, col='red')
+print(reg)
+
+# ###
+
+#plot(cumulative_data$x_mean, (cumulative_data$acf_y_a))
+#reg <- lm((cumulative_data$acf_y_a)~cumulative_data$x_mean)
+#lines(cumulative_data$x_mean, reg$fitted.values, col='red')
+#print(reg)
+#
+#plot(cumulative_data$x_mean, (cumulative_data$acf_y_b))
+#reg <- lm((cumulative_data$acf_y_b)~cumulative_data$x_mean)
+#lines(cumulative_data$x_mean, reg$fitted.values, col='red')
+#print(reg)
+#
+#plot(cumulative_data$x_mean, cumulative_data$acf_y_c)
+#reg <- lm(cumulative_data$acf_y_c~cumulative_data$x_mean)
+#lines(cumulative_data$x_mean, reg$fitted.values, col='red')
+#print(reg)
+
+
+plot(cumulative_data$y_mean, (cumulative_data$acf_y_a))
+reg <- lm((cumulative_data$acf_y_a)~cumulative_data$y_mean)
+lines(cumulative_data$y_mean, reg$fitted.values, col='red')
+print(reg)
+
+plot(cumulative_data$y_mean, (cumulative_data$acf_y_b))
+reg <- lm((cumulative_data$acf_y_b)~cumulative_data$y_mean)
+lines(cumulative_data$y_mean, reg$fitted.values, col='red')
+print(reg)
+
+plot(cumulative_data$y_mean, cumulative_data$acf_y_c)
+reg <- lm(cumulative_data$acf_y_c~cumulative_data$y_mean)
+lines(cumulative_data$y_mean, reg$fitted.values, col='red')
+print(reg)
+
+
+plot(cumulative_data$z_mean, (cumulative_data$acf_z_a))
+reg <- lm((cumulative_data$acf_z_a)~cumulative_data$z_mean)
+lines(cumulative_data$z_mean, reg$fitted.values, col='red')
+print(reg)
+
+plot(cumulative_data$z_mean, (cumulative_data$acf_z_b))
+reg <- lm((cumulative_data$acf_z_b)~cumulative_data$z_mean)
+lines(cumulative_data$z_mean, reg$fitted.values, col='red')
+print(reg)
+
+plot(cumulative_data$z_mean, cumulative_data$acf_z_c)
+reg <- lm(cumulative_data$acf_z_c~cumulative_data$z_mean)
+lines(cumulative_data$z_mean, reg$fitted.values, col='red')
+print(reg)
+
+#plot(cumulative_data$x_mean, (cumulative_data$acf_x_b))
+quit()
#plot(cumulative_data$acf_v_a, cumulative_data$v_mean)
#plot(log(cumulative_data$acf_v_b), cumulative_data$v_mean)
#plot(cumulative_data$acf_v_c, cumulative_data$v_mean)
