Commit 1ba11889 authored by Sebastian Schweer's avatar Sebastian Schweer

Close to first release, API defined, Check not clean yet

parent 261ea9a2
Package: queueingnetworkR
Title: What the Package Does (one line, title case)
Title: Simulation and Estimation of a Geom/Geom/infty Tandem Network
Version: 0.0.1
[email protected]: person("Sebastian", "Schweer", email = "[email protected]", role = c("aut", "cre"))
Description: Simulation and Documentation for the Paper "Nonparametric Estimation of the Service Time Distribution in Discrete-Time Queueing Networks" by Cornelia Wichelhaus and Sebastian Schweer
Depends: R (>= 3.4.2), utils, stats
Depends: R (>= 3.4.2), utils, stats, ggplot2, plyr, reshape2
License: MIT
Encoding: UTF-8
LazyData: true
......
# Generated by roxygen2: do not edit by hand
export(present_estimates)
export(queueing_network_poi_geom)
importFrom(ggplot2,aes)
importFrom(ggplot2,geom_line)
importFrom(ggplot2,ggplot)
importFrom(plyr,adply)
importFrom(reshape2,melt)
importFrom(stats,dgeom)
importFrom(stats,pgeom)
importFrom(stats,rpois)
importFrom(utils,head)
importFrom(utils,setTxtProgressBar)
......
cross_covariance <- function(data, node_1, node_2, lag, nobs){
if(is.null(nobs)){
nobs = nrow(data)
}
if(!is.null(nobs) & nobs > nrow(data)){
nobs = nrow(data)
}
if(lag < nobs){
if(node_1 == 1){
arrival = data$arrival_1
} else if (node_1 ==2){
arrival = data$arrival_2
} else {
message("Arrival node not contained in the model")
}
if(node_2 == 1){
departure = data$departure_1
} else if (node_2 ==2){
departure = data$departure_2
} else {
message("Departure node not contained in the model")
}
arrival_k = head(arrival, (nobs - lag))
departure_k = tail(departure, (nobs - lag))
### Following formula (3) in Schweer/Wichelhaus
estimator = 1/nobs*(sum(arrival_k*departure_k)) - 1/nobs*sum(arrival_k)*mean(departure) - 1/nobs*sum(departure_k)*mean(arrival) + (nobs - lag)/nobs*mean(arrival)*mean(departure)
} else {
estimator = 0
}
return(estimator)
}
# Following Lemma 9 in Edelmann/Wichelhaus
deconvoluting_alpha <- function(number_of_steps, alpha_1, alpha_2, p_12, p_21){
M_n = min(which(alpha_1 > 0))
if(is.finite(M_n)){
length_1 = length(alpha_1) - M_n + 1
length_2 = length(alpha_2) - M_n
max_num_steps = min(length_1, length_2, number_of_steps)
if(max_num_steps > 1){
v_n = numeric(length = max_num_steps)
v_n[1] = (1 - p_12)* alpha_2[M_n + 1] / (alpha_1[M_n] * (1 - p_21) * p_12)
for(idx in c(2:max_num_steps)){
alpha_1_k = alpha_1[(M_n + idx -1):(M_n + 1)]
v_n[idx] = (1 - p_12)* alpha_2[M_n + idx] / (alpha_1[M_n] * (1 - p_21) * p_12) - (sum(alpha_1_k*v_n[1:(idx-1)]))/alpha_1[M_n]
}
return(v_n)
} else {
message("Not enough data points for deconvolution")
}
} else {
message("Alpha Estimator consists of 0 values")
}
cross_covariance <- function(data, node_1, node_2, lag, nobs){
if(is.null(nobs)){
nobs = nrow(data)
}
if(!is.null(nobs) & nobs > nrow(data)){
nobs = nrow(data)
}
if(lag < nobs){
if(node_1 == 1){
arrival = data$arrival_1
} else if (node_1 ==2){
arrival = data$arrival_2
} else {
message("Arrival node not contained in the model")
}
if(node_2 == 1){
departure = data$departure_1
} else if (node_2 ==2){
departure = data$departure_2
} else {
message("Departure node not contained in the model")
}
arrival_k = head(arrival, (nobs - lag))
departure_k = tail(departure, (nobs - lag))
### Following formula (3) in Schweer/Wichelhaus
estimator = 1/nobs*(sum(arrival_k*departure_k)) - 1/nobs*sum(arrival_k)*mean(departure) - 1/nobs*sum(departure_k)*mean(arrival) + (nobs - lag)/nobs*mean(arrival)*mean(departure)
} else {
estimator = 0
}
return(estimator)
}
#' Deconvolute the Estimator to Obtain Service Time Distribution
#'
#' Implementation of the theoretical result of Lemma 9 in Edelmann, Wichelhaus
#'
#' @param number_of_steps Integer, number of steps in deconvolution.
#' @param alpha_1 Numeric vector, estimator "hat" alpha_n^2(...) in Lemma 9.
#' @param alpha_2 Numeric vector, estimator "hat" alpha_n^2(...) in Lemma 9.
#' @param p_12 Numeric, in [0,1], routing probability from node 1 to node 2
#' @param p_21 Numeric, in [0,1], routing probability from node 2 to node 1
#' @return Numeric vector , corresponding to "hat" v_n (...) in Lemma 9.
deconvoluting_alpha <- function(number_of_steps, alpha_1, alpha_2, p_12, p_21){
M_n = min(which(alpha_1 > 0))
if(is.finite(M_n)){
length_1 = length(alpha_1) - M_n + 1
length_2 = length(alpha_2) - M_n
max_num_steps = min(length_1, length_2, number_of_steps)
if(max_num_steps > 1){
v_n = numeric(length = max_num_steps)
v_n[1] = (1 - p_12)* alpha_2[M_n + 1] / (alpha_1[M_n] * (1 - p_21) * p_12)
for(idx in c(2:max_num_steps)){
alpha_1_k = alpha_1[(M_n + idx -1):(M_n + 1)]
v_n[idx] = (1 - p_12)* alpha_2[M_n + idx] / (alpha_1[M_n] * (1 - p_21) * p_12) - (sum(alpha_1_k*v_n[1:(idx-1)]))/alpha_1[M_n]
}
return(v_n)
} else {
message("Not enough data points for deconvolution")
}
} else {
message("Alpha Estimator consists of 0 values")
}
}
\ No newline at end of file
estimate_g2 <- function(max_lag, p_12, p_21, ...){
data = queueing_network_poi_geom(p_12 = p_12, p_21 = p_21, ...)
alpha_12_hat = sapply(c(1:max_lag),
cross_covariance,
data = data, node_1 = 1, node_2 = 2, nobs = NULL)
alpha_11_hat = sapply(c(1:max_lag),
cross_covariance,
data = data, node_1 = 1, node_2 = 1, nobs = NULL)
g_2_est = deconvoluting_alpha(max_lag,
alpha_1 = alpha_11_hat,
alpha_2 = alpha_12_hat,
p_12 = p_12,
p_21 = p_21)
return(g_2_est)
}
#' Helper Function to Apply Theorem 4 in Edelmann, Wichelhaus
#'
#' Creates a cdf estimator from a pgf estimator.
#'
#' @param g2_est Numeric vector, the pgf estimation of the service
#' time distribution
#'
estimate_G2 <- function(g2_est){
trunc_g2_est = ifelse(g2_est > 0, g2_est, 0)
G2_est_hat = cumsum(trunc_g2_est)
trunc_G2_est_hat = ifelse(G2_est_hat < 1, G2_est_hat, 1)
return(trunc_G2_est_hat)
}
\ No newline at end of file
#' Presents the Results of the Cross-Covariance Estimation
#'
#' Given a parametrization of a Tandem Netowrk, this funciton produces a number
#' of repetitions of simulation and estimation of the servive time G2 at the
#' second node, and presents the results as data frame and plot.
#'
#' @param n_reps Integer, number of repetitions of simulation and estimation to
#' be calculated.
#' @param max_lag Integer, length of the estimated service time distribution.
#' @param ... Parameters passed on to the function estimate_g2.
#' @return A list: result_df is the data frame containing the raw estimators with
#' a column rep_id for each repetition, delta_reulst_df contains the same result, with
#' each value normalized by the theoretical g_2, plot_result and plot_delta_result are
#' ggplots for the result data sets.
#' @importFrom ggplot2 ggplot geom_line aes
#' @importFrom plyr adply
#' @importFrom reshape2 melt
#' @importFrom stats dgeom pgeom
#' @export
present_estimates <- function(n_reps, max_lag, G_2, ...){
index_set = c(1:n_reps)
result_set = sapply(index_set,
function(x){cat("\nRun", x, "of", n_reps, "\n")
return(estimate_g2(max_lag, G_2, ...))},
simplify = TRUE)
result_df = as.data.frame(t(result_set))
colnames(result_df) = paste0("est_g2_", c(1:(max_lag - 1)))
result_G2_df = adply(result_df, 1, estimate_G2)[, c((max_lag):(2*(max_lag - 1)))]
colnames(result_G2_df) = paste0("est_G2_", c(1:(max_lag - 1)))
g_2_true = dgeom(c(1:(max_lag - 1)), G_2)
G_2_true = pgeom(c(1:(max_lag - 1)), G_2)
delta_result_df = sweep(result_df, 2, g_2_true)
delta_result_G2_df = sweep(result_G2_df, 2, G_2_true)
colnames(delta_result_df) = paste0("delta_", colnames(result_df))
colnames(delta_result_G2_df) = paste0("delta_", colnames(result_G2_df))
result_df$rep_id = as.factor(index_set)
delta_result_df$rep_id = as.factor(index_set)
result_G2_df$rep_id = as.factor(index_set)
delta_result_G2_df$rep_id = as.factor(index_set)
plot_result_G2_df = melt(result_G2_df, id.vars = 'rep_id')
plot_result_G2 = ggplot(plot_result_G2_df, aes(x = variable, y = value)) + geom_line(aes(color = rep_id, group = rep_id))
plot_delta_result_df = melt(delta_result_df, id.vars = 'rep_id')
plot_delta_result = ggplot(plot_delta_result_df, aes(x = variable, y = value)) + geom_line(aes(color = rep_id, group = rep_id))
return(list(result_df = result_df
, result_G2_df = result_G2_df
, delta_result_df = delta_result_df
, delta_result_G2_df = delta_result_G2_df
, plot_result_G2 = plot_result_G2
, plot_delta_result = plot_delta_result
))
}
\ No newline at end of file
simulation_estimation <- function(max_lag, p_12, p_21, ...){
data = data <- queueing_network_poi_geom(p_12 = p_12, p_21 = p_21, ...)
alpha_12_hat = sapply(c(1:max_lag),
queueingnetworkR:::cross_covariance,
data = data, node_1 = 1, node_2 = 2, nobs = NULL)
alpha_11_hat = sapply(c(1:max_lag),
queueingnetworkR:::cross_covariance,
data = data, node_1 = 1, node_2 = 1, nobs = NULL)
g_2_est = queueingnetworkR:::deconvoluting_alpha(max_lag,
alpha_1 = alpha_11_hat,
alpha_2 = alpha_12_hat,
p_12 = p_12,
p_21 = p_21)
return(g_2_est)
}
#simulation_estimation(10, 0.1, 0.9, n_obs = 100000, burn_in = 1000, lambda_1 = 1, lambda_2 = 0.3, G_1 = 0.3, G_2 = 0.9, progress = TRUE)
#
# p_12 = 0.7
# p_21 = 0.3
# max_lag = 12
#
# data <- queueing_network_poi_geom(p_12 = p_12, p_21 = p_21, n_obs = 10000, burn_in = 1000, lambda_1 = 0.3, lambda_2 = 0.5, G_1 = 0.3, G_2 = 0.1, progress = TRUE)
#
# alpha_12_hat = sapply(c(1:max_lag),
# cross_covariance,
# data = data, node_1 = 1, node_2 = 2, nobs = NULL)
# alpha_11_hat = sapply(c(1:max_lag),
# queueingnetworkR:::cross_covariance,
# data = data, node_1 = 1, node_2 = 1, nobs = NULL)
# g_2_est = deconvoluting_alpha(max_lag,
# alpha_1 = alpha_11_hat,
# alpha_2 = alpha_12_hat,
# p_12 = p_12,
# p_21 = p_21)
\ No newline at end of file
# Simulations and Documentation
## Description
Simulation and Documentation for the Paper "Nonparametric Estimation of the Service Time Distribution in Discrete-Time Queueing Networks"
## List of Files
### queueing<sub>network.R</sub>
Simulation of a two node QN.
### spa<sub>queueingnetwork.RProj</sub>
RProject used in RStudio.
### README.org
# Simulations and Documentation
## Description
Simulation and Documentation for the Paper "Nonparametric Estimation of the Service Time Distribution in Discrete-Time Queueing Networks"
## List of Files
### queueing<sub>network.R</sub>
Simulation of a two node QN.
### spa<sub>queueingnetwork.RProj</sub>
RProject used in RStudio.
### README.org
ORG-File in EMACS, creates README.md.
\ No newline at end of file
---
title: "README"
author: "Sebastian Schweer"
date: "September 15, 2018"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## Live demontration
Here is a running example to try out:
```{r}
p_12 = 0.1
p_21 = 0.2
n_obs = 1000
burn_in = 1000
lambda_1 = 1
lambda_2 = 0.3
G_1 = 0.3
G_2 = 0.9
test <- present_estimates(n_reps = 10, max_lag = 5, p_12 = p_12, p_21 = p_21, n_obs = n_obs, burn_in = burn_in,
lambda_1 = lambda_1, lambda_2 = lambda_2, G_1 = G_1, G_2 = G_2,
progress = TRUE)
```
This leads to plots:
```{r, echo=FALSE}
test$plot_result_G2
```
As well as
```{r, echo=FALSE}
test$plot_delta_result
```
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/cross_covariance_estimator.R
\name{deconvoluting_alpha}
\alias{deconvoluting_alpha}
\title{Deconvolute the Estimator to Obtain Service Time Distribution}
\usage{
deconvoluting_alpha(number_of_steps, alpha_1, alpha_2, p_12, p_21)
}
\arguments{
\item{number_of_steps}{Integer, number of steps in deconvolution.}
\item{alpha_1}{Numeric vector, estimator "hat" alpha_n^2(...) in Lemma 9.}
\item{alpha_2}{Numeric vector, estimator "hat" alpha_n^2(...) in Lemma 9.}
\item{p_12}{Numeric, in [0,1], routing probability from node 1 to node 2}
\item{p_21}{Numeric, in [0,1], routing probability from node 2 to node 1}
}
\value{
Numeric vector , corresponding to "hat" v_n (...) in Lemma 9.
}
\description{
Implementation of the theoretical result of Lemma 9 in Edelmann, Wichelhaus
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/estimate_g2.R
\name{estimate_G2}
\alias{estimate_G2}
\title{Helper Function to Apply Theorem 4 in Edelmann, Wichelhaus}
\usage{
estimate_G2(g2_est)
}
\arguments{
\item{g2_est}{Numeric vector, the pgf estimation of the service
time distribution}
}
\description{
Creates a cdf estimator from a pgf estimator.
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/present_estimates.R
\name{present_estimates}
\alias{present_estimates}
\title{Presents the Results of the Cross-Covariance Estimation}
\usage{
present_estimates(n_reps, max_lag, G_2, ...)
}
\arguments{
\item{n_reps}{Integer, number of repetitions of simulation and estimation to
be calculated.}
\item{max_lag}{Integer, length of the estimated service time distribution.}
\item{...}{Parameters passed on to the function estimate_g2.}
}
\value{
A list: result_df is the data frame containing the raw estimators with
a column rep_id for each repetition, delta_reulst_df contains the same result, with
each value normalized by the theoretical g_2, plot_result and plot_delta_result are
ggplots for the result data sets.
}
\description{
Given a parametrization of a Tandem Netowrk, this funciton produces a number
of repetitions of simulation and estimation of the servive time G2 at the
second node, and presents the results as data frame and plot.
}
context("test-simulation_estimation.R")
p_12 = 0.1
p_21 = 0.2
n_obs = 1000
burn_in = 1000
lambda_1 = 1
lambda_2 = 0.3
G_1 = 0.3
G_2 = 0.9
test_that("Function returns a numeric vector", {
expect_equal(class(estimate_g2(10, p_12 = p_12, p_21 = p_21, n_obs = n_obs, burn_in = burn_in,
lambda_1 = lambda_1, lambda_2 = lambda_2, G_1 = G_1, G_2 = G_2,
progress = FALSE)), "numeric")
expect_equal(class(estimate_g2(20, p_12 = p_12, p_21 = p_21, n_obs = n_obs, burn_in = burn_in,
lambda_1 = lambda_1, lambda_2 = lambda_2, G_1 = G_1, G_2 = G_2,
progress = FALSE)), "numeric")
expect_equal(class(estimate_g2(3, p_12 = p_12, p_21 = p_21, n_obs = n_obs, burn_in = burn_in,
lambda_1 = lambda_1, lambda_2 = lambda_2, G_1 = G_1, G_2 = G_2,
progress = FALSE)), "numeric")
})
test_that("Function returns a numeric vector of expected length" , {
expect_equal(length(estimate_g2(10, p_12 = p_12, p_21 = p_21, n_obs = n_obs, burn_in = burn_in,
lambda_1 = lambda_1, lambda_2 = lambda_2, G_1 = G_1, G_2 = G_2,
progress = FALSE)), 9)
expect_equal(length(estimate_g2(20, p_12 = p_12, p_21 = p_21, n_obs = n_obs, burn_in = burn_in,
lambda_1 = lambda_1, lambda_2 = lambda_2, G_1 = G_1, G_2 = G_2,
progress = FALSE)), 19)
expect_equal(length(estimate_g2(3, p_12 = p_12, p_21 = p_21, n_obs = n_obs, burn_in = burn_in,
lambda_1 = lambda_1, lambda_2 = lambda_2, G_1 = G_1, G_2 = G_2,
progress = FALSE)), 2)
})
context("Estimate G2")
test_that("Correct transformations" , {
expect_equal(estimate_G2(c(0.3, 0.2, 0.4, 0.8, 0.1)), c(0.3, 0.5, 0.9, 1, 1))
expect_equal(estimate_G2(c(0.3, 0.2, - 0.4, 0.8, 0.1)), c(0.3, 0.5, 0.5, 1, 1))
expect_equal(estimate_G2(c(4.3, 0.2, 0.4, 0.8, 0.1)), c(1, 1, 1, 1, 1))
})
\ No newline at end of file
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment