commit 6c53bbcaaf1cd4f409d0ab1c96bd67e352c47140
parent 14baeed498060508184bbf202d31685fecdf27bc
Author: Ivan Gankevich <i.gankevich@spbu.ru>
Date: Tue, 6 Apr 2021 01:42:12 +0300
Direction distribution.
Diffstat:
8 files changed, 241 insertions(+), 58 deletions(-)
diff --git a/Makefile b/Makefile
@@ -18,6 +18,8 @@ build/main.pdf: build/inkscape/anemometer.eps
build/main.pdf: build/daily-stats.eps
build/main.pdf: build/daily-rmse.eps
build/main.pdf: build/gnuplot/velocity-dist.eps
+build/main.pdf: build/gnuplot/direction-dist-max.png
+build/main.pdf: build/gnuplot/direction-dist-min.png
#build/main.pdf: build/gnuplot/anemometer.eps
#build/main.pdf: build/gnuplot/anemometer.svg
build/main.pdf:
@@ -37,6 +39,9 @@ build/iccsa-21-wind-slides.pdf: build/slides.pdf
build/%.eps: build/%.svg
inkscape -z --export-eps=$@ $<
+build/%.png: build/%.svg
+ inkscape -z --export-png=$@ --export-dpi=600 --export-area-drawing $<
+
build/gnuplot/%.svg: gnuplot/%.gnuplot
@mkdir -p build/gnuplot
gnuplot -d $<
@@ -51,6 +56,14 @@ build/daily-stats.svg: gnuplot/daily-stats.gnuplot build/daily-stats
build/daily-rmse.svg: gnuplot/daily-rmse.gnuplot build/daily-rmse
gnuplot -d $<
+build/gnuplot/direction-dist-min.svg: gnuplot/direction-dist.gnuplot
+ gnuplot -d -c gnuplot/direction-dist.gnuplot \
+ build/direction/1616497200-hist-xy build/gnuplot/direction-dist-min.svg
+
+build/gnuplot/direction-dist-max.svg: gnuplot/direction-dist.gnuplot
+ gnuplot -d -c gnuplot/direction-dist.gnuplot \
+ build/direction/1616792400-hist-xy build/gnuplot/direction-dist-max.svg
+
#build/main.zip: build/gnuplot/*.eps main.tex
build/main.zip: main.tex
@mkdir -p build
diff --git a/R/anal.R b/R/anal.R
@@ -4,6 +4,7 @@ library(plotrix)
library(sn)
library(doParallel)
library(foreach)
+library(gplots) # hist2d
registerDoParallel(cores=detectCores())
source("R/common.R")
@@ -24,7 +25,7 @@ time_string_to_timestamp <- function(str) {
as.integer(as.POSIXct(strptime(str, datetime_format)))
}
-speed <- function(x,y,z) { sqrt(x**2 + y**2 + z**2) }
+speed <- function(x,y,z=0) { sqrt(x**2 + y**2 + z**2) }
direction <- function(x,y) { atan2(y,x) }
# retain the original sign
power <- function(x,y) {
@@ -33,8 +34,10 @@ power <- function(x,y) {
#abs(x)**y
}
-rmse <- function(estimated, actual){
- sqrt(mean((estimated - actual)**2))
+# normalized rmse
+rmse <- function(estimated, actual) {
+ s <- max(actual)
+ sqrt(mean(((estimated - actual)/s)**2))
}
dweibull2 <- function (x, a1=1, b1=1, c1=1, a2=1, b2=1, c2=1, g=0) {
@@ -337,7 +340,7 @@ fit_direction <- function (theta, density, plot=TRUE, sign=-1) {
residual=max(abs(pdf_x$density-pdf_x_est)))
}
-fit_direction_bilateral <- function (direction, plot=TRUE) {
+fit_direction_bilateral <- function (direction, plot=TRUE, timestamp=0, dist="vonmises") {
hist <- normalize_hist(hist(direction, plot=FALSE, breaks=1000))
x = hist$x;
v = hist$density;
@@ -348,12 +351,14 @@ fit_direction_bilateral <- function (direction, plot=TRUE) {
#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)
+ if (plot) {
+ 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)
@@ -382,7 +387,7 @@ timestamps <- seq(time_start, time_end, time_delta)
#i <- 0
#for (timestamp in timestamps) {
-cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling='stop') %dopar% {
+cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling='remove') %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),
@@ -407,7 +412,7 @@ cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling
ORDER BY timestamp",
params = list(t0 = timestamp, dt=time_delta,
xmean=mean(velocity$x), ymean=mean(velocity$y), zmean=mean(velocity$z)))
- if (nrow(velocity) == 0) { next }
+ if (nrow(velocity) == 0) { return(NULL) }
filename <- sprintf("build/rplot-%d.pdf", timestamp)
print(filename)
pdf(filename)
@@ -418,8 +423,6 @@ cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling
velocity$z = sample_to_speed(velocity$z, coefficients$z1, coefficients$z2)
x_mode <- getmode(velocity$x)
y_mode <- getmode(velocity$y)
- #direction_mode <- getmode(velocity$direction)
- direction_mode <- direction(x_mode, y_mode)
velocity$speed = speed(velocity$x,velocity$y,velocity$z)
velocity$direction = direction(velocity$x,velocity$y)
new_row = nrow(cumulative_data)+1
@@ -427,7 +430,8 @@ cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling
cumulative_data[new_row, "x_mode"] = x_mode
cumulative_data[new_row, "y_mode"] = y_mode
cumulative_data[new_row, "speed"] = mean(velocity$speed)
- cumulative_data[new_row, "direction"] = direction_mode
+ cumulative_data[new_row, "direction"] = direction(x_mode, y_mode)
+ #cumulative_data[new_row, "direction"] = weighted.mean(velocity$direction, speed(velocity$x,velocity$y))
#fit <- fitdist(velocity$x, "sn", start=list(omega=1,alpha=-1))
#fit <- fitdist(velocity$x, "weibull")
#fit <- fitdist(velocity$x, "expnorm", start=list(a=1,b=1,c=1))
@@ -470,12 +474,34 @@ 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_model <- fit_direction_bilateral(velocity$direction, plot=TRUE)
+ direction_model <- fit_direction_bilateral(velocity$direction, timestamp=timestamp)
cumulative_data[new_row, "direction_rmse"] = direction_model$rmse
old_mar = par("mar")
direction_hist <- hist(velocity$direction, plot=FALSE, breaks=100)
direction_hist <- data.frame(direction=direction_hist$density, angle=direction_hist$mids)
+ path <- file.path('build', 'direction')
+ dir.create(path, recursive=TRUE)
+ filename <- file.path(path, sprintf("%d-hist-xy", timestamp))
+ #
+ hist_xy <- hist2d(velocity$x, velocity$y, same.scale=TRUE, show=TRUE)
+ grid_xy <- expand.grid(x=hist_xy$x,y=hist_xy$y)
+ #
+ # add corner points with mean speed
+ #x_max = max(abs(velocity$x))
+ #y_max = max(abs(velocity$y))
+ x_max = max(abs(grid_xy$x))
+ y_max = max(abs(grid_xy$y))
+ s = mean(velocity$speed)
+ r = sqrt(x_max**2 + y_max**2)
+ #write.table(data.frame(x=c(velocity$x/r,r,-r,-r,r),
+ # y=c(velocity$y/r,r,r,-r,-r),
+ # r=c(velocity$speed,s,s,s,s)),
+ # filename, row.names=FALSE, quote=FALSE)
+ write.table(data.frame(x=c(grid_xy$x/r,r,-r,-r,r),
+ y=c(grid_xy$y/r,r,r,-r,-r),
+ r=c(as.vector(hist_xy$counts),0,0,0,0)),
+ filename, row.names=FALSE, quote=FALSE)
polar.plot(direction_hist$direction,rad2deg(direction_hist$angle),
start=45,clockwise=TRUE,rp.type="polygon")
par(mar=old_mar)
@@ -491,7 +517,7 @@ cumulative_data <- foreach (timestamp=timestamps, .combine=rbind, .errorhandling
}
#close(pb)
-rmse_table <- cumulative_data[c("timestamp","x_rmse","y_rmse","z_rmse","speed","direction_rmse")]
+rmse_table <- cumulative_data[c("timestamp","x_rmse","y_rmse","z_rmse","speed","direction_rmse","direction")]
write.table(rmse_table, "build/daily-rmse", row.names=FALSE, quote=FALSE)
#print(cumulative_data)
#print(summary(cumulative_data))
diff --git a/R/minmax.R b/R/minmax.R
@@ -1,17 +1,38 @@
library(DBI)
+source("R/common.R")
+
data <- read.table("build/daily-rmse", sep="", header=TRUE)
z_rmse_sorted <- sort(data$z_rmse, decreasing=TRUE)
-print("MAX RMSE")
+print("Min. velocity RMSE")
+print(data[data$x_rmse == min(data$x_rmse), ])
+print(data[data$y_rmse == min(data$y_rmse), ])
+print(data[data$z_rmse == min(data$z_rmse), ])
+print("Max. velocity RMSE")
print(data[data$x_rmse == max(data$x_rmse), ])
print(data[data$y_rmse == max(data$y_rmse), ])
print(data[data$z_rmse == z_rmse_sorted[[1]], ])
print(data[data$z_rmse == z_rmse_sorted[[2]], ])
-print("MIN RMSE")
-print(data[data$x_rmse == min(data$x_rmse), ])
-print(data[data$y_rmse == min(data$y_rmse), ])
-print(data[data$z_rmse == min(data$z_rmse), ])
-print("AVG RMSE")
+print("Mean velocity RMSE")
print(summary(data))
print(cor(data))
+with(list(tmp = data[data$x_rmse == min(data$x_rmse), ]),
+ printf("Velocity X min. RMSE = %f (%s), mean speed = %f\n", tmp["x_rmse"], tmp["timestamp"], tmp["speed"]))
+with(list(tmp = data[data$y_rmse == min(data$y_rmse), ]),
+ printf("Velocity Y min. RMSE = %f (%s), mean speed = %f\n", tmp["y_rmse"], tmp["timestamp"], tmp["speed"]))
+with(list(tmp = data[data$z_rmse == min(data$z_rmse), ]),
+ printf("Velocity Z min. RMSE = %f (%s), mean speed = %f\n", tmp["z_rmse"], tmp["timestamp"], tmp["speed"]))
+with(list(tmp = data[data$x_rmse == max(data$x_rmse), ]),
+ printf("Velocity X max. RMSE = %f (%s), mean speed = %f\n", tmp["x_rmse"], tmp["timestamp"], tmp["speed"]))
+with(list(tmp = data[data$y_rmse == max(data$y_rmse), ]),
+ printf("Velocity Y max. RMSE = %f (%s), mean speed = %f\n", tmp["y_rmse"], tmp["timestamp"], tmp["speed"]))
+with(list(tmp = data[data$z_rmse == max(data$z_rmse), ]),
+ printf("Velocity Z max. RMSE = %f (%s), mean speed = %f\n", tmp["z_rmse"], tmp["timestamp"], tmp["speed"]))
#print(data)
+
+direction_rmse_min = data[data$direction_rmse == min(data$direction_rmse), ]
+direction_rmse_max = data[data$direction_rmse == max(data$direction_rmse), ]
+direction_rmse_mean = mean(data$direction_rmse)
+printf("Min. direction RMSE = %f (%s), mean direction = %f\n", direction_rmse_min["direction_rmse"], direction_rmse_min["timestamp"], direction_rmse_min["direction"])
+printf("Max. direction RMSE = %f (%s), mean direction = %f\n", direction_rmse_max["direction_rmse"], direction_rmse_max["timestamp"], direction_rmse_max["direction"])
+printf("Mean direction RMSE = %f\n", direction_rmse_mean)
diff --git a/gnuplot/blues.pal b/gnuplot/blues.pal
@@ -0,0 +1,26 @@
+# line styles for ColorBrewer Blues
+# for use with sequential data
+# provides 8 blue colors of increasing saturation
+# compatible with gnuplot >=4.2
+# author: Anna Schneider
+
+# line styles
+set style line 1 lt 1 lc rgb '#F7FBFF' # very light blue
+set style line 2 lt 1 lc rgb '#DEEBF7' #
+set style line 3 lt 1 lc rgb '#C6DBEF' #
+set style line 4 lt 1 lc rgb '#9ECAE1' # light blue
+set style line 5 lt 1 lc rgb '#6BAED6' #
+set style line 6 lt 1 lc rgb '#4292C6' # medium blue
+set style line 7 lt 1 lc rgb '#2171B5' #
+set style line 8 lt 1 lc rgb '#084594' # dark blue
+
+# palette
+set palette defined ( 0 '#F7FBFF',\
+ 1 '#DEEBF7',\
+ 2 '#C6DBEF',\
+ 3 '#9ECAE1',\
+ 4 '#6BAED6',\
+ 5 '#4292C6',\
+ 6 '#2171B5',\
+ 7 '#084594' )
+
diff --git a/gnuplot/direction-dist.gnuplot b/gnuplot/direction-dist.gnuplot
@@ -0,0 +1,87 @@
+load 'gnuplot/blues.pal'
+set palette defined (0 '#4040c0', 1 '#c0c0c0', 2 '#c04040')
+set terminal svg size 450/2,450/2 font 'Free Serif, 10' enhanced round dynamic
+set xtics nomirror out offset 0,0.25
+set ytics nomirror out offset 0.5,0
+#set border 2+8
+unset border
+unset xtics
+unset ytics
+#set key top center outside Left reverse
+unset key
+
+#set xlabel 'υ_x'
+#set ylabel 'υ_y'
+set cbrange [0:*]
+set lmargin at screen 0.05
+set rmargin at screen 0.85
+set bmargin at screen 0.1
+set tmargin at screen 0.9
+
+positive(x) = x<0 ? (1/0) : x
+negative(x) = x>0 ? (1/0) : x
+round(x) = floor(x+0.5)
+length(x,y) = sqrt(x**2 + y**2)
+angle(x,y) = atan2(y,x)
+circle(x,y,z) = length(x,y)>1 ? NaN : z
+#circle(x,y,z) = z
+max(x,y) = x>y ? x : y
+
+filename = ARG1
+filename_tmp = sprintf('%s.tmp', filename)
+filename_cnt = sprintf('%s.cnt', filename)
+filename_svg = ARG2
+stats filename using 1:3
+x_max = max(abs(STATS_min_x),abs(STATS_max_x))
+y_max = max(abs(STATS_min_y),abs(STATS_max_y))
+print(y_max)
+
+# Min. direction RMSE = 0.024544 (1616497200), mean direction = -2.399689
+# Max. direction RMSE = 0.154791 (1616792400), mean direction = -2.440686
+# Mean direction RMSE = 0.056350ED6',\
+
+# http://www.gnuplotting.org/circular-heat-map/
+set table filename_tmp
+set dgrid3d 100,100
+set xrange [-1:1]
+set yrange [-1:1]
+set colorbox user origin 0.925,0.1 size 0.03,0.8
+splot filename using 1:2:3
+unset table
+
+set table filename_cnt
+set contour base
+set cntrparam levels 5
+set cntrparam order 8 bspline
+#set cntrlabel onecolor format '%8.3g' font ',7' start 25 interval -1
+unset surface
+splot filename using 1:2:3 with lines
+unset table
+
+set output filename_svg
+set multiplot
+plot filename_tmp using 1:2:(circle($1,$2,$3)) with image
+plot filename_cnt using (circle($1,$2,$1)):2 with lines lc '#404040'
+#plot 'build/1616792400-hist-xy.tmp' using 1:2 with points
+
+# now plot the polar grid only
+# https://stackoverflow.com/questions/18792461/gnuplot-2d-polar-plot-with-heatmap-from-3d-dataset-possible
+r = 3
+unset contour
+set border polar lw 5 lc rgb '#404040'
+set surface
+set angles degree
+set style line 11 lc rgb '#404040' lw 1 dashtype 2
+#set grid polar ls 11 nortics
+set polar
+set rrange[0:r]
+unset raxis
+set rtics format '' scale 0
+unset parametric
+set for [i=0:330:30] label at first (r+0.40)*cos(i), first (r+0.40)*sin(i)\
+center sprintf('%d', i)
+plot NaN w l
+unset multiplot
+
+
+#splot 'build/direction/1616497200-hist-xy' using 1:2:3 with pm3d
diff --git a/gnuplot/velocity-dist.gnuplot b/gnuplot/velocity-dist.gnuplot
@@ -1,5 +1,5 @@
load 'gnuplot/paired.pal'
-set terminal svg size 500,500/3*2 font 'Free Serif, 10' enhanced round dynamic
+set terminal svg size 450,450/3*2 font 'Free Serif, 10' enhanced round dynamic
set xtics nomirror out offset 0,0.25
set ytics nomirror out offset 0.5,0
set border 1+2 back
@@ -8,22 +8,6 @@ unset key
set output 'build/gnuplot/velocity-dist.svg'
-# 258 1615914000 0.08700399 0.06517805 0.06370091 1.660843
-# timestamp x_rmse y_rmse z_rmse speed
-# 295 1616187600 0.03527403 0.08648245 0.03217789 1.637611
-# timestamp x_rmse y_rmse z_rmse speed
-# 399 1616958000 0.02882994 0.03539276 0.1744659 0.8539361
-# timestamp x_rmse y_rmse z_rmse speed
-# 292 1616166000 0.06300805 0.02371221 0.094087 2.052675
-# [1] "MIN RMSE"
-# timestamp x_rmse y_rmse z_rmse speed
-# 322 1616382000 0.008555768 0.009697918 0.01404791 3.981209
-# timestamp x_rmse y_rmse z_rmse speed
-# 115 1614877200 0.01305041 0.009301931 0.02022811 4.356801
-# timestamp x_rmse y_rmse z_rmse speed
-# 39 1614330000 0.01867408 0.01648803 0.01133067 4.552804
-
-
set xlabel 'υ_x' offset 0,1.0
set yrange [0:1.2]
set rmargin 0
@@ -32,13 +16,22 @@ set multiplot layout 2,3
positive(x) = x<0 ? (1/0) : x
negative(x) = x>0 ? (1/0) : x
+round(x) = floor(x+0.5)
+
+# Velocity X max. RMSE = 0.123038 (1615230000), mean speed = 2.497133
+# Velocity Y max. RMSE = 0.115007 (1615388400), mean speed = 1.435405
+# Velocity Z max. RMSE = 0.141865 (1616958000), mean speed = 0.853936
+
+###################################
+# first row
+###################################
set lmargin 5
set title offset 0,-1
-set title sprintf('RMSE_x = %f', 0.08700399)
-set xrange [-4:4]
+set title sprintf('RMSE_x = %.0f%%', round(0.123038*100))
+set xrange [-5:3]
plot \
-'build/velocity/x/weibull/1615914000' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_x',\
+'build/velocity/x/weibull/1615230000' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_x',\
'' index 0 using (positive($1)):3 with lines lw 2 lc '#c04040' title 'Estimated PDF of positive υ_x',\
'' index 0 using (negative($1)):3 with lines lw 2 lc '#4040c0' title 'Estimated PDF of negative υ_x'
#'' index 1 using 1:2 with lines dashtype 1 lw 2 lc '#808080' notitle, \
@@ -47,39 +40,42 @@ plot \
set lmargin 5
#set xrange [-3:3]
set xlabel 'υ_y'
-set title sprintf('RMSE_y = %f', 0.08648245)
-set xrange [-3:5]
+set title sprintf('RMSE_y = %.0f%%',round(0.115007*100))
+set xrange [-4:4]
plot \
-'build/velocity/y/weibull/1616187600' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_y',\
+'build/velocity/y/weibull/1615388400' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_y',\
'' index 0 using (positive($1)):3 with lines lw 2 lc '#c04040' title 'Estimated PDF of positive υ_y',\
'' index 0 using (negative($1)):3 with lines lw 2 lc '#4040c0' title 'Estimated PDF of negative υ_y'
#'' index 1 using 1:2 with lines dashtype 1 lw 2 lc '#808080' notitle, \
#unset key
-set rmargin 2
set xrange [-4:4]
#set xtics -2,1,2
set xlabel 'υ_z'
-set title sprintf('RMSE_z = %f', 0.1744659)
+set title sprintf('RMSE_z = %.0f%%', round(0.141865*100))
plot \
'build/velocity/z/weibull/1616958000' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_z',\
'' index 0 using (positive($1)):3 with lines lw 2 lc '#c04040' title 'Estimated PDF of positive υ_z',\
'' index 0 using (negative($1)):3 with lines lw 2 lc '#4040c0' title 'Estimated PDF of negative υ_z'
#'' index 1 using 1:2 with lines dashtype 1 lw 2 lc '#808080' notitle, \
+# Velocity X min. RMSE = 0.016698 (1614135600), mean speed = 2.748603
+# Velocity Y min. RMSE = 0.021981 (1615647600), mean speed = 1.377537
+# Velocity Z min. RMSE = 0.015136 (1615647600), mean speed = 1.377537
+
###################################
# second row
###################################
set bmargin 3
-set xrange [-10:10]
-set yrange [0:0.4]
-set ytics 0,0.1
+set xrange [-4:4]
+set yrange [0:1.2]
+set ytics 0,0.2
#set xtics -10,1
set lmargin 5
-set title sprintf('RMSE_x = %f', 0.008555768)
+set title sprintf('RMSE_x = %.0f%%', round(0.016698*100))
plot \
-'build/velocity/x/weibull/1616382000' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_x',\
+'build/velocity/x/weibull/1614135600' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_x',\
'' index 0 using (positive($1)):3 with lines lw 2 lc '#c04040' title 'Estimated PDF of positive υ_x',\
'' index 0 using (negative($1)):3 with lines lw 2 lc '#4040c0' title 'Estimated PDF of negative υ_x'
#'' index 1 using 1:2 with lines dashtype 1 lw 2 lc '#808080' notitle, \
@@ -87,20 +83,19 @@ plot \
set lmargin 5
#set xrange [-3:3]
set xlabel 'υ_y'
-set title sprintf('RMSE_y = %f', 0.009301931)
+set title sprintf('RMSE_y = %.0f%%', round(0.021981*100))
#set xrange [-3:5]
plot \
-'build/velocity/y/weibull/1614877200' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_y',\
+'build/velocity/y/weibull/1615647600' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_y',\
'' index 0 using (positive($1)):3 with lines lw 2 lc '#c04040' title 'Estimated PDF of positive υ_y',\
'' index 0 using (negative($1)):3 with lines lw 2 lc '#4040c0' title 'Estimated PDF of negative υ_y'
#'' index 1 using 1:2 with lines dashtype 1 lw 2 lc '#808080' notitle, \
-set rmargin 2
#set xtics -2,1,2
set xlabel 'υ_z'
-set title sprintf('RMSE_z = %f', 0.01133067)
+set title sprintf('RMSE_z = %.0f%%', round(0.015136*100))
plot \
-'build/velocity/z/weibull/1614330000' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_z',\
+'build/velocity/z/weibull/1615647600' index 0 using 1:2 with points pt 6 lc '#404040' title 'Observed PDF of υ_z',\
'' index 0 using (positive($1)):3 with lines lw 2 lc '#c04040' title 'Estimated PDF of positive υ_z',\
'' index 0 using (negative($1)):3 with lines lw 2 lc '#4040c0' title 'Estimated PDF of negative υ_z'
#'' index 1 using 1:2 with lines dashtype 1 lw 2 lc '#808080' notitle, \
diff --git a/guix/manifest.scm b/guix/manifest.scm
@@ -6,6 +6,7 @@
(@ (gnu packages statistics) r-plotrix)
(@ (gnu packages statistics) r-foreach)
(@ (gnu packages statistics) r-doparallel)
+ (@ (gnu packages statistics) r-gplots)
(@ (gnu packages cran) r-fitdistrplus)
(@ (gnu packages statistics) r)
(@ (gnu packages maths) gnuplot)
diff --git a/main.tex b/main.tex
@@ -350,6 +350,20 @@ and high wind speeds. TODO
\cite{carta2008}
+\begin{figure}
+ \centering
+ \includegraphics{build/gnuplot/direction-dist-max.png}~
+ \includegraphics{build/gnuplot/direction-dist-min.png}
+ \caption{Per-axis wind velocity distributions fitted into Weibull
+ distribution.
+ Data intervals with the smallest error (the first row),
+ data intervals with the largest error (the second row). Red line shows estimated
+ probability density of positive wind speed projections, blue line shows estimated
+ probability density of negative wind speed projections and circles denote observed
+ probability density of wind speed projections.
+ \label{fig-direction-distributions}}
+\end{figure}
+
\section{Discussion}
One disadvantage of three-axis anemometer is that the arm for the \(z\) axis is
