-
Notifications
You must be signed in to change notification settings - Fork 0
/
Zero2.R
65 lines (37 loc) · 1.26 KB
/
Zero2.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
ll <- function(theta,y){
## evaluates the likelihood function of a
## beta-binomial distribution given theta and y
## note : n is specified inside the function
lastcol <- ncol(y)
y1 <- y[,lastcol]
X <- y[,1:ncol(y)-1]
ones <- rep(1,nrow(y))
X <- cbind(ones,X)
n <- 10
rho <- 1/(1+exp(-theta[length(theta)]))
print(rho)
phi <- (1-rho)/rho
print(phi)
B <- theta[c(1:length(theta)-1)]
Gamma <- B[c(1,2,3)]
Alpha <- B[c(4,5,6)]
Beta <- Gamma + exp(Alpha)
loglik <- vector(mode="numeric", length(y1))
## YOU HAVE TO LOOP OVER THE i VALUES HERE
for (i in 1:length(y1)){
pi <- 1/(1+exp(-X[i,]%*%Beta))
mupi <- 1/(1+exp(-X[i,]%*%Gamma))
mu <- mupi/pi
if(y1[i] != 0){
loglik[i] <- log(pi) + lgamma(phi) - lgamma(phi*mu) -
lgamma(phi*(1-mu)) + lgamma(y1[i]+phi*mu) +
lgamma(n - y1[i] + phi*(1-mu)) - lgamma(n+phi) -
lgamma(y1[i]+1)-lgamma(n-y1[i]+1) + lgamma(n+1)
}else{
loglik[i] <- log(1-pi + choose(n,y1[i])*pi*exp(lbeta(y1[i]+phi*mu,n -y1[i] + phi*(1-mu))-lbeta((phi*mu),phi*(1-mu))))
}
}
final <- sum(loglik)
final1 <- -final
print(final1)
}