library(recurrentR)
set.seed(1)
lambda <- function(x) exp(sin(x))
T_0 <- rpois(1, 40)
curve(lambda, from = 0, to = T_0)
y <- rpois(nrow(iris), T_0)
y <- as.numeric(ifelse(y < T_0, y, T_0))
t <- sapply(y, function(y) {
as.numeric(ceiling(gen_inhomo_poisson(lambda, y - 1, 4)))
})
obj <- new("recurrent-data", iris, y, t, data.frame(), T_0)
rm(list = c("y", "t", "T_0"), envir = globalenv())
F.hat <- obj$F.hat
Lambda.hat <- obj$Lambda.hat
curve(F.hat, 0, obj@T_0)
curve(Lambda.hat, 0, obj@T_0)
curve(lambda, 0, obj@T_0)
Lambda.single <- function(t) integrate(lambda, 0, t)$value
Lambda <- function(T_0) {
return(function(t) sapply(t, Lambda.single))
}
answer <- Lambda(obj@T_0)
curve(answer, 0, obj@T_0, col = 2)
curve(Lambda.hat, 0, obj@T_0, add = TRUE, lty = 1)
x.eval <- seq(from = 0, to = obj@T_0, length = 100)
obj$F.hat(x.eval, TRUE, error.measurement.function = sd)
## $estimate
## [1] 0.00000 0.00000 0.00000 0.03662 0.03662 0.03662 0.09111 0.09111
## [9] 0.09111 0.13393 0.13393 0.14914 0.14914 0.14914 0.15637 0.15637
## [17] 0.15637 0.16804 0.16804 0.16804 0.19742 0.19742 0.25072 0.25072
## [25] 0.25072 0.29871 0.29871 0.29871 0.32071 0.32071 0.32071 0.33031
## [33] 0.33031 0.34094 0.34094 0.34094 0.36354 0.36354 0.36354 0.41079
## [41] 0.41079 0.41079 0.47044 0.47044 0.49879 0.49879 0.49879 0.50854
## [49] 0.50854 0.50854 0.51769 0.51769 0.51769 0.53379 0.53379 0.57735
## [57] 0.57735 0.57735 0.63331 0.63331 0.63331 0.66860 0.66860 0.66860
## [65] 0.68115 0.68115 0.69016 0.69016 0.69016 0.70301 0.70301 0.70301
## [73] 0.73176 0.73176 0.73176 0.78632 0.78632 0.82785 0.82785 0.82785
## [81] 0.84406 0.84406 0.84406 0.85341 0.85341 0.85341 0.86338 0.86338
## [89] 0.88453 0.88453 0.88453 0.93504 0.93504 0.93504 0.98040 0.98040
## [97] 0.98040 1.00000 1.00000 1.00000
##
## $error.measurement
## [1] 0.000000 0.000000 0.000000 0.002201 0.002201 0.002201 0.003545
## [8] 0.003545 0.003545 0.004068 0.004068 0.004532 0.004532 0.004532
## [15] 0.004763 0.004763 0.004763 0.004614 0.004614 0.004614 0.004440
## [22] 0.004440 0.004735 0.004735 0.004735 0.004563 0.004563 0.004563
## [29] 0.005081 0.005081 0.005081 0.004993 0.004993 0.005144 0.005144
## [36] 0.005144 0.005432 0.005432 0.005432 0.005534 0.005534 0.005534
## [43] 0.005478 0.005478 0.006043 0.006043 0.006043 0.006013 0.006013
## [50] 0.006013 0.006050 0.006050 0.006050 0.006063 0.006063 0.006380
## [57] 0.006380 0.006380 0.006043 0.006043 0.006043 0.005930 0.005930
## [64] 0.005930 0.005725 0.005725 0.005573 0.005573 0.005573 0.005797
## [71] 0.005797 0.005797 0.005496 0.005496 0.005496 0.005435 0.005435
## [78] 0.005270 0.005270 0.005270 0.005038 0.005038 0.005038 0.005192
## [85] 0.005192 0.005192 0.005042 0.005042 0.005150 0.005150 0.005150
## [92] 0.003950 0.003950 0.003950 0.002233 0.002233 0.002233 0.000000
## [99] 0.000000 0.000000
obj$F.hat(x.eval, TRUE, error.measurement.function = function(a) c(quantile(a,
0.975), quantile(a, 0.025)))
## $estimate
## [1] 0.00000 0.00000 0.00000 0.03662 0.03662 0.03662 0.09111 0.09111
## [9] 0.09111 0.13393 0.13393 0.14914 0.14914 0.14914 0.15637 0.15637
## [17] 0.15637 0.16804 0.16804 0.16804 0.19742 0.19742 0.25072 0.25072
## [25] 0.25072 0.29871 0.29871 0.29871 0.32071 0.32071 0.32071 0.33031
## [33] 0.33031 0.34094 0.34094 0.34094 0.36354 0.36354 0.36354 0.41079
## [41] 0.41079 0.41079 0.47044 0.47044 0.49879 0.49879 0.49879 0.50854
## [49] 0.50854 0.50854 0.51769 0.51769 0.51769 0.53379 0.53379 0.57735
## [57] 0.57735 0.57735 0.63331 0.63331 0.63331 0.66860 0.66860 0.66860
## [65] 0.68115 0.68115 0.69016 0.69016 0.69016 0.70301 0.70301 0.70301
## [73] 0.73176 0.73176 0.73176 0.78632 0.78632 0.82785 0.82785 0.82785
## [81] 0.84406 0.84406 0.84406 0.85341 0.85341 0.85341 0.86338 0.86338
## [89] 0.88453 0.88453 0.88453 0.93504 0.93504 0.93504 0.98040 0.98040
## [97] 0.98040 1.00000 1.00000 1.00000
##
## $error.measurement
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## 97.5% 0 0 0 0.04168 0.04168 0.04168 0.09739 0.09739 0.09739
## 2.5% 0 0 0 0.03261 0.03261 0.03261 0.08453 0.08453 0.08453
## [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18]
## 97.5% 0.1421 0.1421 0.1567 0.1567 0.1567 0.1645 0.1645 0.1645 0.1753
## 2.5% 0.1270 0.1270 0.1418 0.1418 0.1418 0.1484 0.1484 0.1484 0.1607
## [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27]
## 97.5% 0.1753 0.1753 0.2062 0.2062 0.2607 0.2607 0.2607 0.3084 0.3084
## 2.5% 0.1607 0.1607 0.1904 0.1904 0.2425 0.2425 0.2425 0.2891 0.2891
## [,28] [,29] [,30] [,31] [,32] [,33] [,34] [,35] [,36]
## 97.5% 0.3084 0.3316 0.3316 0.3316 0.3417 0.3417 0.3534 0.3534 0.3534
## 2.5% 0.2891 0.3104 0.3104 0.3104 0.3193 0.3193 0.3286 0.3286 0.3286
## [,37] [,38] [,39] [,40] [,41] [,42] [,43] [,44] [,45]
## 97.5% 0.3752 0.3752 0.3752 0.4249 0.4249 0.4249 0.4839 0.4839 0.5103
## 2.5% 0.3527 0.3527 0.3527 0.3997 0.3997 0.3997 0.4589 0.4589 0.4861
## [,46] [,47] [,48] [,49] [,50] [,51] [,52] [,53] [,54]
## 97.5% 0.5103 0.5103 0.5205 0.5205 0.5205 0.5283 0.5283 0.5283 0.5452
## 2.5% 0.4861 0.4861 0.4975 0.4975 0.4975 0.5080 0.5080 0.5080 0.5241
## [,55] [,56] [,57] [,58] [,59] [,60] [,61] [,62] [,63]
## 97.5% 0.5452 0.5906 0.5906 0.5906 0.6452 0.6452 0.6452 0.6810 0.6810
## 2.5% 0.5241 0.5673 0.5673 0.5673 0.6213 0.6213 0.6213 0.6569 0.6569
## [,64] [,65] [,66] [,67] [,68] [,69] [,70] [,71] [,72]
## 97.5% 0.6810 0.6928 0.6928 0.7005 0.7005 0.7005 0.7130 0.7130 0.7130
## 2.5% 0.6569 0.6696 0.6696 0.6789 0.6789 0.6789 0.6929 0.6929 0.6929
## [,73] [,74] [,75] [,76] [,77] [,78] [,79] [,80] [,81]
## 97.5% 0.7409 0.7409 0.7409 0.7969 0.7969 0.8404 0.8404 0.8404 0.8542
## 2.5% 0.7207 0.7207 0.7207 0.7742 0.7742 0.8171 0.8171 0.8171 0.8335
## [,82] [,83] [,84] [,85] [,86] [,87] [,88] [,89] [,90]
## 97.5% 0.8542 0.8542 0.8638 0.8638 0.8638 0.8736 0.8736 0.8946 0.8946
## 2.5% 0.8335 0.8335 0.8435 0.8435 0.8435 0.8530 0.8530 0.8755 0.8755
## [,91] [,92] [,93] [,94] [,95] [,96] [,97] [,98] [,99] [,100]
## 97.5% 0.8946 0.9417 0.9417 0.9417 0.9855 0.9855 0.9855 1 1 1
## 2.5% 0.8755 0.9279 0.9279 0.9279 0.9762 0.9762 0.9762 1 1 1
result <- obj$Lambda.hat(x.eval, TRUE, error.measurement.function = sd)
curve(answer, 0, obj@T_0, col = 2)
lines(x.eval, result$estimate, col = 1)
lines(x.eval, result$estimate + 2 * result$error.measurement, lty = 2, col = 1)
lines(x.eval, result$estimate - 2 * result$error.measurement, lty = 2, col = 1)
result <- obj$Lambda.hat(x.eval, TRUE, error.measurement.function = function(a) c(quantile(a,
0.975), quantile(a, 0.025)))
curve(answer, 0, obj@T_0, col = 2)
lines(x.eval, result$estimate, col = 1)
lines(x.eval, result$error.measurement[1, ], lty = 2, col = 1)
lines(x.eval, result$error.measurement[2, ], lty = 2, col = 1)
