From cff88c326cb7d222016f1bc7ac6478f903a42a0b Mon Sep 17 00:00:00 2001 From: Mateusz Staniak Date: Wed, 6 Jun 2018 10:17:55 +0000 Subject: [PATCH] version 0.9.2 --- DESCRIPTION | 23 +++ LICENSE | 2 + MD5 | 57 ++++++ NAMESPACE | 35 ++++ NEWS.md | 21 +++ R/barrier_crossing.R | 125 +++++++++++++ R/exact_ladder_moments.R | 86 +++++++++ R/kendallRandomPackage.R | 29 +++ R/simulations_kendall_rw.R | 256 ++++++++++++++++++++++++++ R/stable_kendall_distribution.R | 172 +++++++++++++++++ R/stable_kendall_fit.R | 198 ++++++++++++++++++++ README.md | 35 ++++ build/vignette.rds | Bin 0 -> 248 bytes inst/doc/behaviour.R | 54 ++++++ inst/doc/behaviour.Rmd | 81 ++++++++ inst/doc/behaviour.html | 152 +++++++++++++++ inst/doc/kendall_rws.R | 41 +++++ inst/doc/kendall_rws.Rmd | 67 +++++++ inst/doc/kendall_rws.html | 164 +++++++++++++++++ man/Qn.Rd | 21 +++ man/U.Rd | 19 ++ man/Z.Rd | 19 ++ man/dkend.Rd | 25 +++ man/estimate_stable_alpha.Rd | 17 ++ man/fit_kendall.Rd | 17 ++ man/fit_separate.Rd | 20 ++ man/fit_stable_alpha.Rd | 17 ++ man/full_loglik_gradient.Rd | 17 ++ man/full_minus_loglik.Rd | 17 ++ man/g_function.Rd | 27 +++ man/g_function_single.Rd | 26 +++ man/kendallRandomWalks.Rd | 33 ++++ man/kendall_loglik.Rd | 19 ++ man/ladder_height.Rd | 30 +++ man/ladder_moment.Rd | 30 +++ man/ladder_moment_pmf.Rd | 31 ++++ man/mutate_kendall_rw.Rd | 22 +++ man/pkend.Rd | 25 +++ man/pkendSym.Rd | 25 +++ man/plot.kendall_barrier_crossing.Rd | 19 ++ man/plot.kendall_fit.Rd | 19 ++ man/plot.kendall_simulation.Rd | 26 +++ man/print.kendall_barrier_crossing.Rd | 19 ++ man/print.kendall_simulation.Rd | 16 ++ man/qkend.Rd | 25 +++ man/qkendSym.Rd | 25 +++ man/rkend.Rd | 25 +++ man/simulateOneTrajectory.Rd | 26 +++ man/simulate_kendall_rw.Rd | 40 ++++ man/summarise_kendall_rw.Rd | 19 ++ man/transform_kendall_rw.Rd | 30 +++ tests/testthat.R | 4 + tests/testthat/test_barriers.R | 63 +++++++ tests/testthat/test_dists.R | 48 +++++ tests/testthat/test_fit.R | 38 ++++ tests/testthat/test_simulations.R | 44 +++++ vignettes/behaviour.Rmd | 81 ++++++++ vignettes/kendall_rws.Rmd | 67 +++++++ 58 files changed, 2689 insertions(+) create mode 100644 DESCRIPTION create mode 100644 LICENSE create mode 100644 MD5 create mode 100644 NAMESPACE create mode 100644 NEWS.md create mode 100644 R/barrier_crossing.R create mode 100644 R/exact_ladder_moments.R create mode 100644 R/kendallRandomPackage.R create mode 100644 R/simulations_kendall_rw.R create mode 100644 R/stable_kendall_distribution.R create mode 100644 R/stable_kendall_fit.R create mode 100644 README.md create mode 100644 build/vignette.rds create mode 100644 inst/doc/behaviour.R create mode 100644 inst/doc/behaviour.Rmd create mode 100644 inst/doc/behaviour.html create mode 100644 inst/doc/kendall_rws.R create mode 100644 inst/doc/kendall_rws.Rmd create mode 100644 inst/doc/kendall_rws.html create mode 100644 man/Qn.Rd create mode 100644 man/U.Rd create mode 100644 man/Z.Rd create mode 100644 man/dkend.Rd create mode 100644 man/estimate_stable_alpha.Rd create mode 100644 man/fit_kendall.Rd create mode 100644 man/fit_separate.Rd create mode 100644 man/fit_stable_alpha.Rd create mode 100644 man/full_loglik_gradient.Rd create mode 100644 man/full_minus_loglik.Rd create mode 100644 man/g_function.Rd create mode 100644 man/g_function_single.Rd create mode 100644 man/kendallRandomWalks.Rd create mode 100644 man/kendall_loglik.Rd create mode 100644 man/ladder_height.Rd create mode 100644 man/ladder_moment.Rd create mode 100644 man/ladder_moment_pmf.Rd create mode 100644 man/mutate_kendall_rw.Rd create mode 100644 man/pkend.Rd create mode 100644 man/pkendSym.Rd create mode 100644 man/plot.kendall_barrier_crossing.Rd create mode 100644 man/plot.kendall_fit.Rd create mode 100644 man/plot.kendall_simulation.Rd create mode 100644 man/print.kendall_barrier_crossing.Rd create mode 100644 man/print.kendall_simulation.Rd create mode 100644 man/qkend.Rd create mode 100644 man/qkendSym.Rd create mode 100644 man/rkend.Rd create mode 100644 man/simulateOneTrajectory.Rd create mode 100644 man/simulate_kendall_rw.Rd create mode 100644 man/summarise_kendall_rw.Rd create mode 100644 man/transform_kendall_rw.Rd create mode 100644 tests/testthat.R create mode 100644 tests/testthat/test_barriers.R create mode 100644 tests/testthat/test_dists.R create mode 100644 tests/testthat/test_fit.R create mode 100644 tests/testthat/test_simulations.R create mode 100644 vignettes/behaviour.Rmd create mode 100644 vignettes/kendall_rws.Rmd diff --git a/DESCRIPTION b/DESCRIPTION new file mode 100644 index 0000000..7cc8ea9 --- /dev/null +++ b/DESCRIPTION @@ -0,0 +1,23 @@ +Package: kendallRandomWalks +Title: Simulate and Visualize Kendall Random Walks and Related + Distributions +Version: 0.9.2 +Authors@R: person("Mateusz", "Staniak", email = "mateusz.staniak@math.uni.wroc.pl", role = c("aut", "cre")) +Description: Kendall random walks are a continuous-space Markov chains generated + by the Kendall generalized convolution. This package provides tools + for simulating these random walks and studying distributions + related to them. For more information about Kendall random walks see Jasiulis-Gołdyn (2014) . +Depends: R (>= 3.3) +License: MIT + file LICENSE +Encoding: UTF-8 +LazyData: true +Imports: ggplot2, dplyr, tidyr, EnvStats, tibble, nleqslv +RoxygenNote: 6.0.1 +Suggests: knitr, rmarkdown, testthat, covr +VignetteBuilder: knitr +NeedsCompilation: no +Packaged: 2018-06-05 18:33:35 UTC; mtst +Author: Mateusz Staniak [aut, cre] +Maintainer: Mateusz Staniak +Repository: CRAN +Date/Publication: 2018-06-06 11:17:55 UTC diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..ea94ae4 --- /dev/null +++ b/LICENSE @@ -0,0 +1,2 @@ +YEAR: 2018 +COPYRIGHT HOLDER: Mateusz Staniak diff --git a/MD5 b/MD5 new file mode 100644 index 0000000..c3c7ece --- /dev/null +++ b/MD5 @@ -0,0 +1,57 @@ +802161c45338e0f65405afc5c1e31b88 *DESCRIPTION +127cbfaba59733e6a2c6a1cb068efef0 *LICENSE +bf40ddb040e1475489436cc818624dbd *NAMESPACE +d2b472e46a573de5a3285e5a9ddcc53c *NEWS.md +321cb795d29c29dd544278a672d66d3a *R/barrier_crossing.R +44742a93c96bea8c1e66ebddf3ee3b11 *R/exact_ladder_moments.R +0c0ea4b64e63a12cf1c77d1e192a3586 *R/kendallRandomPackage.R +342bab9561ad2b935f5ef3d3363b1951 *R/simulations_kendall_rw.R +0cdf1ddbbe6db456b3e2761647dbc7cb *R/stable_kendall_distribution.R +ecc63136d6a466f61e1ed491dad20450 *R/stable_kendall_fit.R +b6b60ec945f8a1a226a2bcb0d375a493 *README.md +b84cc913d6ce0d3aaec243382c8dfe9d *build/vignette.rds +a48fdecfe2262ed77d8b61e89c783509 *inst/doc/behaviour.R +7161eac5f4c752b0251bead321523c98 *inst/doc/behaviour.Rmd +95472fdbb4e06f604fc52498a3335f38 *inst/doc/behaviour.html +57b6d03ded5319bbee71626929038290 *inst/doc/kendall_rws.R +a2ba97d109c7e774099396657baf7dee *inst/doc/kendall_rws.Rmd +11a135a403e922a11064872a6a20344e *inst/doc/kendall_rws.html +90cfad0d498cc81140a99dcad617506f *man/Qn.Rd +bac14a58926029a0fea03be2cf7e4b67 *man/U.Rd +dce1d5e75d456212d0000537beb067d8 *man/Z.Rd +1bcbad8ce48ea83faf4981ed33c5f6a9 *man/dkend.Rd +7296860b96d294cfa90a5ddb143e2911 *man/estimate_stable_alpha.Rd +43a4e3c0349028a7207ff90add1f1ed4 *man/fit_kendall.Rd +daf8f31ad99cf97685028d8a8e0dec70 *man/fit_separate.Rd +1e078f529dd4b789c0c4e20b0652a1c8 *man/fit_stable_alpha.Rd +dc478a9804d04e57d5bcbc1b479ce038 *man/full_loglik_gradient.Rd +e966dbd9c7ff22c98413cc5fe52e73f0 *man/full_minus_loglik.Rd +846d3b0a92c23421e9ed338b9f2aa1a3 *man/g_function.Rd +f0682dab3c6e09294d0cd304504fe201 *man/g_function_single.Rd +ee19dae30b9c959f00384f52fd5d4fc5 *man/kendallRandomWalks.Rd +e75c117dd539f486a67685b7f4d8d4d9 *man/kendall_loglik.Rd +741d52e764454a56e6256747236e1220 *man/ladder_height.Rd +00a5d7bd2e162f4dfcde322d8c790fdc *man/ladder_moment.Rd +2a28176cce226256b3ebfb0d38285036 *man/ladder_moment_pmf.Rd +3387928db5e95dac979483e734955cde *man/mutate_kendall_rw.Rd +7266e4a1ff17b60c0ba629a1fe93a026 *man/pkend.Rd +472d9b2fdd3a2ae9e110c8b24924c59f *man/pkendSym.Rd +654edc804cc98fa738f73cd7c25daf28 *man/plot.kendall_barrier_crossing.Rd +f561f98fdea9e12c299cbc836b32d2f1 *man/plot.kendall_fit.Rd +490a4121b7c9eb59a3cf4a46e536288d *man/plot.kendall_simulation.Rd +a89f67567cf8a3c7e8066700c8fb2e5a *man/print.kendall_barrier_crossing.Rd +8f0bcbc88a8135114238114a78888895 *man/print.kendall_simulation.Rd +4f1a7c340be93276c09ad4f942643629 *man/qkend.Rd +1ce650196a895a90eaf4bf0dfe2ca00b *man/qkendSym.Rd +ed3c5ae983d3eb2dfc2683d8757ae478 *man/rkend.Rd +56d91d1a29b9af683a57f03ed9a08fe1 *man/simulateOneTrajectory.Rd +3154f95637a38863cb5786ae6da69b3d *man/simulate_kendall_rw.Rd +d8c7c8c7a35a36d36178bae6c9f7aa90 *man/summarise_kendall_rw.Rd +867c1a92807ac63f896b3832bc347a55 *man/transform_kendall_rw.Rd +1ccf23728997c0389f54ddc50a4ed5da *tests/testthat.R +467ab774ad60f93ebb7667bd63db79fc *tests/testthat/test_barriers.R +a48ae04e756695fa1c0c30d42c222541 *tests/testthat/test_dists.R +4a31b49bed12cd9b0b7ef16aabab43e6 *tests/testthat/test_fit.R +991a9cfa7116d3b333f6289e2e971dee *tests/testthat/test_simulations.R +7161eac5f4c752b0251bead321523c98 *vignettes/behaviour.Rmd +a2ba97d109c7e774099396657baf7dee *vignettes/kendall_rws.Rmd diff --git a/NAMESPACE b/NAMESPACE new file mode 100644 index 0000000..05c47c0 --- /dev/null +++ b/NAMESPACE @@ -0,0 +1,35 @@ +# Generated by roxygen2: do not edit by hand + +S3method(plot,kendall_barrier_crossing) +S3method(plot,kendall_fit) +S3method(plot,kendall_simulation) +S3method(print,kendall_barrier_crossing) +S3method(print,kendall_simulation) +export(dkend) +export(estimate_stable_alpha) +export(fit_kendall) +export(fit_separate) +export(fit_stable_alpha) +export(full_loglik_gradient) +export(full_minus_loglik) +export(g_function) +export(g_function_single) +export(ladder_height) +export(ladder_moment) +export(ladder_moment_pmf) +export(mutate_kendall_rw) +export(pkend) +export(pkendSym) +export(qkend) +export(qkendSym) +export(rkend) +export(simulate_kendall_rw) +export(summarise_kendall_rw) +export(transform_kendall_rw) +importFrom(dplyr,n) +importFrom(stats,integrate) +importFrom(stats,optim) +importFrom(stats,optimize) +importFrom(stats,quantile) +importFrom(stats,rgamma) +importFrom(stats,sd) diff --git a/NEWS.md b/NEWS.md new file mode 100644 index 0000000..f02d52e --- /dev/null +++ b/NEWS.md @@ -0,0 +1,21 @@ +# kendallRandomPackage 0.9.2 + +* Functions for transformations and summaries of simulations of Kendall random walks + (`summarise_kendall_rw` and `mutate_kendall_rw`). + +# kendallRandomPackage 0.9.1 + +* Numerical computation of G function (so also Williamson transform) +* Exact distribution of first ladder moments + +# kendallRandomPackage 0.9.0 + +* Generic function for printing and plotting simulation objects +* Barrier crossing simulations and visualizations + +# kendallRandomPackage 0.0.0.9000 + +* Added a `NEWS.md` file to track changes to the package. + + + diff --git a/R/barrier_crossing.R b/R/barrier_crossing.R new file mode 100644 index 0000000..50e5a60 --- /dev/null +++ b/R/barrier_crossing.R @@ -0,0 +1,125 @@ +#' Estimate the distribution of first ladder moment for given level +#' +#' NA is returned if the level wasn't crossed. Printing the resulting object +#' will give summary of the estimated distribution and information whether +#' level wasn't crossed in some simulations. This information can be used to +#' pick the right trajectory length for the given level. +#' +#' @param simulations kendall_simulation object +#' @param level Positive numeric +#' +#' @return tibble +#' +#' @importFrom dplyr n +#' +#' @export +#' +#' @examples { +#' kendall_rw <- simulate_kendall_rw(100, 100, runif, 0.5) +#' estim_ladder <- ladder_moment(kendall_rw, 1000) +#' estim_ladder +#' } +#' + +ladder_moment <- function(simulations, level) { + if(level < 0) stop("Level must be positive") + sim <- id <- sim_id <- NULL + kendall_rw <- simulations$simulation + kendall_rw <- dplyr::group_by(kendall_rw, sim_id) + kendall_rw <- dplyr::mutate(kendall_rw, id = 1:n()) + all_sims <- dplyr::ungroup(dplyr::distinct(kendall_rw, sim_id)) + kendall_rw <- dplyr::filter(kendall_rw, sim > level) + kendall_rw <- dplyr::summarise(kendall_rw, ladder_moment = min(id)) + kendall_rw <- dplyr::ungroup(kendall_rw) + all_sims <- dplyr::left_join(all_sims, kendall_rw, by = "sim_id") + + class(all_sims) <- c("kendall_barrier_crossing", class(all_sims)) + all_sims +} + + +#' Estimate the distribution of first ladder height for given level +#' +#' NA is returned if the level wasn't crossed. Printing the resulting object +#' will give summary of the estimated distribution and information whether +#' level wasn't crossed in some simulations. This information can be used to +#' pick the right trajectory length for the given level. +#' +#' @param simulations kendall_simulation object +#' @param level Positive numeric +#' +#' @return tibble +#' +#' @importFrom dplyr n +#' +#' @export +#' +#' @examples { +#' kendall_rw <- simulate_kendall_rw(100, 100, runif, 0.5) +#' estim_ladder <- ladder_height(kendall_rw, 1000) +#' estim_ladder +#' } +#' + +ladder_height <- function(simulations, level) { + if(level < 0) stop("Level must be positive") + sim <- id <- sim_id <- NULL + kendall_rw <- simulations$simulation + kendall_rw <- dplyr::group_by(kendall_rw, sim_id) + kendall_rw <- dplyr::mutate(kendall_rw, id = 1:n()) + all_sims <- dplyr::ungroup(dplyr::distinct(kendall_rw, sim_id)) + kendall_rw <- dplyr::filter(kendall_rw, sim > level) + kendall_rw <- dplyr::summarise(kendall_rw, ladder_moment = min(sim)) + kendall_rw <- dplyr::ungroup(kendall_rw) + all_sims <- dplyr::left_join(all_sims, kendall_rw, by = "sim_id") + + class(all_sims) <- c("kendall_barrier_crossing", class(all_sims)) + all_sims +} + +#' Generic function for printing result of ladder_moment function +#' +#' @param x kendall_barrier_crossing object +#' @param ... Additional arguments +#' +#' @return invisible x +#' +#' @importFrom stats quantile sd +#' +#' @export +#' + +print.kendall_barrier_crossing <- function(x, ...) { + quantiles <- quantile(x$ladder_moment, na.rm = T, + probs = seq(0, 1, by = 0.1)) + labels <- names(quantiles) + cat("Mean of the distribution: ", mean(x$ladder_moment, na.rm = T), "\n") + cat("Standard deviation of the distribution: ", sd(x$ladder_moment, na.rm = T), "\n") + cat("Number of observations: ", max(x$sim_id), "\n") + cat("Times the level was not crossed: ", sum(!is.finite(x$ladder_moment)), "\n") + cat("Quantiles of the distribution: \n") + print(quantiles) + invisible(x) +} + + +#' Generic function for plotting results of ladder_moment function. +#' +#' @param x kendall_barrier_crossing object +#' @param ... Additional arguments +#' +#' @return ggplot2 +#' +#' @export +#' + +plot.kendall_barrier_crossing <- function(x, ...) { + mean_value <- mean(x$ladder_moment, na.rm = TRUE) + ggplot2::ggplot(x, ggplot2::aes_string(x = 'ladder_moment')) + + ggplot2::geom_histogram() + + ggplot2::geom_vline(xintercept = mean_value) + + ggplot2::theme_bw() + + ggplot2::xlab("First ladder moments") + + ggplot2::ylab("Count") + # dodac tu ecdf +} diff --git a/R/exact_ladder_moments.R b/R/exact_ladder_moments.R new file mode 100644 index 0000000..bce398d --- /dev/null +++ b/R/exact_ladder_moments.R @@ -0,0 +1,86 @@ +#' Function G(t) - Williamson transform taken at point 1/t. +#' +#' This function return the whole "integrate" object, so precision of the approximation +#' can be checked. +#' +#' @param t Argument to the function. +#' @param alpha Value of the alpha parameter. +#' @param density Density function of the step distribution. +#' +#' @return Object of class "integrate" +#' +#' @importFrom stats integrate +#' +#' @examples +#' g_function_single(5, 0.26, dnorm) +#' +#' @export +#' + +g_function_single <- function(t, alpha, density) { + force(t) + under_integral <- function(x) { + (1 - abs(x/t)^(alpha))*density(x) + } + if(t != 0) { + integrate(under_integral, lower = -abs(t), upper = abs(t)) + } else { + fun <- function(x) x + integrate(fun, lower = -1, upper = 1) + } +} + + +#' Function G(t) - Williamson transform taken at point 1/t. +#' +#' This function return only approximated values. To check their precisions use +#' g_function_single function with an argument of length 1. +#' +#' @param t Argument to the function. +#' @param alpha Value of the alpha parameter. +#' @param density Density function of the step distribution. +#' +#' @return Object of class "integrate" +#' +#' @importFrom stats integrate +#' +#' @export +#' +#' @examples +#' g_function(1:5, 0.75, dnorm) +#' +#' + +g_function <- function(t, alpha, density) { + sapply(t, function(x) g_function_single(x, alpha, density)$value) +} + + +#' Distribution of the first ladder moment. +#' +#' @param n Argument to the PDF. +#' @param level Level a to be crossed. +#' @param alpha Alpha parameter of Kendall random walk. +#' @param step_cdf CDF of the step distribution. +#' @param step_pdf PDF of the step distribution. +#' +#' @return Value of PMF of the distribution of first ladder moment +#' +#' @export +#' +#' @examples +#' prob <- ladder_moment_pmf(10, 1000, 0.5, pnorm, dnorm) +#' prob +#' +#' + +ladder_moment_pmf <- function(n, level, alpha, step_cdf, step_pdf) { + Ga <- g_function(level, alpha, step_pdf) + Fa <- step_cdf(level) + Ha <- 2*Fa - 1 - Ga + A <- 1 + Ha/((2*Ga - 1)^2) - Ga/(2*Ga - 1) + B <- Ha/((1 - Ga)*(2*Ga - 1)) + C <- Ga/(2*Ga - 1) - (Ha*Ga)/(((2*Ga - 1)^2)*(1 - Ga)) + + A*(2^(-n)) + B*n*((1 - Ga)^2)*(Ga^(n - 1)) + C*(1 - Ga)*(Ga^(n - 1)) +} diff --git a/R/kendallRandomPackage.R b/R/kendallRandomPackage.R new file mode 100644 index 0000000..55ce818 --- /dev/null +++ b/R/kendallRandomPackage.R @@ -0,0 +1,29 @@ +#' kendallRandomWalks: explore and visualize Kendall random walks. +#' +#' Kendall random walks are Markov processes generated by the Kendall convolution. +#' This package helps simulate and visualize these random walks and associated +#' distributions. It also provides function to fit these distributions to data. +#' +#' @section Important functions: +#' \code{\link{simulate_kendall_rw}} simulates Kendall random walks. +#' +#' \code{\link{transform_kendall_rw}} applies scaling and shift to simulated Kendall r.w-s. +#' +#' \code{\link{ladder_moment}} estimates the distribution of first ladder moment +#' by simulating Kendall random walks. +#' +#' \code{\link{ladder_height}} estimates the distribution of first ladder height +#' by simulating Kendall random walks. +#' +#' \code{\link{ladder_moment_pmf}} computes the PMF of the distribution of first +#' ladder moment. +#' +#' \code{\link{g_function}} Finds the value of G(t) numerically. +#' +#' \code{\link{pkend}}, \code{\link{dkend}}, \code{\link{qkend}}, \code{\link{rkend}} +#' give CDF, PDF, quantile function and random numbers from stable Kendall distribution. +#' +#' +#' @docType package +#' @name kendallRandomWalks +NULL diff --git a/R/simulations_kendall_rw.R b/R/simulations_kendall_rw.R new file mode 100644 index 0000000..b8d3772 --- /dev/null +++ b/R/simulations_kendall_rw.R @@ -0,0 +1,256 @@ +#' Helper function: min/max +#' +#' @param x numeric +#' @param y numeric +#' +#' @return min of arguments divided by max of arguments +#' + +Z <- function(x, y){ + min(x, y)/max(x, y) +} + +#' Helper function +#' +#' @param x numeric +#' @param y numeric +#' @param alpha numeric, parameter of Kendall random walk +#' +#' @return 0 or 1 with probability depending on x, y, alpha +#' + +Qn <- function(x, y, alpha){ + p <- Z(abs(x), abs(y))^alpha + sample(c(0,1), 1, prob=c(1-p, p)) +} + +#' Helper function +#' +#' @param x numeric +#' @param y numeric +#' +#' @return sign of the argument whose abs. val. is bigger +#' + +U <- function(x, y){ + if(abs(x) > abs(y)) sign(x) + else sign(y) +} + +#' Simulate one trajectory ofa Kendall random walk +#' +#' @param trajectory_length Number of samples to simulate. +#' @param step_dist Function that returns random numbers from step distribution. +#' @param alpha Alpha parameter of the random walk +#' @param symmetric If TRUE, random walk on the whole real line will be simulated. +#' @param ... Additional parameters to step distribution. +#' +#' @return Generated path of the random walk. +#' +#' + +simulateOneTrajectory <- function(trajectory_length, step_dist, + alpha, symmetric = FALSE, ...) { + Y <- step_dist(trajectory_length, ...) + if(symmetric) { + theta <- EnvStats::rpareto(trajectory_length, + 1, + 2*alpha)*sample(c(-1, 1), + trajectory_length, + prob = c(0.5, 0.5), + replace = TRUE) + } else { + theta <- EnvStats::rpareto(trajectory_length, 1, 2*alpha) + } + + Xn <- vector("numeric", trajectory_length ) + Xn[1:2] <- c(0, Y[1]) + + for(i in 2:(trajectory_length - 1)) { + Xn[i + 1] <- max(abs(Xn[i]), + abs(Y[i + 1]))*theta[i]^Qn(Xn[i], + Y[i+1], + alpha)*U(Xn[i], + Y[i + 1]) + } + Xn +} + +#' Simulate multiple trajectories of Kendall random walk +#' +#' Object returned by this has print and plot methods. +#' +#' @param number_of_simulations number of trajectories to generate. +#' @param trajectory_length length of trajectories. +#' @param step_dist function returning random numbers from step dist. +#' @param alpha alpha parameter. +#' @param symmetric If TRUE, random walk on the whole real line will be simulated. +#' @param ... parameters for step distribution. +#' +#' @return Object of class kendall_simulation. It is a list that consists of +#' \item{simulation}{Tibble with simulation id and simulated values,} +#' \item{step_distribution}{Name of the step distribution,} +#' \item{alpha}{Value of alpha parameter,} +#' \item{is_symmetric}{Logical value indicating if this is a symmetric Kendall R.W.} +#' +#' @export +#' +#' @examples +#' kendall_simulations <- simulate_kendall_rw(10, 1000, runif, 0.5) +#' # Kendall R.W. on positive half-line with uniform step distribution - 10 trajectories. +#' only_simulations <- kendall_simulations$simulation # tibble with simulated values +#' kendall_simulations +#' +#' + +simulate_kendall_rw <- function(number_of_simulations, trajectory_length, + step_dist, alpha, symmetric = FALSE, ...) { + + listTmp <- as.list(1:number_of_simulations) + tmp <- lapply(listTmp, function(l) + simulateOneTrajectory(trajectory_length, step_dist, + alpha, symmetric, ...)) + + result <- dplyr::bind_rows(lapply(listTmp, + function(x) tibble::tibble(sim_id = x, sim = tmp[[x]]))) + simulation <- list(simulation = result, + step_distribution = as.character(substitute(step_dist))[1], + alpha = alpha, + is_symmetric = symmetric) + class(simulation) <- c("kendall_simulation", "list") + simulation +} + + +#' Transforming (scaling and shifting) Kendall random walks +#' +#' If one trajectory has length n, an_seq and bn_seq arguments should be sequnces of length n. +#' Object returned by this function has plot and print methods. +#' +#' +#' @param simulations tibble returned by simulation function +#' @param an_seq sequence that the trajectories will be multiplied by +#' @param bn_seq sequence that will be substracted from scaled trajectory +#' +#' @return List like in simulate_kendall_rw function after transforming trajectories. +#' +#' @export +#' +#' @examples +#' kendall_simulations <- simulate_kendall_rw(10, 1000, runif, 0.5) +#' scaled_kendall <- transform_kendall_rw(kendall_simulations, (1:1000)^(-2)) +#' scaled_kendall # kendall random walked scaled by the sequence n^(-1/alpha) +#' scaled_data <- scaled_kendall$simulation # simulated values +#' plot(scaled_kendall) +#' + +transform_kendall_rw <- function(simulations, an_seq = 1, bn_seq = 0) { + sim_id <- sim <- NULL + result <- dplyr::mutate(simulations$simulation, + simNo = as.factor(as.character(sim_id))) + result <- dplyr::group_by(result, sim_id) + result <- dplyr::mutate(result, sim = an_seq*sim - bn_seq) + + simulations$simulation <- result + simulations +} + +#' Generic function that draws simulated trajectories of Kendall random walk +#' +#' @param x object returned by normalising_sequences function. +#' @param max_x maximum value on x axis. +#' @param max_id Number of trajectories to plot. If NULL, all paths will be plotted. +#' @param level Y-axis value which will be marked (level to be crossed). +#' @param ... Other arguments +#' +#' @return ggplot2 object +#' +#' @export +#' + +plot.kendall_simulation <- function(x, max_x = NULL, max_id = NULL, level = NULL, ...) { + sim_id <- NULL + n_sim <- max(unique(as.integer(as.character(x$simulation$sim_id)))) + trajectory_length <- dim(x$simulation)[1]/n_sim + if(is.null(max_x)) max_x <- trajectory_length + if(is.null(max_id)) max_id <- n_sim + + to_plot <- dplyr::ungroup(x$simulation) + to_plot <- dplyr::mutate(to_plot, x = rep(1:trajectory_length, n_sim)) + to_plot <- dplyr::filter(to_plot, x <= max_x, sim_id <= max_id) + plot_result <- ggplot2::ggplot(to_plot, + ggplot2::aes_string(x = 'x', y = 'sim', + group = 'sim_id')) + + ggplot2::geom_line() + + ggplot2::theme_bw() + + ggplot2::xlab("") + + ggplot2::ylab("") + + ggplot2::guides(color = "none") + if(!is.null(level)) + plot_result <- plot_result + ggplot2::geom_abline(slope = 0, + intercept = level, + color = "red", + size = 1.2) + plot_result +} + +#' Generic function that prints information about simulated Kendall random walk +#' +#' @param x Object returned by simulate_kendall_rw or transform_kendall_rw function. +#' @param ... Other arguments. +#' +#' @export +#' + +print.kendall_simulation <- function(x, ...) { + simulations_number <- max(as.numeric(as.character(x$simulation$sim_id))) + cat("Simulations of Kendall random walk \n") + cat("Number of simulations: ", simulations_number, "\n") + cat("Length of a single simulation: ", nrow(x$simulation)/simulations_number, "\n") + cat("Step distribution: ", x$step_distribution, "\n") + cat("Alpha parameter: ", x$alpha) + invisible(x) +} + +#' Calculate some characteristic for every simulated instance. +#' +#' @param simulations Object of class kendall_simulation. +#' @param summary_function Function that will be applied to each trajectory. +#' +#' @return data frame or a list (of class kendall_simulation) +#' +#' @export +#' + +summarise_kendall_rw <- function(simulations, summary_function) { + with(simulations$simulation, + dplyr::group_by(simulations$simulation, sim_id) %>% + dplyr::summarise(aggregated = summary_function(sim))) +} + + +#' Mutate each trajectory. +#' +#' @param simulations Object of class kendall_simulation. +#' @param mutate_function Function that will be applied to each trajectory. +#' @param df If TRUE, a d.f will be returned, if FALSE, simulations in the kendall_simulation +#' object passed in simulations argument will be replaced by the result of mutate_function. +#' +#' @return data frame or a list (of class kendall_simulation) +#' +#' @export +#' + +mutate_kendall_rw <- function(simulations, mutate_function, df = T) { + tmp <- with(simulations$simulation, + dplyr::group_by(simulations$simulation, sim_id) %>% + dplyr::mutate(sim = mutate_function(sim))) + if(df) { + tmp + } else { + simulations$simulation <- tmp + simulations + } + +} + diff --git a/R/stable_kendall_distribution.R b/R/stable_kendall_distribution.R new file mode 100644 index 0000000..79cb980 --- /dev/null +++ b/R/stable_kendall_distribution.R @@ -0,0 +1,172 @@ +#' PDF of Kendall stable distribution +#' +#' @param m_alpha function giving moments of order alpha of step dist. +#' +#' @return function that returns values of the PDF +#' +#' @export +#' +#' @examples +#' dKend <- dkend(function(x) 1) +#' # Step distribution: delta_{1} +#' dKendall <- dKend(1:10, 0.5) +#' # Values of PDF for arguments 1:10 and alpha = 0.5 +#' +#' + +dkend <- function(m_alpha) { + force(m_alpha) + function(x, alpha, mu = 0, sigma = 1) { + sapply(x, function(y) { + if(!is.finite(y)) return(NA) + else if((y - mu)/sigma <= 0) 0 + else (alpha/sigma)*(m_alpha(alpha)^2)*((y - mu)/sigma)^(-(2*alpha + 1))*exp(-m_alpha(alpha)*((y - mu)/sigma)^(-alpha)) + }) + } +} + + +#' CDF of Kendall stable distribution +#' +#' @param m_alpha function giving moments of order alpha of step dist. +#' +#' @return function function giving values of CDF of Kendall stable distribution +#' +#' @export +#' +#' @examples +#' pKend <- pkend(function(x) 1) +#' # Step distribution: delta_{1} +#' pKendall <- pKend(1:10, 0.5) +#' # Values of CDF for arguments 1:10 and alpha = 0.5 +#' +#' + +pkend <- function(m_alpha) { + force(m_alpha) + function(x, alpha, mu = 0, sigma = 1) { + sapply(x, function(y) { + if(!is.finite(y)) return(NA) + else { + if((y - mu)/sigma <= 0) 0 + else (1 + m_alpha(alpha)*((y - mu)/sigma)^(-alpha))*exp(-m_alpha(alpha)*((y - mu)/sigma)^(-alpha)) + } + }) + } +} + + +#' Quantiles of Kendall stable distribution +#' +#' @param m_alpha function giving moments of order alpha of step dist. +#' +#' @return function function returning quantiles of given orders +#' +#' @export +#' +#' @examples +#' qKend <- qkend(function(x) 1) +#' # Step distribution: delta_{1} +#' qKendall <- qKend(c(0.1, 0.9), 0.5) +#' # Quantiles of order 0.1 and 0.9 for alpha = 0.5 +#' +#' + +qkend<- function(m_alpha) { + force(m_alpha) + function(p, alpha, mu = 0, sigma = 1) { + oCDF <- function(x) pkend(m_alpha)(x, alpha, mu, sigma) + sapply(p, function(q) { + if(!is.finite(q)) return(NA) + else stats::uniroot({function(x) oCDF(x) - q}, lower = 0, upper = 10^80)$root + }) + } +} + + +#' Pseudo-random number from Kendall stable distribution +#' +#' @param m_alpha function giving moments of order alpha of step dist. +#' +#' @return function return n numbers genereted from Kendall stable dist. +#' +#' @importFrom stats rgamma +#' +#' @export +#' +#' @examples +#' rKend <- rkend(function(x) 1) +#' # Step distribution: delta_{1} +#' rKendall <- rKend(10, 0.5) +#' # Ten random number from stable Kendall distribution with alpha = 0.5 +#' +#' + + +rkend <- function(m_alpha) { + force(m_alpha) + function(n, alpha, mu = 0, sigma = 1) { + (sigma*rgamma(n, shape = 2, rate = m_alpha(alpha)) + mu)^(-1/alpha) + } +} + + +#' CDF of symmetrical Kendall stable distribution +#' +#' @param m_alpha function giving moments of order alpha of step dist. +#' +#' @return function function giving values of CDF of Kendall stable distribution +#' +#' @export +#' +#' @examples +#' pKend <- pkendSym(function(x) 1) +#' # Step distribution: delta_{1} +#' pKendall <- pKend(1:10, 0.5) +#' # Values of CDF for arguments 1:10 and alpha = 0.5 +#' +#' + +pkendSym <- function(m_alpha) { + force(m_alpha) + function(x, alpha, mu = 0, sigma = 1) { + Ft <- function(y) { + 0.5*(1 + m_alpha(alpha)*(y^((-1)*alpha)) + exp(m_alpha(alpha)*(y^((-1)*alpha))))*exp((-1)*m_alpha(alpha)*(y^((-1)*alpha))) + } + sapply(x, function(y) { + if(!is.finite(y)) return(NA) + else if(y == 0) 0.5 + else if(y > 0) Ft(y) + else 1 - Ft((-1)*y) + }) + } +} + + +#' Quantiles of symmetrical Kendall stable distribution +#' +#' @param m_alpha function giving moments of order alpha of step dist. +#' +#' @return function function returning quantiles of given orders +#' +#' @export +#' +#' @examples +#' qKend <- qkendSym(function(x) 1) +#' # Step distribution: delta_{1} +#' qKendall <- qKend(c(0.1, 0.9), 0.5) +#' # Quantiles of order 0.1 and 0.9 for alpha = 0.5 +#' +#' + +qkendSym <- function(m_alpha) { + force(m_alpha) + function(p, alpha, mu = 0, sigma = 1) { + oCDF <- function(x) pkendSym(m_alpha)(x, alpha) + sapply(p, function(q) { + if(!is.finite(q)) return(NA) + else if(q >= 0.5) stats::uniroot({function(x) oCDF(x) - q}, lower = 0, upper = 10^80)$root + else (-1)*stats::uniroot({function(x) oCDF(x) - (1 - q)}, lower = 0, upper = 10^80)$root + }) + } +} diff --git a/R/stable_kendall_fit.R b/R/stable_kendall_fit.R new file mode 100644 index 0000000..ab7a480 --- /dev/null +++ b/R/stable_kendall_fit.R @@ -0,0 +1,198 @@ +#' Log-likelihood for stable kendall distribution with m_alpha = 1 +#' +#' @param alpha alpha parameter of the Kendall random walk +#' @param x numeric vector of observations +#' +#' @return numeric +#' + +kendall_loglik <- function(alpha, x) { + x <- x[is.finite(x) & !is.na(x)] + length(x)*log(alpha) - (2*alpha + 1)*sum(log(x)) - sum(x^(-alpha)) +} + +#' Fit alpha parameter using MLE for distribution with one parameter +#' +#' @param data Numeric vector +#' +#' @importFrom stats optimize +#' +#' @return list with optimal parameter and loglikelihood value +#' +#' @export +#' + +estimate_stable_alpha <- function(data) { + to_optimize <- function(x) kendall_loglik(x, data) + optimize(to_optimize, c(0, 1), maximum = TRUE) +} + + +#' Fit stable Kendall distribution with one parameter (alpha) +#' +#' @param data Numeric vector of data. +#' +#' @return list of type kendall_fit +#' +#' @export +#' + +fit_stable_alpha <- function(data) { + alpha_param <- estimate_stable_alpha(data)$maximum + qKend <- qkend(function(x) 1) + theoretical_quantiles <- qKend((1:length(data) - 0.5)/length(data), alpha_param) + fitted_list <- list(fit = tibble::tibble(observed_quantiles = sort(data), + fitted_quantiles = sort(theoretical_quantiles)), + params = list(alpha = alpha_param)) + class(fitted_list) <- c("kendall_fit", "list") + fitted_list +} + +#' Negative loglikelihood for stable Kendall distr. with 3 parameters. +#' +#' @param data Dataset for which the loglikelihood will be calculated. +#' +#' @return numeric, value of loglikelihood +#' +#' @export +#' + +full_minus_loglik <- function(data) { + force(data) + n <- length(data) + function(parameters) { + # data2 <- data + alpha <- parameters[1] + location <- parameters[2] + scale <- parameters[3] + scaled <- (data - location)/scale + (2*alpha + 1)*sum(log(scaled)) + sum(scaled^(-alpha)) + n*log(scale) - n*log(alpha) + } +} + + +#' Gradient of minus loglikelihood for stable Kendall distribution with 3 parameters. +#' +#' @param data numeric vector of observation. +#' +#' @return Function of one argument of length 3 (alpha, location, scale). +#' +#' @export +#' + +full_loglik_gradient <- function(data) { + function(parameters) { + n <- length(data) + (-1)*c(n/parameters[1] + + sum(log((data - parameters[2])/parameters[3])* + (((data - parameters[2])/parameters[3])^(-parameters[1]) - 2)), + (2*parameters[1] + 1)*sum(1/(data - parameters[2])) - + (parameters[3]^(parameters[1]))* + parameters[1]*sum((data - parameters[2])^(-(parameters[1] + 1))), + (2*n*parameters[1])/parameters[3] - + parameters[1]*((parameters[3])^(parameters[1] - 1)) - + sum((data - parameters[2])^(-parameters[1])) + ) + } +} + + +#' Fit stable Kendall distribution for given data and m_alpha function. +#' +#' @param data Numeric vector of observation to which the distribution will be fitted. +#' +#' @importFrom stats optim +#' +#' @return fitted quantiles +#' +#' @export +#' + +fit_kendall <- function(data) { + alpha_start <- estimate_stable_alpha(data)$maximum + loglik_data <- full_minus_loglik(data) + loglik_grad_data <- full_loglik_gradient(data) + to_estimate <- optim(c(alpha_start, 0.1, 0.1), loglik_data, loglik_grad_data, + lower = c(0.01, -100, 0.001), + upper = c(0.99, min(data)-0.1, 100), + method = "L-BFGS-B") + alpha <- to_estimate$par[1] + loc <- to_estimate$par[2] + scale <- to_estimate$par[3] + quantiles_function <- qkend(function(x) 1) + fitted_quantiles <- quantiles_function((1:length(data) - 0.5)/length(data), + alpha, loc, scale) + result <- list(fit = tibble::tibble(fitted_quantiles = fitted_quantiles, + observed_quantiles = sort(data)), + params = list(estimated_alpha = alpha, + estimated_location = loc, + estimated_scale = scale)) + class(result) <- c("kendall_fit", "list") + result +} + + +#' Function for fitting stable Kendall distribution separately to two parts of data +#' +#' @param data Numeric vector. Observation to which the distribution will be fitted. +#' @param separation_point Order above which data (quantiles) will be separated. +#' +#' @return List of class kendall_fit with estimated and theoretical quantiles +#' and estimated parameters. +#' +#' @export +#' + +fit_separate <- function(data, separation_point) { + data <- data[is.finite(data) & !is.na(data)] + n_obs <- length(data) + if(n_obs == 0) { + stop("No legal observations were provided") + } + kendall_quantiles <- qkend(function(x) 1) + separate <- stats::quantile(data, separation_point) + all_quantiles <- (1:n_obs - 0.5)/n_obs + lower_quantiles <- all_quantiles[all_quantiles <= separation_point] + upper_quantiles <- all_quantiles[all_quantiles > separation_point] + lower_dataset <- data[data <= separate] + upper_dataset <- data[data > separate] + all_lower_fit <- fit_kendall(lower_dataset) + all_upper_fit <- fit_kendall(upper_dataset) + lower_fitted<- unlist(all_lower_fit[[1]]$fitted_quantiles, + use.names = F) + upper_fitted <- unlist(all_upper_fit[[1]]$fitted_quantiles, + use.names = F) + quantiles <- unname(c(lower_fitted, upper_fitted)) + result <- list(fit = tibble::tibble(observed_quantiles = sort(data), + fitted_quantiles = quantiles), + params = list(alpha_lower = all_lower_fit[[2]]$estimated_alpha, + location_lower = all_lower_fit[[2]]$estimated_location, + scale_lower = all_lower_fit[[2]]$estimated_scale, + alpha_upper = all_upper_fit[[2]]$estimated_alpha, + location_upper = all_upper_fit[[2]]$estimated_location, + scale_upper = all_upper_fit[[2]]$estimated_scale)) + class(result) <- c("kendall_fit", "list") + result +} + + +#' QQ-plot for the result of fitting stable Kendall distribtion. +#' +#' @param x List returned by fit_separate or fit_kendall function. +#' @param ... Aditional arguments. +#' +#' @return ggplot2 object +#' +#' @export +#' + +plot.kendall_fit <- function(x, ...) { + observed_vs_fitted <- x$fit + ggplot2::ggplot(observed_vs_fitted, ggplot2::aes_string(x = 'observed_quantiles', + y = 'fitted_quantiles')) + + ggplot2::theme_bw() + + ggplot2::xlab("Empirical") + + ggplot2::ylab("Theoretical quantile") + + ggplot2::geom_point() +} + diff --git a/README.md b/README.md new file mode 100644 index 0000000..b3fc6a2 --- /dev/null +++ b/README.md @@ -0,0 +1,35 @@ +# kendallRandomPackage + +[![Travis-CI Build Status](https://travis-ci.org/mstaniak/kendallRandomPackage.svg?branch=master)](https://travis-ci.org/mstaniak/kendallRandomPackage) +[![Coverage Status](https://img.shields.io/codecov/c/github/mstaniak/kendallRandomPackage/master.svg)](https://codecov.io/github/mstaniak/kendallRandomPackage?branch=master) + + + +Simulations and distributions related to Kendall random walks: +[visit dedicated project page on Researchgate](https://www.researchgate.net/project/First-order-Kendall-maximal-autoregressive-processes-and-their-applications) + +Install the newest version: + +>devtools::install_github("mstaniak/kendallRandomPackage") + +Help: +>?kendallRandomPackage + +Main functionalities: + + * `simulate_kendall_rw` functions simulates Kendall random walks for a given step distribution. The simulation can be then plotted using generic `plot` function. + * `transform_kendall_rw` function allows user to play with different scalings and transformations of Kendall random walks to study its properties related to convergence. + * `ladder_moment` and `ladder_height` functions help study distribution of first ladder moment and first ladder height empirically. + * `ladder_moment_pmf` gives exact PMF for first ladder moment. + * `pkend`, `dkend`, `qkend` and `rkend` are typical functions related to the stable Kendall distribution. + * `g_function` calculates Williamson transform numerically. + +See vignette for examples. + +If you have a feature request or you found a bug, please leave an issue. +The goal of this package is to let anyone interested in Kendall convolution/generalized convolutions get familiar with them through visual means and experimentation. + + +Acknowledgement + +This work is a part of project "First order Kendall maximal autoregressive processes and their applications", which is carried out within the POWROTY/REINTEGRATION programme of the Foundation for Polish Science co-financed by the European Union under the European Regional Development Fund. diff --git a/build/vignette.rds b/build/vignette.rds new file mode 100644 index 0000000000000000000000000000000000000000..698bf33b3499801eb0b9780268bc70817fa1c1ff GIT binary patch literal 248 zcmVKWu!2kb_*$L5>SYcM72cYf^I$^jO(#t|HMdoR z#tddGYXSh%`u$6X8v7saJy}Z^Yv&|m#%f1HZ=P^w9;_4fPTWN}*T% + mutate(id = 1:n()) +lengths <- diffs %>% + filter(sim != 0) %>% + mutate(previous = ifelse(is.na(lag(id)), 0, lag(id))) %>% + mutate(length = id - previous) +ggplot(subset(lengths, sim_id < 5), + aes(x = length, fill = as.factor(sim_id), group = as.factor(sim_id))) + + geom_density() + + theme_bw() + + ggtitle("Distribution of time with no state-change (by simulation)") +ggplot(lengths, aes(x = length)) + + geom_density() + + theme_bw() + + ggtitle("Distribution of time with no state-change (aggregated)") +ggplot(subset(lengths, sim_id < 10), aes(x = id, y = length, color = as.factor(sim_id))) + + geom_point() + + theme_bw() + + geom_line() + + ggtitle("Time with no state-change in time") + diff --git a/inst/doc/behaviour.Rmd b/inst/doc/behaviour.Rmd new file mode 100644 index 0000000..3fc1692 --- /dev/null +++ b/inst/doc/behaviour.Rmd @@ -0,0 +1,81 @@ +--- +title: "Studying the behaviour of Kendall random walks" +author: "Mateusz Staniak" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Studying the behaviour of Kendall random walks} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + warning = F, + message = F, + collapse = TRUE, + comment = "#>" +) +``` + +We will approximate the distribution of moments when the random walk changes state through simulations. + +First, we simulate many paths of Kendall random walk with normal step distribution. + +```{r sim} +library(kendallRandomWalks) +library(dplyr) +library(ggplot2) +set.seed(17) +walks <- simulate_kendall_rw(1000, 1000, rnorm, 0.5, T) +``` + +Example trajectory + +```{r ex1} +plot(walks, max_id = 1) +``` + +Number of unique states + +```{r unique} +ggplot(summarise_kendall_rw(walks, n_distinct), aes(x = aggregated)) + + # geom_histogram() + + theme_bw() + + geom_density() +``` + +Jumps + +```{r jumps} +diffs <- mutate_kendall_rw(walks, function(x) x - lag(x)) +diffs +diffs2 <- mutate_kendall_rw(walks, function(x) x - lag(x), F) +plot(diffs2, max_id = 10) +``` + +Time with no change of state + +```{r skoki_dlugosc} +diffs3 <- mutate_kendall_rw(diffs2, function(x) as.numeric(x != 0), F) +diffs <- diffs3$simulation +diffs <- group_by(diffs, sim_id) %>% + mutate(id = 1:n()) +lengths <- diffs %>% + filter(sim != 0) %>% + mutate(previous = ifelse(is.na(lag(id)), 0, lag(id))) %>% + mutate(length = id - previous) +ggplot(subset(lengths, sim_id < 5), + aes(x = length, fill = as.factor(sim_id), group = as.factor(sim_id))) + + geom_density() + + theme_bw() + + ggtitle("Distribution of time with no state-change (by simulation)") +ggplot(lengths, aes(x = length)) + + geom_density() + + theme_bw() + + ggtitle("Distribution of time with no state-change (aggregated)") +ggplot(subset(lengths, sim_id < 10), aes(x = id, y = length, color = as.factor(sim_id))) + + geom_point() + + theme_bw() + + geom_line() + + ggtitle("Time with no state-change in time") +``` diff --git a/inst/doc/behaviour.html b/inst/doc/behaviour.html new file mode 100644 index 0000000..5712caa --- /dev/null +++ b/inst/doc/behaviour.html @@ -0,0 +1,152 @@ + + + + + + + + + + + + + + + +Studying the behaviour of Kendall random walks + + + + + + + + + + + + + + + + + +

Studying the behaviour of Kendall random walks

+

Mateusz Staniak

+ + + +

We will approximate the distribution of moments when the random walk changes state through simulations.

+

First, we simulate many paths of Kendall random walk with normal step distribution.

+
library(kendallRandomWalks)
+library(dplyr)
+library(ggplot2)
+set.seed(17)
+walks <- simulate_kendall_rw(1000, 1000, rnorm, 0.5, T)
+

Example trajectory

+
plot(walks, max_id = 1)
+

+

Number of unique states

+
ggplot(summarise_kendall_rw(walks, n_distinct), aes(x = aggregated)) +
+  # geom_histogram() +
+  theme_bw() +
+  geom_density()
+

+

Jumps

+
diffs <- mutate_kendall_rw(walks, function(x) x - lag(x))
+diffs
+#> # A tibble: 1,000,000 x 2
+#> # Groups:   sim_id [1,000]
+#>    sim_id     sim
+#>     <int>   <dbl>
+#>  1      1  NA    
+#>  2      1  -1.02 
+#>  3      1   0.   
+#>  4      1  -0.478
+#>  5      1   3.88 
+#>  6      1  -5.04 
+#>  7      1  26.0  
+#>  8      1   0.   
+#>  9      1   0.   
+#> 10      1   0.   
+#> # ... with 999,990 more rows
+diffs2 <- mutate_kendall_rw(walks, function(x) x - lag(x), F)
+plot(diffs2, max_id = 10)
+

