-
Notifications
You must be signed in to change notification settings - Fork 0
/
bernoulli_funs.R
99 lines (90 loc) · 3.32 KB
/
bernoulli_funs.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
logist <- function(t){1 / (1 + exp(-t))}
f_cond <- function(Y, w, mu, lam){
# Y and w should be be n x r and n x 1, respectively
# mu and lam should be scalars
r <- ncol(Y)
y <- as.vector(Y)
P <- matrix(dbinom(y,
size = 1,
prob = logist(mu + lam * rep(w, r))),
ncol = r)
# Return
apply(P, 1, prod)
}
f_marg <- function(Y, mu, lam, num_nodes, dist){
# Y should be n x r matrix; mu and lam scalar
quad <- statmod::gauss.quad.prob(n = num_nodes, dist = dist)
quad$nodes <- quad$nodes - ifelse(dist == "gamma", 1, 0)
# Return
apply(Y, 1, function(y){
sum(quad$weights * f_cond(Y = matrix(rep(y, each = num_nodes),
nrow = num_nodes),
w = quad$nodes,
mu = mu,
lam = lam))
})
}
score <- function(Y, mu, lam, num_nodes = 10, dist, modified = FALSE,
tol = 1e-12){
# Y should be n x r matrix
# mu should be scalar
n <- nrow(Y); r <- ncol(Y)
quad <- statmod::gauss.quad.prob(n = num_nodes, dist = dist)
quad$nodes <- quad$nodes - ifelse(dist == "gamma", 1, 0)
score_mu <- apply(Y, 1, function(y){
-sum(quad$weights * (r * logist(mu + quad$nodes * lam) - sum(y)) *
f_cond(Y = matrix(rep(y, each = num_nodes),
nrow = num_nodes),
w = quad$nodes,
mu = mu,
lam = lam))
})
if((abs(lam) > tol) | (!modified)){
score_lam <- apply(Y, 1, function(y){
-sum(quad$weights * quad$nodes *
(r * logist(mu + quad$nodes * lam) - sum(y)) *
f_cond(Y = matrix(rep(y, each = num_nodes), nrow = num_nodes),
w = quad$nodes, mu = mu, lam = lam))
})
} else{
score_lam <- apply(Y, 1, function(y){
prob <- logist(mu + quad$nodes * lam)
-sum(quad$weights * quad$nodes^2 *
((r * prob - sum(y))^2 - r * prob * (1 - prob)) *
f_cond(Y = matrix(rep(y, each = num_nodes),
nrow = num_nodes), w = quad$nodes, mu = mu,
lam = lam))
})
}
# Return
cbind(score_mu, score_lam) * (1 / f_marg(Y, mu, lam, num_nodes, dist = dist))
}
finf <- function(mu, lam, r, num_nodes, dist, modified = FALSE, tol = 1e-12)
{
Y <- matrix(rep(0, r), nrow = 1)
I <- tcrossprod(c(score(Y = Y, mu = mu, lam = lam,
num_nodes = num_nodes, dist = dist,
modified = modified, tol = tol))) *
f_marg(Y = Y, mu = mu, lam = lam, num_nodes = 10, dist = dist)
for(ii in 1:r){
Y[1, ii] <- 1
I <- I + choose(r, ii) * tcrossprod(c(score(Y = Y, mu = mu, lam = lam,
num_nodes = num_nodes,
dist = dist,
modified = modified,
tol = tol))) *
f_marg(Y = Y, mu = mu, lam = lam, num_nodes = 10, dist = dist)
}
#Return
I
}
gen_mvber <- function(n, mu, lam, r, dist){
if(dist == "normal"){
w <- rnorm(n)
} else{
w <- rexp(n) - 1
}
lin_pred <- matrix(lam * rep(w, r), ncol = r) + mu
# Return
matrix(rbinom(n = r * n, size = 1, prob = logist(as.vector(lin_pred))), ncol = r)
}