-
Notifications
You must be signed in to change notification settings - Fork 1
/
stats_bernoulli.R
120 lines (109 loc) · 3.75 KB
/
stats_bernoulli.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
#' @title Bernoulli Trials
#' @description Conduct bernoulli trials
#' @param x predictor variable(s), Default: NULL
#' @param x.names optional names for predictor variable(s), Default: NULL
#' @param DF data for analysis
#' @param params define parameters to observe, Default: NULL
#' @param initial.list initial values for analysis, Default: list()
#' @param ... further arguments passed to or from other methods
#' @examples
#' ## Create coin toss data: heads = 50 and tails = 50
#' #fair.coin<- as.matrix(c(rep("Heads",50),rep("Tails",50)))
#' #colnames(fair.coin) <- "X"
#' #fair.coin <- bfw(project.data = fair.coin,
#' # x = "X",
#' # saved.steps = 50000,
#' # jags.model = "bernoulli",
#' # jags.seed = 100,
#' # ROPE = c(0.4,0.6),
#' # silent = TRUE)
#' #fair.coin.freq <- binom.test( 50000 * 0.5, 50000)
#'
#' ## Create coin toss data: heads = 20 and tails = 80
#' #biased.coin <- as.matrix(c(rep("Heads",20),rep("Tails",80)))
#' #colnames(biased.coin) <- "X"
#' #biased.coin <- bfw(project.data = biased.coin,
#' # x = "X",
#' # saved.steps = 50000,
#' # jags.model = "bernoulli",
#' # jags.seed = 101,
#' # initial.list = list(theta = 0.7),
#' # ROPE = c(0.4,0.6),
#' # silent = TRUE)
#' #biased.coin.freq <- binom.test( 50000 * 0.8, 50000)
#'
#' ## Print Bayesian and frequentist results of fair coin
#' #fair.coin$summary.MCMC[,c(3:6,9:12)]
#'
#' ## Mode ESS HDIlo HDIhi ROPElo ROPEhi ROPEin n
#' ## 0.505 50480.000 0.405 0.597 2.070 2.044 95.886 100.00
#'
#' #sprintf("Frequentist: %.3f [%.3f , %.3f], p = %.3f" ,
#' # fair.coin.freq$estimate ,
#' # fair.coin.freq$conf.int[1] ,
#' # fair.coin.freq$conf.int[2] ,
#' # fair.coin.freq$p.value)
#'
#' ## [1] "Frequentist: 0.500 [0.496 , 0.504], p = 1.000"
#'
#' ## Print Bayesian and frequentist results of biased coin
#' #biased.coin$summary.MCMC[,c(3:6,9:12)]
#'
#' ## Mode ESS HDIlo HDIhi ROPElo ROPEhi ROPEin n
#' ## 0.803 50000.000 0.715 0.870 0.000 99.996 0.004 100.000
#'
#' #sprintf("Frequentist: %.3f [%.3f , %.3f], p = %.3f" ,
#' # biased.coin.freq$estimate ,
#' # biased.coin.freq$conf.int[1] ,
#' # biased.coin.freq$conf.int[2] ,
#' # biased.coin.freq$p.value)
#'
#' ## [1] "Frequentist: 0.800 [0.796 , 0.803], p = 0.000"
#' @seealso
#' \code{\link[stats]{complete.cases}}
#' @rdname StatsBernoulli
#' @export
#' @importFrom stats complete.cases
StatsBernoulli <- function(x = NULL,
x.names = NULL,
DF,
params = NULL,
initial.list = list(),
...
) {
# Create crosstable for x parameters
n.data <- as.data.frame(table(DF[, x]))
# Fetch x parameter and exclude missing observations
x <- DF[stats::complete.cases(DF[, x]), x]
# Final name list
job.names <- paste(levels(x)[[1]],levels(x)[[2]], sep=" vs. ")
# Convert x to numeric vector (start at zero)
x <- as.numeric(as.factor(x)) - 1
# Number of observations
n <- length(x)
# Paramter(s) of interest
if (is.null(params)) params <- c("theta")
# Create data for Jags
data.list <- list(
x = x,
n = n
)
# Priors
if (!length(initial.list)) {
initial.list <- list(
theta = 0.5
)
}
# Create name list
name.list <- list(
job.names = job.names
)
# Return data list
return ( list (
params = params,
data.list = data.list,
name.list = name.list,
initial.list = initial.list,
n.data = n.data
))
}