+

Time with no change of state

+
diffs3 <- mutate_kendall_rw(diffs2, function(x) as.numeric(x != 0), F)
+diffs <- diffs3$simulation
+diffs <- group_by(diffs, sim_id) %>%
+  mutate(id = 1:n())
+lengths <- diffs %>%
+  filter(sim != 0) %>%
+  mutate(previous = ifelse(is.na(lag(id)), 0, lag(id))) %>%
+  mutate(length = id - previous)
+ggplot(subset(lengths, sim_id < 5), 
+       aes(x = length, fill = as.factor(sim_id), group = as.factor(sim_id))) +
+  geom_density() +
+  theme_bw() +
+  ggtitle("Distribution of time with no state-change (by simulation)")
+

+
ggplot(lengths, aes(x = length)) +
+  geom_density() +
+  theme_bw() +
+  ggtitle("Distribution of time with no state-change (aggregated)")
+

+
ggplot(subset(lengths, sim_id < 10), aes(x = id, y = length, color = as.factor(sim_id))) +
+  geom_point() +
+  theme_bw() +
+  geom_line() +
+  ggtitle("Time with no state-change in time")
+

+ + + + + + + + diff --git a/inst/doc/kendall_rws.R b/inst/doc/kendall_rws.R new file mode 100644 index 0000000..4c2ecc6 --- /dev/null +++ b/inst/doc/kendall_rws.R @@ -0,0 +1,41 @@ +## ----setup, include = FALSE---------------------------------------------- +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) + +## ----example1------------------------------------------------------------ +# uzupelnic fig_captions, zeby dalo sie je zrobic +library(kendallRandomWalks) +kendall_rws <- simulate_kendall_rw(10, 100, runif, 0.25) +kendall_rws +plot(kendall_rws) +plot(simulate_kendall_rw(10, 100, rnorm, 0.76), level = 300) + +## ----example1.1---------------------------------------------------------- +kendall_rws_sym <- simulate_kendall_rw(10, 100, rnorm, 0.76, T) +kendall_rws_sym +plot(kendall_rws_sym) + +## ----example2------------------------------------------------------------ +kendall_rws2 <- simulate_kendall_rw(1000, 100, runif, 0.25) +ladder_moments <- ladder_moment(kendall_rws2, 1000) +ladder_moments +plot(ladder_moments) + +ladder_heights <- ladder_height(kendall_rws2, 2000) +ladder_heights +plot(ladder_heights) + +## ----example3------------------------------------------------------------ +y <- seq(10, 10000, 50) +ladders <- sapply(y, + function(x) + ladder_moment_pmf(10, x, 0.5, pnorm, dnorm)) +plot(y, ladders) + +y <- seq(2000, 2200, 1) +plot(y, g_function(y, 0.1, dnorm)) + +plot(seq(0, 400, by = 1), g_function(seq(0, 400, by = 1), 0.5, dunif)) + diff --git a/inst/doc/kendall_rws.Rmd b/inst/doc/kendall_rws.Rmd new file mode 100644 index 0000000..b7cc72b --- /dev/null +++ b/inst/doc/kendall_rws.Rmd @@ -0,0 +1,67 @@ +--- +title: "Kendall random walks" +author: "Mateusz Staniak" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Kendall random walks} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + + +## Simulate and plot + +```{r example1} +# uzupelnic fig_captions, zeby dalo sie je zrobic +library(kendallRandomWalks) +kendall_rws <- simulate_kendall_rw(10, 100, runif, 0.25) +kendall_rws +plot(kendall_rws) +plot(simulate_kendall_rw(10, 100, rnorm, 0.76), level = 300) +``` + +Symmetric + +```{r example1.1} +kendall_rws_sym <- simulate_kendall_rw(10, 100, rnorm, 0.76, T) +kendall_rws_sym +plot(kendall_rws_sym) +``` + +## Barrier crossing + +```{r example2} +kendall_rws2 <- simulate_kendall_rw(1000, 100, runif, 0.25) +ladder_moments <- ladder_moment(kendall_rws2, 1000) +ladder_moments +plot(ladder_moments) + +ladder_heights <- ladder_height(kendall_rws2, 2000) +ladder_heights +plot(ladder_heights) +``` + +Exact first ladder moments distribution with $G(a)$ computed numerically. + +```{r example3} +y <- seq(10, 10000, 50) +ladders <- sapply(y, + function(x) + ladder_moment_pmf(10, x, 0.5, pnorm, dnorm)) +plot(y, ladders) + +y <- seq(2000, 2200, 1) +plot(y, g_function(y, 0.1, dnorm)) + +plot(seq(0, 400, by = 1), g_function(seq(0, 400, by = 1), 0.5, dunif)) +``` + + diff --git a/inst/doc/kendall_rws.html b/inst/doc/kendall_rws.html new file mode 100644 index 0000000..87538f8 --- /dev/null +++ b/inst/doc/kendall_rws.html @@ -0,0 +1,164 @@ + + + + + + + + + + + + + + + + +Kendall random walks + + + + + + + + + + + + + + + + + +

