-
Notifications
You must be signed in to change notification settings - Fork 0
/
exponential_class.R
135 lines (127 loc) · 4.82 KB
/
exponential_class.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
## Survival objects for simulations
## by JJAV 20220420
#' Factory of SURVIVAL objects with Exponential distributions
#'
#' Creates a SURVIVAL object with an Exponential distribution.
#'
#' @section Parameters:
#'
#' To create an exponential survival object the following
#' options are available:
#'
#' _`lambda`_ to specify the canonical parameter of the distribution, or
#'
#' _`surv`_ and _`t`_ for the proportion surviving (no events) at time t, or
#'
#' _`fail`_ and _`t`_ for the proportion failing (events) at time t
#'
#' lambda = -log(surv)/t
#'
#' lambda = -log(1-fail)/t
#'
#' The parameters should be spell correctly as partial matching is not available
#'
#' @param ... Parameters to define the distribution. See the Parameters for details
#' @return a SURVIVAL object of the exponential distribution family. See the
#' documentation of `s_factory` for the methods available for SURVIVAL objects
#' @importFrom stats runif
#' @export
#' @examples
#' s_exponential(lambda = 3)
#' s_exponential(surv = 0.4, t = 2)
#' s_exponential(fail = 0.6, t = 2)
s_exponential <- function(...) {
params <- list(...)
nparam <- names(params)
# This function is the factory of the class
.factory_exponential <- function(lambda) {
iCum_Hfx <- function(H){
stopifnot("Must be positive number" = all(H >= 0))
H/lambda
}
structure(
list(
distribution = "EXPONENTIAL",
params = list(lambda = lambda),
sfx = function(t) {
stopifnot("t must be numeric" = is.numeric(t))
stopifnot("t must be positive number" = all(t >= 0))
exp(-lambda*t)
},
hfx = function(t) {
stopifnot("t must be numeric" = is.numeric(t))
stopifnot("t must be positive number" = all(t >= 0))
rep(lambda,length(t))
},
Cum_Hfx = function(t) {
stopifnot("t must be numeric" = is.numeric(t))
stopifnot("t must be positive number" = all(t >= 0))
lambda*t
},
invCum_Hfx=iCum_Hfx,
rsurv = function(n){
stopifnot("n must be numeric" = is.numeric(n))
stopifnot("n must be positive number" = all(n > 0))
iCum_Hfx(-log(runif(n)))
},
rsurvhr = function(hr){
stopifnot("hr must be numeric" = is.numeric(hr))
stopifnot("hr must be positive numbers > 0" = all(hr > 0))
# Following Bender, Augustin and Blettner 2005
iCum_Hfx(-log(runif(length(hr)))/hr)
},
rsurvaft = function(aft){
stopifnot("aft must be numeric" = is.numeric(aft))
stopifnot("aft must be positive numbers > 0" = all(aft > 0))
iCum_Hfx(-log(runif(length(aft))))/aft
},
rsurvah = function(aft,hr){
stopifnot("aft must be numeric" = is.numeric(aft))
stopifnot("hr must be numeric" = is.numeric(hr))
stopifnot("aft and hr must be of the same length" = length(aft)==length(hr) )
stopifnot("aft must be positive numbers > 0" = all(aft > 0))
stopifnot("hr must be positive numbers > 0" = all(hr > 0))
iCum_Hfx(-log(runif(length(aft)))/hr)/aft
}
),
class = c("SURVIVAL")
)
}
# Definition based on lambda
if (length(nparam == 1) &
("lambda" %in% nparam)) {
stopifnot("lambda should be a single number" = is_single_number(params$lambda) )
stopifnot("lambda should be greater than 0" = params$lambda > 0 )
return(.factory_exponential(params$lambda))
}
# Definition based in proportion surviving and time
if(
length(nparam == 2) &
all(c("surv","t") %in% nparam)) {
stopifnot("surv must be a single number" = is_single_number(params$surv))
stopifnot("surv must be greater than 0" = params$surv > 0)
stopifnot("surv must be smaller than 1" = params$surv < 1)
stopifnot("t must be a single number" = is_single_number(params$t))
stopifnot("t must be greater than 0" = params$t > 0)
lambda = -log(params$surv) / params$t
return(.factory_exponential(lambda))
}
# Definition based on proportion failing and time
if(
length(nparam == 2) &
all(c("fail","t") %in% nparam)) {
stopifnot("fail must be a single number" = is_single_number(params$fail))
stopifnot("fail must be greater than 0" = params$fail > 0)
stopifnot("fail must be lower than 1" = params$fail < 1)
stopifnot("t must be a single number" = is_single_number(params$t))
stopifnot("t must be greater than 0" = params$t > 0)
lambda = -log(1 -params$fail)/params$t
return(.factory_exponential(lambda))
}
message(
"Valid parameters to define a Exponential distribution are: \n",
"lambda: for the canonical parameter of the distribution, or\n",
"surv, t: for the surviving proportion (no events) at time t, or\n",
"fail, t: for the failure proportion (events) at time t \n")
stop("Not valid parameters")
}