-
Notifications
You must be signed in to change notification settings - Fork 1
/
chisq.R
82 lines (73 loc) · 2.38 KB
/
chisq.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
## --##--##--##--##--##--##--##--##--##--##--##--##--##--##
## Functions related to the Chi Square distribution ##
## --##--##--##--##--##--##--##--##--##--##--##--##--##--##
#' @param df degrees of freedom for "parent" distribution
#' @rdname rtrunc
#' @export
rtruncchisq <- rtrunc.chisq <- function(n, df, a = 0, b = Inf, faster = FALSE) {
class(n) <- "trunc_chisq"
if (faster) {
family <- gsub("trunc_", "", class(n))
parms <- mget(ls())[grep("^faster$|^n$|^family$", ls(), invert = TRUE)]
return(rtrunc_direct(n, family, parms, a, b))
} else {
parms <- mget(ls())[grep("^faster$", ls(), invert = TRUE)]
return(sampleFromTruncated(parms))
}
}
#' @export
dtrunc.trunc_chisq <- function(y, df, eta, a = 0, b = Inf, ...) {
if (missing(eta)) {
eta <- parameters2natural.parms_chisq(c("df" = df))
}
df <- natural2parameters.parms_chisq(eta)
dens <- rescaledDensities(y, a, b, dchisq, pchisq, df)
return(dens)
}
#' @rdname dtrunc
#' @export
dtruncchisq <- dtrunc.trunc_chisq
#' @export
empiricalParameters.trunc_chisq <- function(y, ...) {
# Returns empirical parameter estimate for df
parms <- c("df" = mean(y))
class(parms) <- "parms_chisq"
return(parms)
}
#' @method sufficientT trunc_chisq
sufficientT.trunc_chisq <- function(y) {
return(suff.T = log(y))
}
#' @export
natural2parameters.parms_chisq <- function(eta, ...) {
# eta: The natural parameters in a Chi Square distribution
# returns df
if (length(eta) != 1) stop("Eta must be one single number")
df <- c(df = 2 * (eta[[1]] + 1))
class(df) <- class(eta)
return(df)
}
#' @export
parameters2natural.parms_chisq <- function(parms, ...) {
# parms: The parameter lambda in a Chi Square distribution
# returns the natural parameters
eta <- prepEta(parms / 2 - 1, class(parms))
return(eta)
}
#' @method getGradETinv parms_chisq
getGradETinv.parms_chisq <- function(eta, ...) {
# eta: Natural parameter
# return the inverse of E.T differentiated with respect to eta
return(A = 1 / sum(1 / (as.vector(eta) + (1:1e6))^2))
}
#' @method getYseq trunc_chisq
getYseq.trunc_chisq <- function(y, y.min = 0, y.max, n = 100) {
mean <- mean(y, na.rm = TRUE)
var.y <- var(y, na.rm = TRUE)
lo <- max(round(y.min), 0)
hi <- min(y.max, round(mean + 10 * sqrt(var.y)))
out <- seq(lo, hi, length = n)
out <- out[out > 0] # prevents NaN as sufficient statistics
class(out) <- class(y)
return(out)
}