/
master.R
85 lines (72 loc) · 4.12 KB
/
master.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
####################blackrhino CD4 #############
#setwd("~/Lynn/cts_measurements_v1.2")
set.seed(984)
source("repmat.R")
dyn.load("model.so")
source("realdata_rv144.R")
T = as.integer(40000);
source("init_realdata_rv144.R")
ttt=as.integer(2000);
n_ss = as.integer(n_s); n_uu = as.integer(n_u)
ybars = apply(ybar_s, 3, function(x) as.vector(x))
ybaru = apply(ybar_u, 3, function(x) as.vector(x))
ys2s = apply(ys2_s, 3, function(x) as.vector(x))
ys2u = apply(ys2_u, 3, function(x) as.vector(x))
gammat = apply(gamma,3, function(x) as.integer(x))
##########################case-control CD4 ##############
set.seed(536)
source("repmat.R")
dyn.load("model.so")
source("realdata_Casecontrol.R")
T = as.integer(4000);
source("init_realdata_Casecontrol.R")
ttt=as.integer(2000);
n_ss = as.integer(n_s); n_uu = as.integer(n_u)
ybars = apply(ybar_s, 3, function(x) as.vector(x))
ybaru = apply(ybar_u, 3, function(x) as.vector(x))
ys2s = apply(ys2_s, 3, function(x) as.vector(x))
ys2u = apply(ys2_u, 3, function(x) as.vector(x))
rm(ybar_s,ybar_u,ys2_u,ys2_s)
gammat = apply(gamma,3, function(x) as.integer(x))
################## parameter tuning ######################
ptm <- proc.time()
result<-.Call("model",T=T,I=I, K=K, M=M, ttt=ttt, SS=as.integer(1),alpha_u=alpha_u, alpha_s=alpha_s, mu_u=mu_u, mu_s=mu_s, alpha=alpha,
beta=beta, gamma=gammat,n_s=n_ss, n_u=n_uu,varp_u=varp_u, lambda_u=lambda_u, indi=indi, d=d, ybar_s = ybars, ybar_u=ybaru,
ys2_s = ys2s, ys2_u = ys2u, a=a, b=b,lambda=lambda,mk = mk, Istar = Istar, mKstar = mKstar, pp = pp, pb1 = pb1,
pb2 = pb2,lambda_s = lambda_s, var_1 = varp_s1,var_2 = varp_s2,p_var = pvar_s, p_vars = pvar_mus,var_1s = sqrt_mus1, var_2s=sqrt_mus2,m_s=m_s,Sigma_s =Sigma_s,
p_varu=pvar_muu,var_1u=sqrt_muu1,var_2u=sqrt_muu2,m_u=m_u, Sigma_u=Sigma_u,p_vara=pvar_alpha, var_1a=sqrt_alpha1,var_2a=sqrt_alpha2,sig_alpha1=sig_alpha1,alpha1=alpha1,
p_varb=pvar_beta,var_1b=sqrt_beta1,var_2b=sqrt_beta2,sig_beta1 = sig_beta1, beta1= beta1, A_alphau=A_alphau,A_alphas=A_alphas, A_gm=A_gm, A_mus=A_mus, A_muu=A_muu,
A_alpha=A_alpha,A_beta=A_beta, Tune = as.integer(1), pgamma = 0.4)
proc.time() - ptm
Istar = result$Istar
mk = result$mk
mKstar = result$mKstar
var_1a = result$var_1a
var_1b = result$var_1b
var_2a = result$var_2a
var_2b = result$var_2b
alpha_u[1,] = alpha_u[T,]; alpha_s[1,] = alpha_s[T,]; gammat[,1] = gammat[,T];
mu_s[1,] = mu_s[T,]; mu_u[1,]=mu_u[T,]; alpha[1] = alpha[T]; beta[1] = beta[T];
sTT =25;
#######################################################
for (stt in 1:sTT) {
result<-.Call("model",T=T,I=I, K=K, M=M, ttt=ttt, SS=as.integer(1),alpha_u=alpha_u, alpha_s=alpha_s, mu_u=mu_u, mu_s=mu_s, alpha=alpha,
beta=beta, gamma=gammat,n_s=n_ss, n_u=n_uu,varp_u=varp_u, lambda_u=lambda_u, indi=indi, d=d, ybar_s = ybars, ybar_u=ybaru,
ys2_s = ys2s, ys2_u = ys2u, a=a, b=b,lambda=lambda,mk = mk, Istar = Istar, mKstar = mKstar, pp = pp, pb1 = pb1,
pb2 = pb2,lambda_s = lambda_s, var_1 = varp_s1,var_2 = varp_s2,p_var = pvar_s, p_vars = pvar_mus,var_1s = sqrt_mus1, var_2s=sqrt_mus2,m_s=m_s,Sigma_s =Sigma_s,
p_varu=pvar_muu,var_1u=sqrt_muu1,var_2u=sqrt_muu2,m_u=m_u, Sigma_u=Sigma_u,p_vara=pvar_alpha, var_1a=sqrt_alpha1,var_2a=sqrt_alpha2,sig_alpha1=sig_alpha1,alpha1=alpha1,
p_varb=pvar_beta,var_1b=sqrt_beta1,var_2b=sqrt_beta2,sig_beta1 = sig_beta1, beta1= beta1, A_alphau=A_alphau,A_alphas=A_alphas, A_gm=A_gm, A_mus=A_mus, A_muu=A_muu,
A_alpha=A_alpha,A_beta=A_beta, Tune = as.integer(0),pgamma = 0.5)
if (stt == sTT) {break}
Istar = result$Istar
mk = result$mk
mKstar = result$mKstar
alpha_u[1,] = alpha_u[T,]; alpha_s[1,] = alpha_s[T,]; gammat[,1] = gammat[,T];
mu_s[1,] = mu_s[T,]; mu_u[1,]=mu_u[T,]; alpha[1] = alpha[T]; beta[1] = beta[T];
cat("iteration = ", iter <- stt, "\n")
}
##########################################################
for ( tt in 1:T) {
gamma[,,tt]=matrix(gammat[,tt], nrow=I)
}
Mgamma = apply(gamma,c(1,2),sum)/T