Nothing
X_update_d_alphaco_i_new2 <- function(){
##
sub_freq_locus <- BBB2$freq_locus
sub_freq_pop <- BBB2$freq_pop
old_alpha <- BBB2$d_alpha
sub_sample_size <- BBB2$sample_size
sub_GROUP <- BBB2$GROUP
locnum2 <- length(BBB2$d_alpha)
##
old_theta <- outer(old_alpha,BBB2$d_beta,"+")
old_theta <- exp(-(old_theta))
old1 <- lgamma(old_theta)-lgamma(sub_sample_size + old_theta)
val <- matrix(0,locnum2,dim(BBB2$freq_locus)[2])
for(xx in 1:BBB2$popnum){
xyz <- old_theta[,xx]*sub_freq_locus
val <- val + lgamma(sub_freq_pop[[xx]] + xyz) - lgamma(xyz)
}
old1 <- rowSums(old1)
oldL <- old1 + rowSums(val,na.rm=TRUE)
#oldL <- tapply(oldL,BBB2$GROUP,sum)
e <- rnorm(locnum2,0,BBB2$var_prop_alpha)
#e <- rnorm(locnum2,0,1)
new_alpha <- old_alpha + e
new_theta <- outer(new_alpha,BBB2$d_beta,"+")
new_theta <- exp(-(new_theta))
new1 <- lgamma(new_theta)-lgamma(sub_sample_size+new_theta)
val <- matrix(0,locnum2,dim(BBB2$freq_locus)[2])
for(xx in 1:BBB2$popnum){
xyz <- new_theta[,xx]*sub_freq_locus
val <- val + lgamma(sub_freq_pop[[xx]] + xyz) - lgamma(xyz)
}
new1 <- rowSums(new1)
newL <- new1 + rowSums(val,na.rm=TRUE)
#newL <- tapply(newL,BBB2$GROUP,sum)
#newL <- newL[BBB2$GROUP]
#oldL <- oldL[BBB2$GROUP]
A <- -oldL + newL + (old_alpha^2-new_alpha^2)/(2*BBB2$sd_prior_alpha^2)
r <- runif(locnum2)
check <- log(r) < A
BBB2$d_alpha[check] <- new_alpha[check]
BBB2$acc_alpha[check] <- BBB2$acc_alpha[check] + 1
}
# END of Function
X_jump_model_codominant2 <- function(){
old_alpha <- BBB2$d_alpha
excluded <- which(!BBB2$alpha_included)
BBB2$d_alpha[BBB2$alpha_included] <- 0
if(length(excluded)!=0){
# choose random SNP of each group excluded
#groups <- BBB2$GROUP[excluded]
#ids <- tapply(excluded,groups,function(x){rand <- sample(1:length(x),1); return(x[rand]) })
#e_group <- rnorm(length(ids),mean=BBB2$mean_alpha[ids],sd=sqrt(BBB2$var_alpha[ids]))
#BBB2$d_alpha[ids] <- e_group
e_group <- rnorm(length(excluded),mean=BBB2$mean_alpha[excluded],sd=sqrt(BBB2$var_alpha[excluded]))
BBB2$d_alpha[excluded] <- e_group
}
BBB2$alpha_included <- !BBB2$alpha_included
## old_theta
old_theta <- outer(old_alpha,BBB2$d_beta,"+")
old_theta <- exp(-(old_theta))
old1 <- lgamma(old_theta)-lgamma(BBB2$sample_size+old_theta)
val <- BBB2$GLOBAL_INIT1
for(xx in 1:BBB2$popnum){
xyz <- old_theta[,xx]*BBB2$freq_locus
val <- val + lgamma(BBB2$freq_pop[[xx]]+xyz) - lgamma(xyz)
}
old1 <- rowSums(old1)
oldL <- old1 + rowSums(val,na.rm=TRUE)
oldL <- tapply(oldL,BBB2$GROUP,sum)
## new_theta
new_theta <- outer(BBB2$d_alpha,BBB2$d_beta,"+")
new_theta <- exp(-(new_theta))
new1 <- lgamma(new_theta)-lgamma(BBB2$sample_size+new_theta)
val <- BBB2$GLOBAL_INIT1
for(xx in 1:BBB2$popnum){
xyz <- new_theta[,xx]*BBB2$freq_locus
val <- val + lgamma(BBB2$freq_pop[[xx]]+xyz) - lgamma(xyz)
}
new1 <- rowSums(new1)
newL <- new1 + rowSums(val,na.rm=TRUE)
newL <- tapply(newL,BBB2$GROUP,sum)
A_group <- -oldL + newL
A <- A_group[BBB2$GROUP]
included <- BBB2$alpha_included #==TRUE
excluded <- !included
A_included <- 0
A_excluded <- 0
if(any(included)){
A_included <- -log(BBB2$tempting)
}
if(any(excluded)){
A_excluded <- log(BBB2$tempting)
}
A[included] <- A[included] + A_included
A[excluded] <- A[excluded] + A_excluded
A <- tapply(A,BBB2$GROUP,unique)
write.table(A,"Jump_A")
r <- runif(BBB2$N.REGIONS)
check <- log(r) > A
group_id <- which(check)
if(length(group_id)>0){
.Call("update_group_Jump2",
BBB2$d_alpha,
old_alpha,
as.integer(BBB2$GROUP),
as.integer(group_id),
BBB2$alpha_included
)
}
}# End of Function
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.