/
dWeibull.R
207 lines (189 loc) · 8.4 KB
/
dWeibull.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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
#' Probability calculations for Weibull count models
#'
#' Probability computations for the univariate Weibull count process. Several
#' methods are provided.
#' \code{dWeibullCount} computes probabilities.
#'
#' Argument \code{method} can be used to specify the desired method, as follows:
#' \describe{
#' \item{\code{"series_mat"}}{- series expansion using matrix techniques,}
#' \item{\code{"series_acc"}}{- Euler-van Wijngaarden accelerated series expansion (default),}
#' \item{\code{"conv_direc"t}}{- direct convolution method of section 2,}
#' \item{\code{"conv_naive"}}{- naive convolurion described in section 3.1,}
#' \item{\code{"conv_dePril"}}{- dePril convolution described in section 3.2.}
#' }
#' The arguments have sensible default values.
#'
#' @param scale numeric (length 1), scale parameter of the Weibull count.
#' @param shape numeric (length 1), shape parameter of the Weibull count.
#' @param x integer (vector), the desired count values.
#' @param time double, length of the observation window (defaults to 1).
#' @param log logical, if TRUE, the log of the probability will be returned.
#' @param method character, one of the available methods, see details.
#' @param conv_steps numeric, number of steps used for the extrapolation.
#' @param conv_extrap logical, should Richardson extrappolation be applied ?
#' @param series_terms numeric, number of terms in the series expansion.
#' @param series_acc_niter numeric, number of iterations in the
#' Euler-van Wijngaarden algorithm.
#' @param series_acc_eps numeric, tolerance of convergence in the
#' Euler-van Wijngaarden algorithm.
#' @return for \code{dWeibullCount}, a vector of probabilities
#' \eqn{P(x(i)), i = 1, \dots n}, where \eqn{n =} \code{length(x)}.
#' @export
dWeibullCount <- function(x, shape, scale,
method = c("series_acc", "series_mat",
"conv_direct", "conv_naive", "conv_dePril"),
time = 1, log = FALSE, conv_steps = 100,
conv_extrap = TRUE, series_terms = 50,
series_acc_niter = 300, series_acc_eps = 1e-10) {
method <- match.arg(method)
.dist <- function(i, scale, shape)
list(scale = scale[i], shape = shape[i])
vec <- FALSE
if (length(x) == 1)
distPars <- list(scale = scale, shape = shape)
else {
if (length(shape) == 1 & length(scale) == 1)
distPars <- list(scale = scale, shape = shape)
else if (length(shape) == 1 & length(scale) > 1) {
shape <- rep(shape, length(scale))
distPars <- lapply(1:length(shape), .dist,
scale = scale, shape = shape)
vec <- TRUE
} else if (length(shape) > 1 & length(scale) == 1) {
scale <- rep(scale, length(shape))
distPars <- lapply(1:length(shape), .dist,
scale = scale, shape = shape)
vec <- TRUE
} else {
distPars <- lapply(1:length(shape), .dist,
scale = scale, shape = shape)
vec <- TRUE
}
}
dist <- "weibull"
switch(method,
"conv_direct" = {
if (vec)
return(dCount_allProbs_vec_bi(x, distPars, dist,
conv_steps, time, conv_extrap,
log)
)
else
return(dCount_allProbs_bi(x, distPars, dist,
conv_steps, time, conv_extrap, log)
)
},
"conv_naive" = {
if (vec)
return(dCount_naive_vec_bi(x, distPars, dist,
conv_steps, time, conv_extrap,
log))
else
return(dCount_naive_bi(x, distPars, dist,
conv_steps, time, conv_extrap, log))
},
"conv_dePril" = {
if (vec)
return(dCount_dePril_vec_bi(x, distPars, dist,
conv_steps, time, conv_extrap,
log))
else
return(dCount_dePril_bi(x, distPars, dist,
conv_steps, time, conv_extrap,
log))
},
"series_mat" = {
if (vec)
return(dWeibullCount_mat_vec(x, shape, scale, time, log,
series_terms))
else
return(dWeibullCount_mat(x, shape, scale, time, log,
series_terms))
},
"series_acc" = {
if (vec)
return(dWeibullCount_acc_vec(x, shape, scale, time, log,
series_terms, series_acc_niter,
series_acc_eps)
)
else
return(dWeibullCount_acc(x, shape, scale, time, log,
series_terms, series_acc_niter,
series_acc_eps)
)
}
)
}
#' % Log-likelihood of the the Weibull count process
#' % Compute the log-likelihood of the univariate Weibull count process. Several
#' % methods are provided.
#'
#' \code{dWeibullCount_loglik} computes the log-likelihood.
#'
#' @param na.rm logical, if TRUE \code{NA}'s (produced by taking the log of very
#' small probabilities) will be replaced by the smallest allowed probaility;
#' default is \code{TRUE}.
#' @param weights numeric, vector of weights to apply. If \code{NULL}, a vector
#' of one's will be applied.
#'
#' @return for \code{dWeibullCount_loglik},
#' a double, the log-likelihood of the count process.
#' @rdname dWeibullCount
#' @export
dWeibullCount_loglik <- function(x, shape, scale,
method = c("series_acc", "series_mat",
"conv_direct", "conv_naive", "conv_dePril"),
time = 1, na.rm = TRUE, conv_steps = 100,
conv_extrap = TRUE, series_terms = 50,
series_acc_niter = 300,
series_acc_eps = 1e-10, weights = NULL) {
if (is.null(weights))
weights <- rep(1, length(x))
method <- match.arg(method)
pbs <- dWeibullCount(x, shape, scale, method,
time = time, log = TRUE, conv_steps = conv_steps,
conv_extrap = conv_extrap, series_terms = series_terms,
series_acc_niter = series_acc_niter,
series_acc_eps = series_acc_eps)
pbs <- pbs * weights
if (na.rm)
pbs[is.na(pbs)] <- .logNaReplace()
sum(pbs)
}
#' % Weibull count expected value and variance
#' % Expected value and variance of the weibull count process
#'
#' \code{evWeibullCount} computes the expected value and variance.
#'
#' @param xmax unsigned integer, maximum count to be used.
#' @return for \code{evWeibullCount}, a list with components:
#' \item{ExpectedValue}{expected value,}
#' \item{Variance}{variance.}
#' @rdname dWeibullCount
#' @export
evWeibullCount <- function(xmax, shape, scale,
method = c("series_acc", "series_mat",
"conv_direct", "conv_naive", "conv_dePril"),
time = 1, conv_steps = 100,
conv_extrap = TRUE, series_terms = 50,
series_acc_niter = 300,
series_acc_eps = 1e-10) {
method <- match.arg(method)
x <- 1:xmax
px <- dWeibullCount(x, shape, scale, method, time, FALSE,
conv_steps, conv_extrap, series_terms,
series_acc_niter, series_acc_eps)
## adjust probability
px[px < 0] <- 0
px <- zapsmall(px)
##px <- px / sum(px)
ev <- sum(x * px)
ev2 <- sum(x^2 * px)
var <- ev2 - ev^2
if (var < 0) {
var <- NA
warning("failed to compute variance! NA returned")
}
list(ExpectedValue = ev, Variance = var)
}