Kendall random walks

+

Mateusz Staniak

+

2018-06-05

+ + + +
+

Simulate and plot

+
# uzupelnic fig_captions, zeby dalo sie je zrobic
+library(kendallRandomWalks)
+kendall_rws <- simulate_kendall_rw(10, 100, runif, 0.25)
+kendall_rws
+#> Simulations of Kendall random walk 
+#> Number of simulations:  10 
+#> Length of a single simulation:  100 
+#> Step distribution:  runif 
+#> Alpha parameter:  0.25
+plot(kendall_rws)
+

+
plot(simulate_kendall_rw(10, 100, rnorm, 0.76), level = 300)
+

+

Symmetric

+
kendall_rws_sym <- simulate_kendall_rw(10, 100, rnorm, 0.76, T)
+kendall_rws_sym
+#> Simulations of Kendall random walk 
+#> Number of simulations:  10 
+#> Length of a single simulation:  100 
+#> Step distribution:  rnorm 
+#> Alpha parameter:  0.76
+plot(kendall_rws_sym)
+

+
+
+

Barrier crossing

+
kendall_rws2 <- simulate_kendall_rw(1000, 100, runif, 0.25)
+ladder_moments <- ladder_moment(kendall_rws2, 1000)
+ladder_moments
+#> Mean of the distribution:  15.726 
+#> Standard deviation of the distribution:  9.839397 
+#> Number of observations:  1000 
+#> Times the level was not crossed:  0 
+#> Quantiles of the distribution: 
+#>   0%  10%  20%  30%  40%  50%  60%  70%  80%  90% 100% 
+#>    3    6    8   10   11   13   15   19   23   29   64
+plot(ladder_moments)
+#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
+

