commit 27bc5fc97c7fc11e28920896e575238285121b40
parent 171827e17c4cad8461de940d9759bb0895fe9ffe
Author: Ivan Gankevich <igankevich@ya.ru>
Date: Sat, 14 Jan 2017 15:46:23 +0300
Skew normal distribution graphs in R.
Diffstat:
| R/common.R | | | 35 | +++++++++++++++++++++++++++++++++++ |
| R/transform.R | | | 23 | +++++++++++++++++++++++ |
| phd-diss-ru.org | | | 60 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-- |
| phd-diss.org | | | 60 | ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-- |
4 files changed, 174 insertions(+), 4 deletions(-)
diff --git a/R/common.R b/R/common.R
@@ -1,4 +1,5 @@
source(file.path("R", "waves.R"))
+source(file.path("R", "transform.R"))
arma.qqplot_grid <- function (dir, params) {
wave_params <- arma.load_wave_parameters(dir, params)
@@ -30,3 +31,37 @@ arma.wavy_plot <- function (data, t, ...) {
facetcol <- cut(zfacet, nbcol)
persp(x, y, z, phi=30, theta=30, col=color[facetcol], ...)
}
+
+arma.skew_normal_1_plot <- function(x, params) {
+ data <- mapply(
+ function (s, k) arma.skew_normal_1(x, s, k),
+ params$skewness,
+ params$kurtosis
+ )
+ plot.new()
+ plot.window(xlim=range(x),ylim=range(data))
+ axis(1); axis(2); box()
+ for (i in seq_len(ncol(data))) {
+ d <- data[,i]
+ lines(x, d, lty=paste(params$linetypes[[i]]))
+ }
+}
+
+
+arma.skew_normal_2_plot <- function(x, params) {
+ data <- mapply(
+ function (a) arma.skew_normal_2(x, a),
+ params$alpha
+ )
+ plot.new()
+ plot.window(xlim=range(x),ylim=range(data))
+ axis(1); axis(2); box()
+ for (i in seq_len(ncol(data))) {
+ d <- data[,i]
+ lines(x, d, lty=paste(params$linetypes[[i]]))
+ }
+}
+
+arma.fmt <- function(x, ndigits) {
+ format(round(x, ndigits), nsmall=ndigits)
+}
diff --git a/R/transform.R b/R/transform.R
@@ -0,0 +1,23 @@
+arma.skew_normal_1 <- function(z, skewness, kurtosis) {
+ (24 + 4*z*(-3 + z**2)* skewness + (3 - 6*z**2 + z**4)* kurtosis)/(24*(exp((z**2)/2)*sqrt(2*pi)))
+}
+
+arma.bits.erf <- function(x) {
+ 2 * pnorm(x * sqrt(2)) - 1
+}
+
+arma.bits.erfc <- function(x) {
+ 2 * pnorm(x * sqrt(2), lower = FALSE)
+}
+
+arma.bits.skewness_2 <- function(alpha) {
+ (sqrt(2) * (4 - pi) * (alpha**3)) / (sqrt(pi + (pi-2) * (alpha**2))**3)
+}
+
+arma.bits.kurtosis_2 <- function(alpha) {
+ (8*(-3 + pi)*(alpha**4))/ ((pi + (-2 + pi)*(alpha**2))**2)
+}
+
+arma.skew_normal_2 <- function(z, alpha) {
+ arma.bits.erfc( -((z*alpha)/sqrt(2)) ) / (exp((z**2)/2)*sqrt(2*pi))
+}
diff --git a/phd-diss-ru.org b/phd-diss-ru.org
@@ -1288,8 +1288,34 @@ $\gamma_2$ --- эксцесс, $f$ --- ФПР, $F$ --- функция распр
распределения при различных параметрах показано на [[fig:skew-normal-1]].
#+name: fig:skew-normal-1
+#+begin_src R :results output graphics :exports results :file build/skew-normal-1.pdf
+source(file.path("R", "common.R"))
+x <- seq(-3, 3, length.out=100)
+params <- data.frame(
+ skewness = c(0.00, 0.52, 0.00, 0.52),
+ kurtosis = c(0.00, 0.00, 0.70, 0.70),
+ linetypes = c("solid", "dashed", "dotdash", "dotted")
+)
+arma.skew_normal_1_plot(x, params)
+legend(
+ "topleft",
+ mapply(
+ function (s, k) {
+ as.expression(bquote(list(
+ gamma[1] == .(arma.fmt(s, 2)),
+ gamma[2] == .(arma.fmt(k, 2))
+ )))
+ },
+ params$skewness,
+ params$kurtosis
+ ),
+ lty = paste(params$linetypes)
+)
+#+end_src
+
#+caption: Вид плотности распределения eqref:eq:skew-normal-1 аппликат взволнованной морской поверхности при различных значениях асимметрии $\gamma_1$ и эксцесса $\gamma_2$.
-[[file:build/skew-normal-1.eps]]
+#+RESULTS: fig:skew-normal-1
+[[file:build/skew-normal-1.pdf]]
**** Асимметричное нормальное распределение.
Альтернативной аппроксимацией распределения волновых аппликат служит формула
@@ -1309,8 +1335,38 @@ $\gamma_2$ --- эксцесс, $f$ --- ФПР, $F$ --- функция распр
на [[fig:skew-normal-2]].
#+name: fig:skew-normal-2
+#+begin_src R :results output graphics :exports results :file build/skew-normal-2.pdf
+source(file.path("R", "common.R"))
+x <- seq(-3, 3, length.out=100)
+alpha <- c(0.00, 0.87, 2.25, 4.90)
+params <- data.frame(
+ alpha = alpha,
+ skewness = arma.bits.skewness_2(alpha),
+ kurtosis = arma.bits.kurtosis_2(alpha),
+ linetypes = c("solid", "dashed", "dotdash", "dotted")
+)
+arma.skew_normal_2_plot(x, params)
+legend(
+ "topleft",
+ mapply(
+ function (a, s, k) {
+ as.expression(bquote(list(
+ alpha == .(arma.fmt(a, 2)),
+ gamma[1] == .(arma.fmt(s, 2)),
+ gamma[2] == .(arma.fmt(k, 2))
+ )))
+ },
+ params$alpha,
+ params$skewness,
+ params$kurtosis
+ ),
+ lty = paste(params$linetypes)
+)
+#+end_src
+
#+caption: Вид плотности распределения eqref:eq:skew-normal-2 волновых аппликат при различных значениях коэффициента асимметрии $\alpha$.
-[[file:build/skew-normal-2.eps]]
+#+RESULTS: fig:skew-normal-2
+[[file:build/skew-normal-2.pdf]]
**** Тестирование.
Решение уравнения eqref:eq:distribution-transformation с выбранной функцией
diff --git a/phd-diss.org b/phd-diss.org
@@ -1375,8 +1375,34 @@ $0.1\leq\gamma_2\leq{0.7}$. Family of probability density functions for
different parameters is shown in [[fig:skew-normal-1]].
#+name: fig:skew-normal-1
+#+begin_src R :results output graphics :exports results :file build/skew-normal-1.pdf
+source(file.path("R", "common.R"))
+x <- seq(-3, 3, length.out=100)
+params <- data.frame(
+ skewness = c(0.00, 0.52, 0.00, 0.52),
+ kurtosis = c(0.00, 0.00, 0.70, 0.70),
+ linetypes = c("solid", "dashed", "dotdash", "dotted")
+)
+arma.skew_normal_1_plot(x, params)
+legend(
+ "topleft",
+ mapply(
+ function (s, k) {
+ as.expression(bquote(list(
+ gamma[1] == .(arma.fmt(s, 2)),
+ gamma[2] == .(arma.fmt(k, 2))
+ )))
+ },
+ params$skewness,
+ params$kurtosis
+ ),
+ lty = paste(params$linetypes)
+)
+#+end_src
+
#+caption: Probability density function eqref:eq:skew-normal-1 of ocean wavy surface elevation for different values of skewness $\gamma_1$ and kurtosis $\gamma_2$.
-[[file:build/skew-normal-1.eps]]
+#+RESULTS: fig:skew-normal-1
+[[file:build/skew-normal-1.pdf]]
**** Skew-normal distribution.
Alternative approach is to approximate distribution of ocean wavy surface
@@ -1396,8 +1422,38 @@ and mathematical libraries. Its graph for different values of $\alpha$ is shown
in [[fig:skew-normal-2]].
#+name: fig:skew-normal-2
+#+begin_src R :results output graphics :exports results :file build/skew-normal-2.pdf
+source(file.path("R", "common.R"))
+x <- seq(-3, 3, length.out=100)
+alpha <- c(0.00, 0.87, 2.25, 4.90)
+params <- data.frame(
+ alpha = alpha,
+ skewness = arma.bits.skewness_2(alpha),
+ kurtosis = arma.bits.kurtosis_2(alpha),
+ linetypes = c("solid", "dashed", "dotdash", "dotted")
+)
+arma.skew_normal_2_plot(x, params)
+legend(
+ "topleft",
+ mapply(
+ function (a, s, k) {
+ as.expression(bquote(list(
+ alpha == .(arma.fmt(a, 2)),
+ gamma[1] == .(arma.fmt(s, 2)),
+ gamma[2] == .(arma.fmt(k, 2))
+ )))
+ },
+ params$alpha,
+ params$skewness,
+ params$kurtosis
+ ),
+ lty = paste(params$linetypes)
+)
+#+end_src
+
#+caption: Probability density function eqref:eq:skew-normal-2 of ocean wavy surface for different values of skewness coefficient $\alpha$.
-[[file:build/skew-normal-2.eps]]
+#+RESULTS: fig:skew-normal-2
+[[file:build/skew-normal-2.pdf]]
**** Evaluation.
Equation eqref:eq:distribution-transformation with selected wave elevation
