-
Notifications
You must be signed in to change notification settings - Fork 0
/
Pinitialize.R
154 lines (100 loc) · 4.19 KB
/
Pinitialize.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
Pinitialize = function(Y,D,N,K,time,censoring){
Time <- cbind(time, censoring)
## HYPER PRIORS
## Hyper parameters of the DP
shape.alpha <- 2
rate.alpha <- 1
## Hyperparameters for the GMM
beta = D+1
ro = 0.5
source('rchinese.R')
alpha = rgamma(1, shape = shape.alpha, rate = rate.alpha )
c <- rchinese(N,alpha)
f <- table(factor(c, levels = 1:max(c)))
## Empirical Bayes Estimate of the Hyperparameters
epsilon = as.vector(apply(Y,2,mean))
W = diag(diag(cov(Y)))
## Initialization of the parameters for Gaussian Mixture
mu = matrix(data = NA, nrow = K, ncol = D)
S = array(data = NA, dim =c(K,D,D))
#Sparsity controlling hyperparameter of the BAYESIAN LASSO MODEL
r =1
si = 1.78
##Actual parameters
lambda2 <- numeric(K)
tau2 = matrix(data = NA, nrow = K, ncol = D)
betahat = matrix(data = NA, nrow = K, ncol = D)
sigma2 <- rep(NA, K)
beta0 <- rep(NA, K)
That <- numeric(N)
## Fitting a linear model to the whole model
Ysc <- scale(Y[1:N,1:D], center = TRUE, scale =TRUE)
lm.data <- lm(time ~ Ysc)
sig2.dat <- var(lm.data$residuals)
## Set Some Initial Values for the Cluster Parameters
source('priorPARAMETERS.R')
disclass <- table(factor(c, levels = 1:K))
activeclass <- which(disclass!=0)
for ( j in 1:length(activeclass)){
priorone <- priordraw(beta, W, epsilon, ro, r, si, N, D, sig2.dat)
mu[activeclass[j],] <- (priorone$mu)
S[activeclass[j],1:D,1:D] <- priorone$Sigma
beta0[activeclass[j]] <- priorone$beta0
sigma2[activeclass[j]] <- priorone$sigma2
betahat[activeclass[j],1:D] <- priorone$betahat
lambda2[activeclass[j]] <- priorone$lambda2
tau2[activeclass[j], 1:D] <- priorone$tau2
}
# The Time has to be initialized
source('posteriorCensoredTime.R')
ti <- updatetime(c, Y, Time,That, beta0, betahat, sigma2)
That <- ti$time
################# K-Means BLASSO INITIALIZATION ############################################
G <- F
k.data <- kmeans(Y,F,nstart =10)
c <- k.data$cluster
c.kmeans <- c
#### Under special cases
###### c <- c.true
prior.numclust <- table(factor(c, levels = 1:K))
prior.activeclass<- which(prior.numclust!=0)
### The means are set using the k-means
for ( i in 1:length(prior.activeclass)){
mu[prior.activeclass[i],1:D] <- k.data$centers[i,1:D]
S[prior.activeclass[i],1:D,1:D] <- priordraw(beta, W, epsilon, ro, r, si,N,D, sig2.dat)$Sigma
lclust <- which(c == prior.activeclass[i])
reg.blas <- 0
sum <- c(0)
coeff <- 0
Ytemp <- matrix(NA, nrow = length(lclust), ncol = D)
Ytemp <- scale(Y[lclust,1:D], center = TRUE, scale = TRUE)
### Part where I use the MONOMVN PACKAGE
Ttemp <- as.vector(That[lclust])
ntemp <- length(lclust)
reg.blas <- blasso(Ytemp, Ttemp, T = 300,thin = 50, RJ = TRUE, mprior = 0.0 ,normalize = TRUE, verb = 0)
sum <- summary(reg.blas, burnin= 100)
## Selecting those features which are relevant
coeff <- unlist(lapply(strsplit(sum$coef[3,], split = ":"), function(x) as.numeric(unlist(x)[2])))
beta0[prior.activeclass[1]] <- coeff[1]
indexplusone <- D+1
ind <- 2:indexplusone
betahat[prior.activeclass[i], ] <- coeff[ind]
ta <- unlist(lapply(strsplit(sum$tau2i[3,], split = ":"), function(x) as.numeric(unlist(x)[2])))
tau2[prior.activeclass[i],] <- ta
sigma2[prior.activeclass[i]] <- sum$s2[3]
lambda2[prior.activeclass[i]] <- sum$lambda2[3]
}
## Deleting those values which are no longer relevant
g <- table(factor(c, levels = 1:K))
inactive <- which(g==0)
for ( i in 1:length(inactive)){
mu[inactive[i],1:D] <- NA
S[inactive[i],1:D,1:D] <- NA
beta0[inactive[i]] <- NA
sigma2[inactive[i]] <- NA
betahat[inactive[i],1:D] <- NA
lambda2[inactive[i]] <- NA
tau2[inactive[i], 1:D] <- NA
}
return(list(Time = 'Time', 'alpha' = alpha,'That' = That, 'c'= c,'c.kmeans' = c.kmeans, 'epsilon'= epsilon,'W' = W,'mu'= mu,'S'= S,'beta0'= beta0,'betahat'= betahat,'sigma2'= sigma2, 'lambda2'= lambda2,'tau2'= tau2,'sig2.dat'=sig2.dat))
}