+

+ladder_heights <- ladder_height(kendall_rws2, 2000)
+ladder_heights
+#> Mean of the distribution:  8438489 
+#> Standard deviation of the distribution:  181627127 
+#> Number of observations:  1000 
+#> Times the level was not crossed:  0 
+#> Quantiles of the distribution: 
+#>           0%          10%          20%          30%          40% 
+#> 2.001583e+03 2.374325e+03 3.013069e+03 3.928392e+03 5.260592e+03 
+#>          50%          60%          70%          80%          90% 
+#> 7.180237e+03 1.095814e+04 1.862697e+04 4.335170e+04 1.705180e+05 
+#>         100% 
+#> 5.583728e+09
+plot(ladder_heights)
+#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
+

+

Exact first ladder moments distribution with \(G(a)\) computed numerically.

+
y <- seq(10, 10000, 50)
+ladders <- sapply(y, 
+                  function(x) 
+                  ladder_moment_pmf(10, x, 0.5, pnorm, dnorm))
+plot(y, ladders)
+

+

+y <- seq(2000, 2200, 1)
+plot(y, g_function(y, 0.1, dnorm))
+

+

+plot(seq(0, 400, by = 1), g_function(seq(0, 400, by = 1), 0.5, dunif))
+

+
+ + + + + + + + diff --git a/man/Qn.Rd b/man/Qn.Rd new file mode 100644 index 0000000..d316914 --- /dev/null +++ b/man/Qn.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/simulations_kendall_rw.R +\name{Qn} +\alias{Qn} +\title{Helper function} +\usage{ +Qn(x, y, alpha) +} +\arguments{ +\item{x}{numeric} + +\item{y}{numeric} + +\item{alpha}{numeric, parameter of Kendall random walk} +} +\value{ +0 or 1 with probability depending on x, y, alpha +} +\description{ +Helper function +} diff --git a/man/U.Rd b/man/U.Rd new file mode 100644 index 0000000..a44e79d --- /dev/null +++ b/man/U.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/simulations_kendall_rw.R +\name{U} +\alias{U} +\title{Helper function} +\usage{ +U(x, y) +} +\arguments{ +\item{x}{numeric} + +\item{y}{numeric} +} +\value{ +sign of the argument whose abs. val. is bigger +} +\description{ +Helper function +} diff --git a/man/Z.Rd b/man/Z.Rd new file mode 100644 index 0000000..fe570f5 --- /dev/null +++ b/man/Z.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/simulations_kendall_rw.R +\name{Z} +\alias{Z} +\title{Helper function: min/max} +\usage{ +Z(x, y) +} +\arguments{ +\item{x}{numeric} + +\item{y}{numeric} +} +\value{ +min of arguments divided by max of arguments +} +\description{ +Helper function: min/max +} diff --git a/man/dkend.Rd b/man/dkend.Rd new file mode 100644 index 0000000..a8b2d07 --- /dev/null +++ b/man/dkend.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stable_kendall_distribution.R +\name{dkend} +\alias{dkend} +\title{PDF of Kendall stable distribution} +\usage{ +dkend(m_alpha) +} +\arguments{ +\item{m_alpha}{function giving moments of order alpha of step dist.} +} +\value{ +function that returns values of the PDF +} +\description{ +PDF of Kendall stable distribution +} +\examples{ +dKend <- dkend(function(x) 1) +# Step distribution: delta_{1} +dKendall <- dKend(1:10, 0.5) +# Values of PDF for arguments 1:10 and alpha = 0.5 + + +} diff --git a/man/estimate_stable_alpha.Rd b/man/estimate_stable_alpha.Rd new file mode 100644 index 0000000..7d768ac --- /dev/null +++ b/man/estimate_stable_alpha.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stable_kendall_fit.R +\name{estimate_stable_alpha} +\alias{estimate_stable_alpha} +\title{Fit alpha parameter using MLE for distribution with one parameter} +\usage{ +estimate_stable_alpha(data) +} +\arguments{ +\item{data}{Numeric vector} +} +\value{ +list with optimal parameter and loglikelihood value +} +\description{ +Fit alpha parameter using MLE for distribution with one parameter +} diff --git a/man/fit_kendall.Rd b/man/fit_kendall.Rd new file mode 100644 index 0000000..c74a0da --- /dev/null +++ b/man/fit_kendall.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stable_kendall_fit.R +\name{fit_kendall} +\alias{fit_kendall} +\title{Fit stable Kendall distribution for given data and m_alpha function.} +\usage{ +fit_kendall(data) +} +\arguments{ +\item{data}{Numeric vector of observation to which the distribution will be fitted.} +} +\value{ +fitted quantiles +} +\description{ +Fit stable Kendall distribution for given data and m_alpha function. +} diff --git a/man/fit_separate.Rd b/man/fit_separate.Rd new file mode 100644 index 0000000..8124eed --- /dev/null +++ b/man/fit_separate.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stable_kendall_fit.R +\name{fit_separate} +\alias{fit_separate} +\title{Function for fitting stable Kendall distribution separately to two parts of data} +\usage{ +fit_separate(data, separation_point) +} +\arguments{ +\item{data}{Numeric vector. Observation to which the distribution will be fitted.} + +\item{separation_point}{Order above which data (quantiles) will be separated.} +} +\value{ +List of class kendall_fit with estimated and theoretical quantiles +and estimated parameters. +} +\description{ +Function for fitting stable Kendall distribution separately to two parts of data +} diff --git a/man/fit_stable_alpha.Rd b/man/fit_stable_alpha.Rd new file mode 100644 index 0000000..d96e9c3 --- /dev/null +++ b/man/fit_stable_alpha.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stable_kendall_fit.R +\name{fit_stable_alpha} +\alias{fit_stable_alpha} +\title{Fit stable Kendall distribution with one parameter (alpha)} +\usage{ +fit_stable_alpha(data) +} +\arguments{ +\item{data}{Numeric vector of data.} +} +\value{ +list of type kendall_fit +} +\description{ +Fit stable Kendall distribution with one parameter (alpha) +} diff --git a/man/full_loglik_gradient.Rd b/man/full_loglik_gradient.Rd new file mode 100644 index 0000000..ad62478 --- /dev/null +++ b/man/full_loglik_gradient.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stable_kendall_fit.R +\name{full_loglik_gradient} +\alias{full_loglik_gradient} +\title{Gradient of minus loglikelihood for stable Kendall distribution with 3 parameters.} +\usage{ +full_loglik_gradient(data) +} +\arguments{ +\item{data}{numeric vector of observation.} +} +\value{ +Function of one argument of length 3 (alpha, location, scale). +} +\description{ +Gradient of minus loglikelihood for stable Kendall distribution with 3 parameters. +} diff --git a/man/full_minus_loglik.Rd b/man/full_minus_loglik.Rd new file mode 100644 index 0000000..e408bec --- /dev/null +++ b/man/full_minus_loglik.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stable_kendall_fit.R +\name{full_minus_loglik} +\alias{full_minus_loglik} +\title{Negative loglikelihood for stable Kendall distr. with 3 parameters.} +\usage{ +full_minus_loglik(data) +} +\arguments{ +\item{data}{Dataset for which the loglikelihood will be calculated.} +} +\value{ +numeric, value of loglikelihood +} +\description{ +Negative loglikelihood for stable Kendall distr. with 3 parameters. +} diff --git a/man/g_function.Rd b/man/g_function.Rd new file mode 100644 index 0000000..bea7551 --- /dev/null +++ b/man/g_function.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/exact_ladder_moments.R +\name{g_function} +\alias{g_function} +\title{Function G(t) - Williamson transform taken at point 1/t.} +\usage{ +g_function(t, alpha, density) +} +\arguments{ +\item{t}{Argument to the function.} + +\item{alpha}{Value of the alpha parameter.} + +\item{density}{Density function of the step distribution.} +} +\value{ +Object of class "integrate" +} +\description{ +This function return only approximated values. To check their precisions use +g_function_single function with an argument of length 1. +} +\examples{ +g_function(1:5, 0.75, dnorm) + + +} diff --git a/man/g_function_single.Rd b/man/g_function_single.Rd new file mode 100644 index 0000000..ee096a0 --- /dev/null +++ b/man/g_function_single.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/exact_ladder_moments.R +\name{g_function_single} +\alias{g_function_single} +\title{Function G(t) - Williamson transform taken at point 1/t.} +\usage{ +g_function_single(t, alpha, density) +} +\arguments{ +\item{t}{Argument to the function.} + +\item{alpha}{Value of the alpha parameter.} + +\item{density}{Density function of the step distribution.} +} +\value{ +Object of class "integrate" +} +\description{ +This function return the whole "integrate" object, so precision of the approximation +can be checked. +} +\examples{ +g_function_single(5, 0.26, dnorm) + +} diff --git a/man/kendallRandomWalks.Rd b/man/kendallRandomWalks.Rd new file mode 100644 index 0000000..d473fab --- /dev/null +++ b/man/kendallRandomWalks.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/kendallRandomPackage.R +\docType{package} +\name{kendallRandomWalks} +\alias{kendallRandomWalks} +\alias{kendallRandomWalks-package} +\title{kendallRandomWalks: explore and visualize Kendall random walks.} +\description{ +Kendall random walks are Markov processes generated by the Kendall convolution. +This package helps simulate and visualize these random walks and associated +distributions. It also provides function to fit these distributions to data. +} +\section{Important functions}{ + +\code{\link{simulate_kendall_rw}} simulates Kendall random walks. + +\code{\link{transform_kendall_rw}} applies scaling and shift to simulated Kendall r.w-s. + +\code{\link{ladder_moment}} estimates the distribution of first ladder moment +by simulating Kendall random walks. + +\code{\link{ladder_height}} estimates the distribution of first ladder height +by simulating Kendall random walks. + +\code{\link{ladder_moment_pmf}} computes the PMF of the distribution of first +ladder moment. + +\code{\link{g_function}} Finds the value of G(t) numerically. + +\code{\link{pkend}}, \code{\link{dkend}}, \code{\link{qkend}}, \code{\link{rkend}} +give CDF, PDF, quantile function and random numbers from stable Kendall distribution. +} + diff --git a/man/kendall_loglik.Rd b/man/kendall_loglik.Rd new file mode 100644 index 0000000..b0efb37 --- /dev/null +++ b/man/kendall_loglik.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stable_kendall_fit.R +\name{kendall_loglik} +\alias{kendall_loglik} +\title{Log-likelihood for stable kendall distribution with m_alpha = 1} +\usage{ +kendall_loglik(alpha, x) +} +\arguments{ +\item{alpha}{alpha parameter of the Kendall random walk} + +\item{x}{numeric vector of observations} +} +\value{ +numeric +} +\description{ +Log-likelihood for stable kendall distribution with m_alpha = 1 +} diff --git a/man/ladder_height.Rd b/man/ladder_height.Rd new file mode 100644 index 0000000..8f56a78 --- /dev/null +++ b/man/ladder_height.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/barrier_crossing.R +\name{ladder_height} +\alias{ladder_height} +\title{Estimate the distribution of first ladder height for given level} +\usage{ +ladder_height(simulations, level) +} +\arguments{ +\item{simulations}{kendall_simulation object} + +\item{level}{Positive numeric} +} +\value{ +tibble +} +\description{ +NA is returned if the level wasn't crossed. Printing the resulting object +will give summary of the estimated distribution and information whether +level wasn't crossed in some simulations. This information can be used to +pick the right trajectory length for the given level. +} +\examples{ +{ + kendall_rw <- simulate_kendall_rw(100, 100, runif, 0.5) + estim_ladder <- ladder_height(kendall_rw, 1000) + estim_ladder +} + +} diff --git a/man/ladder_moment.Rd b/man/ladder_moment.Rd new file mode 100644 index 0000000..16cc58a --- /dev/null +++ b/man/ladder_moment.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/barrier_crossing.R +\name{ladder_moment} +\alias{ladder_moment} +\title{Estimate the distribution of first ladder moment for given level} +\usage{ +ladder_moment(simulations, level) +} +\arguments{ +\item{simulations}{kendall_simulation object} + +\item{level}{Positive numeric} +} +\value{ +tibble +} +\description{ +NA is returned if the level wasn't crossed. Printing the resulting object +will give summary of the estimated distribution and information whether +level wasn't crossed in some simulations. This information can be used to +pick the right trajectory length for the given level. +} +\examples{ +{ +kendall_rw <- simulate_kendall_rw(100, 100, runif, 0.5) +estim_ladder <- ladder_moment(kendall_rw, 1000) +estim_ladder +} + +} diff --git a/man/ladder_moment_pmf.Rd b/man/ladder_moment_pmf.Rd new file mode 100644 index 0000000..882a270 --- /dev/null +++ b/man/ladder_moment_pmf.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/exact_ladder_moments.R +\name{ladder_moment_pmf} +\alias{ladder_moment_pmf} +\title{Distribution of the first ladder moment.} +\usage{ +ladder_moment_pmf(n, level, alpha, step_cdf, step_pdf) +} +\arguments{ +\item{n}{Argument to the PDF.} + +\item{level}{Level a to be crossed.} + +\item{alpha}{Alpha parameter of Kendall random walk.} + +\item{step_cdf}{CDF of the step distribution.} + +\item{step_pdf}{PDF of the step distribution.} +} +\value{ +Value of PMF of the distribution of first ladder moment +} +\description{ +Distribution of the first ladder moment. +} +\examples{ +prob <- ladder_moment_pmf(10, 1000, 0.5, pnorm, dnorm) +prob + + +} diff --git a/man/mutate_kendall_rw.Rd b/man/mutate_kendall_rw.Rd new file mode 100644 index 0000000..caacf2f --- /dev/null +++ b/man/mutate_kendall_rw.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/simulations_kendall_rw.R +\name{mutate_kendall_rw} +\alias{mutate_kendall_rw} +\title{Mutate each trajectory.} +\usage{ +mutate_kendall_rw(simulations, mutate_function, df = T) +} +\arguments{ +\item{simulations}{Object of class kendall_simulation.} + +\item{mutate_function}{Function that will be applied to each trajectory.} + +\item{df}{If TRUE, a d.f will be returned, if FALSE, simulations in the kendall_simulation +object passed in simulations argument will be replaced by the result of mutate_function.} +} +\value{ +data frame or a list (of class kendall_simulation) +} +\description{ +Mutate each trajectory. +} diff --git a/man/pkend.Rd b/man/pkend.Rd new file mode 100644 index 0000000..a614965 --- /dev/null +++ b/man/pkend.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stable_kendall_distribution.R +\name{pkend} +\alias{pkend} +\title{CDF of Kendall stable distribution} +\usage{ +pkend(m_alpha) +} +\arguments{ +\item{m_alpha}{function giving moments of order alpha of step dist.} +} +\value{ +function function giving values of CDF of Kendall stable distribution +} +\description{ +CDF of Kendall stable distribution +} +\examples{ +pKend <- pkend(function(x) 1) +# Step distribution: delta_{1} +pKendall <- pKend(1:10, 0.5) +# Values of CDF for arguments 1:10 and alpha = 0.5 + + +} diff --git a/man/pkendSym.Rd b/man/pkendSym.Rd new file mode 100644 index 0000000..cdda53b --- /dev/null +++ b/man/pkendSym.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stable_kendall_distribution.R +\name{pkendSym} +\alias{pkendSym} +\title{CDF of symmetrical Kendall stable distribution} +\usage{ +pkendSym(m_alpha) +} +\arguments{ +\item{m_alpha}{function giving moments of order alpha of step dist.} +} +\value{ +function function giving values of CDF of Kendall stable distribution +} +\description{ +CDF of symmetrical Kendall stable distribution +} +\examples{ +pKend <- pkendSym(function(x) 1) +# Step distribution: delta_{1} +pKendall <- pKend(1:10, 0.5) +# Values of CDF for arguments 1:10 and alpha = 0.5 + + +} diff --git a/man/plot.kendall_barrier_crossing.Rd b/man/plot.kendall_barrier_crossing.Rd new file mode 100644 index 0000000..fc075d1 --- /dev/null +++ b/man/plot.kendall_barrier_crossing.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/barrier_crossing.R +\name{plot.kendall_barrier_crossing} +\alias{plot.kendall_barrier_crossing} +\title{Generic function for plotting results of ladder_moment function.} +\usage{ +\method{plot}{kendall_barrier_crossing}(x, ...) +} +\arguments{ +\item{x}{kendall_barrier_crossing object} + +\item{...}{Additional arguments} +} +\value{ +ggplot2 +} +\description{ +Generic function for plotting results of ladder_moment function. +} diff --git a/man/plot.kendall_fit.Rd b/man/plot.kendall_fit.Rd new file mode 100644 index 0000000..6df090d --- /dev/null +++ b/man/plot.kendall_fit.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stable_kendall_fit.R +\name{plot.kendall_fit} +\alias{plot.kendall_fit} +\title{QQ-plot for the result of fitting stable Kendall distribtion.} +\usage{ +\method{plot}{kendall_fit}(x, ...) +} +\arguments{ +\item{x}{List returned by fit_separate or fit_kendall function.} + +\item{...}{Aditional arguments.} +} +\value{ +ggplot2 object +} +\description{ +QQ-plot for the result of fitting stable Kendall distribtion. +} diff --git a/man/plot.kendall_simulation.Rd b/man/plot.kendall_simulation.Rd new file mode 100644 index 0000000..c92ed17 --- /dev/null +++ b/man/plot.kendall_simulation.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/simulations_kendall_rw.R +\name{plot.kendall_simulation} +\alias{plot.kendall_simulation} +\title{Generic function that draws simulated trajectories of Kendall random walk} +\usage{ +\method{plot}{kendall_simulation}(x, max_x = NULL, max_id = NULL, + level = NULL, ...) +} +\arguments{ +\item{x}{object returned by normalising_sequences function.} + +\item{max_x}{maximum value on x axis.} + +\item{max_id}{Number of trajectories to plot. If NULL, all paths will be plotted.} + +\item{level}{Y-axis value which will be marked (level to be crossed).} + +\item{...}{Other arguments} +} +\value{ +ggplot2 object +} +\description{ +Generic function that draws simulated trajectories of Kendall random walk +} diff --git a/man/print.kendall_barrier_crossing.Rd b/man/print.kendall_barrier_crossing.Rd new file mode 100644 index 0000000..18d32a6 --- /dev/null +++ b/man/print.kendall_barrier_crossing.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/barrier_crossing.R +\name{print.kendall_barrier_crossing} +\alias{print.kendall_barrier_crossing} +\title{Generic function for printing result of ladder_moment function} +\usage{ +\method{print}{kendall_barrier_crossing}(x, ...) +} +\arguments{ +\item{x}{kendall_barrier_crossing object} + +\item{...}{Additional arguments} +} +\value{ +invisible x +} +\description{ +Generic function for printing result of ladder_moment function +} diff --git a/man/print.kendall_simulation.Rd b/man/print.kendall_simulation.Rd new file mode 100644 index 0000000..8601b9c --- /dev/null +++ b/man/print.kendall_simulation.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/simulations_kendall_rw.R +\name{print.kendall_simulation} +\alias{print.kendall_simulation} +\title{Generic function that prints information about simulated Kendall random walk} +\usage{ +\method{print}{kendall_simulation}(x, ...) +} +\arguments{ +\item{x}{Object returned by simulate_kendall_rw or transform_kendall_rw function.} + +\item{...}{Other arguments.} +} +\description{ +Generic function that prints information about simulated Kendall random walk +} diff --git a/man/qkend.Rd b/man/qkend.Rd new file mode 100644 index 0000000..12603fb --- /dev/null +++ b/man/qkend.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stable_kendall_distribution.R +\name{qkend} +\alias{qkend} +\title{Quantiles of Kendall stable distribution} +\usage{ +qkend(m_alpha) +} +\arguments{ +\item{m_alpha}{function giving moments of order alpha of step dist.} +} +\value{ +function function returning quantiles of given orders +} +\description{ +Quantiles of Kendall stable distribution +} +\examples{ +qKend <- qkend(function(x) 1) +# Step distribution: delta_{1} +qKendall <- qKend(c(0.1, 0.9), 0.5) +# Quantiles of order 0.1 and 0.9 for alpha = 0.5 + + +} diff --git a/man/qkendSym.Rd b/man/qkendSym.Rd new file mode 100644 index 0000000..058b6ec --- /dev/null +++ b/man/qkendSym.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stable_kendall_distribution.R +\name{qkendSym} +\alias{qkendSym} +\title{Quantiles of symmetrical Kendall stable distribution} +\usage{ +qkendSym(m_alpha) +} +\arguments{ +\item{m_alpha}{function giving moments of order alpha of step dist.} +} +\value{ +function function returning quantiles of given orders +} +\description{ +Quantiles of symmetrical Kendall stable distribution +} +\examples{ +qKend <- qkendSym(function(x) 1) +# Step distribution: delta_{1} +qKendall <- qKend(c(0.1, 0.9), 0.5) +# Quantiles of order 0.1 and 0.9 for alpha = 0.5 + + +} diff --git a/man/rkend.Rd b/man/rkend.Rd new file mode 100644 index 0000000..353cb82 --- /dev/null +++ b/man/rkend.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stable_kendall_distribution.R +\name{rkend} +\alias{rkend} +\title{Pseudo-random number from Kendall stable distribution} +\usage{ +rkend(m_alpha) +} +\arguments{ +\item{m_alpha}{function giving moments of order alpha of step dist.} +} +\value{ +function return n numbers genereted from Kendall stable dist. +} +\description{ +Pseudo-random number from Kendall stable distribution +} +\examples{ +rKend <- rkend(function(x) 1) +# Step distribution: delta_{1} +rKendall <- rKend(10, 0.5) +# Ten random number from stable Kendall distribution with alpha = 0.5 + + +} diff --git a/man/simulateOneTrajectory.Rd b/man/simulateOneTrajectory.Rd new file mode 100644 index 0000000..a01f209 --- /dev/null +++ b/man/simulateOneTrajectory.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/simulations_kendall_rw.R +\name{simulateOneTrajectory} +\alias{simulateOneTrajectory} +\title{Simulate one trajectory ofa Kendall random walk} +\usage{ +simulateOneTrajectory(trajectory_length, step_dist, alpha, symmetric = FALSE, + ...) +} +\arguments{ +\item{trajectory_length}{Number of samples to simulate.} + +\item{step_dist}{Function that returns random numbers from step distribution.} + +\item{alpha}{Alpha parameter of the random walk} + +\item{symmetric}{If TRUE, random walk on the whole real line will be simulated.} + +\item{...}{Additional parameters to step distribution.} +} +\value{ +Generated path of the random walk. +} +\description{ +Simulate one trajectory ofa Kendall random walk +} diff --git a/man/simulate_kendall_rw.Rd b/man/simulate_kendall_rw.Rd new file mode 100644 index 0000000..127c70d --- /dev/null +++ b/man/simulate_kendall_rw.Rd @@ -0,0 +1,40 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/simulations_kendall_rw.R +\name{simulate_kendall_rw} +\alias{simulate_kendall_rw} +\title{Simulate multiple trajectories of Kendall random walk} +\usage{ +simulate_kendall_rw(number_of_simulations, trajectory_length, step_dist, alpha, + symmetric = FALSE, ...) +} +\arguments{ +\item{number_of_simulations}{number of trajectories to generate.} + +\item{trajectory_length}{length of trajectories.} + +\item{step_dist}{function returning random numbers from step dist.} + +\item{alpha}{alpha parameter.} + +\item{symmetric}{If TRUE, random walk on the whole real line will be simulated.} + +\item{...}{parameters for step distribution.} +} +\value{ +Object of class kendall_simulation. It is a list that consists of +\item{simulation}{Tibble with simulation id and simulated values,} +\item{step_distribution}{Name of the step distribution,} +\item{alpha}{Value of alpha parameter,} +\item{is_symmetric}{Logical value indicating if this is a symmetric Kendall R.W.} +} +\description{ +Object returned by this has print and plot methods. +} +\examples{ +kendall_simulations <- simulate_kendall_rw(10, 1000, runif, 0.5) +# Kendall R.W. on positive half-line with uniform step distribution - 10 trajectories. +only_simulations <- kendall_simulations$simulation # tibble with simulated values +kendall_simulations + + +} diff --git a/man/summarise_kendall_rw.Rd b/man/summarise_kendall_rw.Rd new file mode 100644 index 0000000..638a9f1 --- /dev/null +++ b/man/summarise_kendall_rw.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/simulations_kendall_rw.R +\name{summarise_kendall_rw} +\alias{summarise_kendall_rw} +\title{Calculate some characteristic for every simulated instance.} +\usage{ +summarise_kendall_rw(simulations, summary_function) +} +\arguments{ +\item{simulations}{Object of class kendall_simulation.} + +\item{summary_function}{Function that will be applied to each trajectory.} +} +\value{ +data frame or a list (of class kendall_simulation) +} +\description{ +Calculate some characteristic for every simulated instance. +} diff --git a/man/transform_kendall_rw.Rd b/man/transform_kendall_rw.Rd new file mode 100644 index 0000000..a79aa68 --- /dev/null +++ b/man/transform_kendall_rw.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/simulations_kendall_rw.R +\name{transform_kendall_rw} +\alias{transform_kendall_rw} +\title{Transforming (scaling and shifting) Kendall random walks} +\usage{ +transform_kendall_rw(simulations, an_seq = 1, bn_seq = 0) +} +\arguments{ +\item{simulations}{tibble returned by simulation function} + +\item{an_seq}{sequence that the trajectories will be multiplied by} + +\item{bn_seq}{sequence that will be substracted from scaled trajectory} +} +\value{ +List like in simulate_kendall_rw function after transforming trajectories. +} +\description{ +If one trajectory has length n, an_seq and bn_seq arguments should be sequnces of length n. +Object returned by this function has plot and print methods. +} +\examples{ +kendall_simulations <- simulate_kendall_rw(10, 1000, runif, 0.5) +scaled_kendall <- transform_kendall_rw(kendall_simulations, (1:1000)^(-2)) +scaled_kendall # kendall random walked scaled by the sequence n^(-1/alpha) +scaled_data <- scaled_kendall$simulation # simulated values +plot(scaled_kendall) + +} diff --git a/tests/testthat.R b/tests/testthat.R new file mode 100644 index 0000000..125e27b --- /dev/null +++ b/tests/testthat.R @@ -0,0 +1,4 @@ +library(testthat) +library(kendallRandomWalks) + +test_check("kendallRandomWalks") diff --git a/tests/testthat/test_barriers.R b/tests/testthat/test_barriers.R new file mode 100644 index 0000000..2b88559 --- /dev/null +++ b/tests/testthat/test_barriers.R @@ -0,0 +1,63 @@ +context("Test functions related to barrier crossing") + +set.seed(17) +kendall_rw <- simulate_kendall_rw(1000, 100, runif, 0.75) +symmetric_kendall_rw <- simulate_kendall_rw(1000, 100, rnorm, 0.25) + +moment <- ladder_moment(kendall_rw, 100) +moment2 <- ladder_moment(kendall_rw, 200) +symmetric_moment <- ladder_moment(symmetric_kendall_rw, 1000) + +height <- ladder_height(kendall_rw, 100) +height2 <- ladder_height(kendall_rw, 200) +symmetric_height <- ladder_height(symmetric_kendall_rw, 1000) + +testthat::test_that("Level must be positive", { + testthat::expect_error(ladder_moment(kendall_rw, -10)) + testthat::expect_error(ladder_height(kendall_rw, -10)) +}) + +testthat::test_that("Object have the right type", { + testthat::expect_is(symmetric_moment, "tbl_df") + testthat::expect_is(symmetric_moment, "kendall_barrier_crossing") + + testthat::expect_is(moment, "tbl_df") + testthat::expect_is(moment, "kendall_barrier_crossing") + + testthat::expect_is(symmetric_height, "tbl_df") + testthat::expect_is(symmetric_height, "kendall_barrier_crossing") + + testthat::expect_is(height, "tbl_df") + testthat::expect_is(height, "kendall_barrier_crossing") +}) + +testthat::test_that("Object have the right size", { + testthat::expect_equal(nrow(moment), 1000) + testthat::expect_equal(nrow(symmetric_moment), 1000) + testthat::expect_equal(nrow(height), 1000) + testthat::expect_equal(nrow(symmetric_height), 1000) +}) + +testthat::test_that("Bigger level implies greater mean of the distribution", { + testthat::expect_gte(mean(moment2$ladder_moment, na.rm = T), + mean(moment$ladder_moment, na.rm = T)) + testthat::expect_gte(mean(height2$ladder_moment, na.rm = T), + mean(height$ladder_moment, na.rm = T)) + +}) + +testthat::test_that("S3 methods are fine", { + testthat::expect_silent(plot(moment)) + testthat::expect_output(print(height)) +}) + +testthat::test_that("G function types are correct", { + testthat::expect_is(g_function_single(0.25, 0.77, dnorm), "integrate") + testthat::expect_is(g_function(1:5, 0.24, dnorm), "numeric") +}) + +testthat::test_that("Ladder moment PMF is okay based on level = 0", { + testthat::expect_equal(ladder_moment_pmf(2, 0, 0.75, pnorm, dnorm), 0.25) +}) + + diff --git a/tests/testthat/test_dists.R b/tests/testthat/test_dists.R new file mode 100644 index 0000000..0a8eb69 --- /dev/null +++ b/tests/testthat/test_dists.R @@ -0,0 +1,48 @@ +context("Check function related to stable Kendall distribution") + +pKend <- pkend(function(x) 1) +qKend <- qkend(function(x) 1) +pKendSym <- pkendSym(function(x) 1) +qKendSym <- qkendSym(function(x) 1) +rKend <- rkend(function(x) 1) +simulated <- rKend(10, 0.5) +dKend <- dkend(function(x) 1) +dens <- dKend(runif(10, 0, 10), 0.56) + +testthat::test_that("Support is on positive numbers", { + testthat::expect_equal(pKend(c(-20, -10, 0), 0.5), rep(0, 3)) +}) + +testthat::test_that("Quantiles are computed correctly", { + testthat::expect_equal(pKend(qKend(c(0.1, 0.5, 0.9), 0.5), 0.5), c(0.1, 0.5, 0.9), + tolerance = 1e-4) +}) + +testthat::test_that("Symmetric distribution is in fact symmetric", { + testthat::expect_equal(pKendSym(-2, 0.5), 1 - pKendSym(2, 0.5), + tolerance = 1e-4) +}) + + +testthat::test_that("Quantiles of symmetric distribution are computed correctly", { + testthat::expect_equal(pKendSym(qKendSym(c(0.1, 0.5, 0.9), 0.5), 0.5), + c(0.1, 0.5, 0.9), + tolerance = 1e-4) + +}) + +testthat::test_that("Numbers are simulated", { + testthat::expect_equal(length(simulated), 10) + testthat::expect_equal(sum(is.na(simulated)), 0) +}) + +testthat::test_that("Density is okay", { + testthat::expect_equal(sum(is.na(dens)), 0) + testthat::expect_equal(dKend(0, 0.76), 0) + testthat::expect_equal(dKend(1, 0.8), 0.8/exp(1)) +}) + +testthat::test_that("Non-finite values are treated correctly", { + testthat::expect_true(all(is.na(c(pKend(NA, 0.5), qKend(NA, 0.5), pKendSym(NA, 0.5), + qKendSym(NA, 0.5), dKend(NA, 0.5))))) +}) diff --git a/tests/testthat/test_fit.R b/tests/testthat/test_fit.R new file mode 100644 index 0000000..8e3a2a5 --- /dev/null +++ b/tests/testthat/test_fit.R @@ -0,0 +1,38 @@ +context("Test fit-related functions") + +set.seed(17) +rKend <- rkend(function(x) 1) +some_data <- rKend(100, 0.7) + +one_param <- estimate_stable_alpha(some_data) +one_param_fit <- fit_stable_alpha(some_data) + +loglik_more <- full_minus_loglik(some_data) +loglik_grad <- full_loglik_gradient(some_data) + +full_fit <- fit_kendall(some_data) + +testthat::test_that("Fit for one parameter is okay", { + testthat::expect_is(one_param$maximum, "numeric") + testthat::expect_equal(one_param$maximum, 0.7, tolerance = 0.1) +}) + +testthat::test_that("QQ-plot works", { + testthat::expect_output(plot(one_param_fit), regexp = NA) +}) + +testthat::test_that("Full minus loglik does not fail", { + testthat::expect_output(loglik_more(c(0.7, 0, 1)), regexp = NA) + testthat::expect_equal(is.na(loglik_more(c(0.7, 0, 1))), FALSE) + testthat::expect_equal(any(is.na(loglik_grad(c(0.7, 0, 1)))), FALSE) +}) + +testthat::test_that("Fit happened", { + testthat::expect_is(full_fit, "kendall_fit") + testthat::expect_equal(all(!is.na(unlist(full_fit$params))), TRUE) +}) + +testthat::test_that("Partial fit works", { + testthat::expect_output(fit_separate(some_data, 0.5), regexp = NA) + testthat::expect_error(fit_separate(rep(NA, 5))) +}) diff --git a/tests/testthat/test_simulations.R b/tests/testthat/test_simulations.R new file mode 100644 index 0000000..1e72c9f --- /dev/null +++ b/tests/testthat/test_simulations.R @@ -0,0 +1,44 @@ +context("Check function related to Kendall random walk simulations") + +kendall_rw <- simulate_kendall_rw(10, 100, runif, 0.75) +scaled_kendall <- transform_kendall_rw(kendall_rw, (1:100)^(-1/0.75)) + +symmetric_kendall_rw <- simulate_kendall_rw(10, 100, rnorm, 0.25, symmetric = T) +scaled_symmetric <- transform_kendall_rw(symmetric_kendall_rw, (1:100)^(-1/0.25)) + +testthat::test_that("Objects have the rights class", { + testthat::expect_is(kendall_rw, "kendall_simulation") + testthat::expect_is(scaled_kendall, "kendall_simulation") + testthat::expect_is(scaled_kendall$simulation, "tbl_df") +}) + +testthat::test_that("Number of simulations and observations is correct", { + testthat::expect_equal(max(kendall_rw$simulation$sim_id), 10) + testthat::expect_equal(nrow(kendall_rw$simulation), 10*100) +}) + +testthat::test_that("R.W. on positive half-line is positive, symmetric is not", { + testthat::expect_equal(all(kendall_rw$simulation$sim >= 0), TRUE) + testthat::expect_equal(any(symmetric_kendall_rw$simulation$sim < 0), TRUE) +}) + +testthat::test_that("Scaling does the right thing", { + testthat::expect_equal(any(kendall_rw$simulation$sim != scaled_kendall$simulation$sim), TRUE) + testthat::expect_equal(scaled_kendall$simulation$sim/rep((1:100)^(-1/0.75), times = 10), kendall_rw$simulation$sim) +}) + + +testthat::test_that("S3 methods are fine", { + testthat::expect_silent(plot(symmetric_kendall_rw)) + testthat::expect_silent(plot(symmetric_kendall_rw, level = 100)) + testthat::expect_output(print(symmetric_kendall_rw)) + testthat::expect_silent(plot(symmetric_kendall_rw, max_x = 100)) + testthat::expect_silent(plot(symmetric_kendall_rw, max_id = 1)) +}) + +testthat::test_that("Summary/mutate function does the right thing", { + testthat::expect_equal(unique(summarise_kendall_rw(symmetric_kendall_rw, length)$aggregated), 100) + testthat::expect_silent(mutate_kendall_rw(symmetric_kendall_rw, function(x) x^2)) + testthat::expect_is(mutate_kendall_rw(symmetric_kendall_rw, function(x) x^2, F), "kendall_simulation") +}) + diff --git a/vignettes/behaviour.Rmd b/vignettes/behaviour.Rmd new file mode 100644 index 0000000..3fc1692 --- /dev/null +++ b/vignettes/behaviour.Rmd @@ -0,0 +1,81 @@ +--- +title: "Studying the behaviour of Kendall random walks" +author: "Mateusz Staniak" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Studying the behaviour of Kendall random walks} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + warning = F, + message = F, + collapse = TRUE, + comment = "#>" +) +``` + +We will approximate the distribution of moments when the random walk changes state through simulations. + +First, we simulate many paths of Kendall random walk with normal step distribution. + +```{r sim} +library(kendallRandomWalks) +library(dplyr) +library(ggplot2) +set.seed(17) +walks <- simulate_kendall_rw(1000, 1000, rnorm, 0.5, T) +``` + +Example trajectory + +```{r ex1} +plot(walks, max_id = 1) +``` + +Number of unique states + +```{r unique} +ggplot(summarise_kendall_rw(walks, n_distinct), aes(x = aggregated)) + + # geom_histogram() + + theme_bw() + + geom_density() +``` + +Jumps + +```{r jumps} +diffs <- mutate_kendall_rw(walks, function(x) x - lag(x)) +diffs +diffs2 <- mutate_kendall_rw(walks, function(x) x - lag(x), F) +plot(diffs2, max_id = 10) +``` + +Time with no change of state + +```{r skoki_dlugosc} +diffs3 <- mutate_kendall_rw(diffs2, function(x) as.numeric(x != 0), F) +diffs <- diffs3$simulation +diffs <- group_by(diffs, sim_id) %>% + mutate(id = 1:n()) +lengths <- diffs %>% + filter(sim != 0) %>% + mutate(previous = ifelse(is.na(lag(id)), 0, lag(id))) %>% + mutate(length = id - previous) +ggplot(subset(lengths, sim_id < 5), + aes(x = length, fill = as.factor(sim_id), group = as.factor(sim_id))) + + geom_density() + + theme_bw() + + ggtitle("Distribution of time with no state-change (by simulation)") +ggplot(lengths, aes(x = length)) + + geom_density() + + theme_bw() + + ggtitle("Distribution of time with no state-change (aggregated)") +ggplot(subset(lengths, sim_id < 10), aes(x = id, y = length, color = as.factor(sim_id))) + + geom_point() + + theme_bw() + + geom_line() + + ggtitle("Time with no state-change in time") +``` diff --git a/vignettes/kendall_rws.Rmd b/vignettes/kendall_rws.Rmd new file mode 100644 index 0000000..b7cc72b --- /dev/null +++ b/vignettes/kendall_rws.Rmd @@ -0,0 +1,67 @@ +--- +title: "Kendall random walks" +author: "Mateusz Staniak" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Kendall random walks} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + + +## Simulate and plot + +```{r example1} +# uzupelnic fig_captions, zeby dalo sie je zrobic +library(kendallRandomWalks) +kendall_rws <- simulate_kendall_rw(10, 100, runif, 0.25) +kendall_rws +plot(kendall_rws) +plot(simulate_kendall_rw(10, 100, rnorm, 0.76), level = 300) +``` + +Symmetric + +```{r example1.1} +kendall_rws_sym <- simulate_kendall_rw(10, 100, rnorm, 0.76, T) +kendall_rws_sym +plot(kendall_rws_sym) +``` + +## Barrier crossing + +```{r example2} +kendall_rws2 <- simulate_kendall_rw(1000, 100, runif, 0.25) +ladder_moments <- ladder_moment(kendall_rws2, 1000) +ladder_moments +plot(ladder_moments) + +ladder_heights <- ladder_height(kendall_rws2, 2000) +ladder_heights +plot(ladder_heights) +``` + +Exact first ladder moments distribution with $G(a)$ computed numerically. + +```{r example3} +y <- seq(10, 10000, 50) +ladders <- sapply(y, + function(x) + ladder_moment_pmf(10, x, 0.5, pnorm, dnorm)) +plot(y, ladders) + +y <- seq(2000, 2200, 1) +plot(y, g_function(y, 0.1, dnorm)) + +plot(seq(0, 400, by = 1), g_function(seq(0, 400, by = 1), 0.5, dunif)) +``` + +