commit 14baeed498060508184bbf202d31685fecdf27bc
parent b0f5903d3345be34cb8288e6c9a06945f057931e
Author: Ivan Gankevich <i.gankevich@spbu.ru>
Date: Thu, 1 Apr 2021 17:30:07 +0300
wind direction bilateral distribution
Diffstat:
| R/anal.R | | | 73 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------- |
1 file changed, 54 insertions(+), 19 deletions(-)
diff --git a/R/anal.R b/R/anal.R
@@ -308,34 +308,69 @@ fit_velocity_acf <- function (x, plot=TRUE, axis="x") {
c
}
-fit_direction <- function (x, x_mode, plot=TRUE) {
- pdf_x <- hist(x-x_mode, plot=FALSE, breaks=1000)
- pdf_x <- data.frame(density=pdf_x$density/max(abs(pdf_x$density)),
- x=pdf_x$mids/max(abs(pdf_x$mids)))
+fit_direction <- function (theta, density, plot=TRUE, sign=-1) {
+ pdf_x <- data.frame(x=theta,density=density);
+ #print(pdf_x)
+ scale <- max(abs(theta))
c <- tryCatch({
- model <- nls(density ~ dvonmises(x,0,a,b,c),
+ model <- nls(density ~ dvonmises(2*pi*x/scale-m,0,a,b,c),
data=pdf_x,
- start=list(a=1,b=1,c=1),
+ start=list(a=1,b=1,c=1,m=0),
control=list(warnOnly=TRUE),
algorithm="port",
- lower=c(0,1,0), upper=c(1000,1000,1000))
+ lower=c(0,1,0,-10),
+ upper=c(1000,1000,1000,10))
c <- summary(model)$coefficients
},error=function(cond) {
message(cond)
c(1,1,1)
})
+ pdf_x_est <- dvonmises(2*pi*pdf_x$x/scale,c[[4]],c[[1]],c[[2]],c[[3]])
if (plot) {
- pdf_x_est <- sapply(pdf_x$x, function (x) { dvonmises(x,0,c[[1]],c[[2]],c[[3]]) })
- plot(pdf_x$x+x_mode, pdf_x$density, xlab='Direction', ylab='PDF')
- lines(pdf_x$x+x_mode, pdf_x_est, col='red')
+ plot(pdf_x$x, pdf_x$density, xlab='Direction', ylab='PDF')
+ lines(pdf_x$x, pdf_x_est, col='red')
}
- c
+ list(coefficients=c,
+ actual_density=pdf_x$density,
+ estimated_density=pdf_x_est,
+ x=pdf_x$x,
+ residual=max(abs(pdf_x$density-pdf_x_est)))
+}
+
+fit_direction_bilateral <- function (direction, plot=TRUE) {
+ hist <- normalize_hist(hist(direction, plot=FALSE, breaks=1000))
+ x = hist$x;
+ v = hist$density;
+ d = density(direction,adjust=0.5)
+ dist_left <- fit_direction(x[x<=0], v[x<=0], sign=-1, plot=FALSE)
+ dist_right <- fit_direction(x[x>=0], v[x>=0], sign=1, plot=FALSE)
+ #dist_left <- fit_direction(x[x>=0], v[x>=0], sign=1, plot=FALSE)
+ #x <- c(dist_left$x, dist_right$x)
+ actual_density <- c(dist_left$actual_density, dist_right$actual_density)
+ estimated_density <- c(dist_left$estimated_density, dist_right$estimated_density)
+ plot(x, actual_density, xlab='Direction', ylab='PDF')
+ #hist(direction, xlab=paste('Velocity', axis), ylab='PDF', freq=FALSE, breaks=100)
+ lines(d, col='green', lwd=2)
+ #lines(x, estimated_density, col='red', lwd=2)
+ lines(dist_left$x, dist_left$estimated_density, col='red', lwd=2)
+ lines(dist_right$x, dist_right$estimated_density, col='blue', lwd=2)
+ #printf("RMSE %f\n", rmse(estimated_density, actual_density))
+ path <- file.path('build', 'direction', dist)
+ dir.create(path, recursive=TRUE)
+ filename <- file.path(path, sprintf("%d", timestamp))
+ write.table(data.frame(x=x,actual=actual_density,estimated=estimated_density),
+ filename, row.names=FALSE, quote=FALSE)
+ write("\n",file=filename,append=TRUE)
+ write.table(data.frame(x=d$x,density=d$y), filename,
+ row.names=FALSE, quote=FALSE, append=TRUE)
+ list(coefficients=dist_left$coefficients,
+ rmse=rmse(estimated_density, actual_density))
}
-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_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
print(time_start)
@@ -347,7 +382,7 @@ timestamps <- seq(time_start, time_end, time_delta)
#i <- 0
#for (timestamp in timestamps) {
-cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling='remove') %dopar% {
+cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling='stop') %dopar% {
db <- dbConnect(RSQLite::SQLite(), "samples/load-cell.sqlite3")
cumulative_data <- data.frame(speed=numeric(0),direction=numeric(0),
x_mode=numeric(0),abs_x_mode=numeric(0),x_sd=numeric(0),
@@ -435,8 +470,8 @@ cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling
cumulative_data[new_row, "acf_y_b"] = y_coefs[[2]]
cumulative_data[new_row, "acf_y_c"] = y_coefs[[3]]
- direction_coefs <- fit_direction(velocity$direction, direction_mode, plot=TRUE)
- #print(direction_coefs)
+ direction_model <- fit_direction_bilateral(velocity$direction, plot=TRUE)
+ cumulative_data[new_row, "direction_rmse"] = direction_model$rmse
old_mar = par("mar")
direction_hist <- hist(velocity$direction, plot=FALSE, breaks=100)
@@ -456,7 +491,7 @@ cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling
}
#close(pb)
-rmse_table <- cumulative_data[c("timestamp","x_rmse","y_rmse","z_rmse","speed")]
+rmse_table <- cumulative_data[c("timestamp","x_rmse","y_rmse","z_rmse","speed","direction_rmse")]
write.table(rmse_table, "build/daily-rmse", row.names=FALSE, quote=FALSE)
#print(cumulative_data)
#print(summary(cumulative_data))
