R/update_slim.R
In pievos101/BlockFeST: Bayesian Calculation of Region-Specific Fixation Index to Detect Local Adaptation
Defines functions
X_update_d_alphaco_i_new_id
X_update_d_alphaco_i_new
X_jump_model_codominant
update_freq_codominant
update_d_betaco
X_update_d_alphaco_i
#source("math.R")
#constants
#e_f <- 0.05
#sd_prior_beta <- 1 #3
#mean_prior_beta <- -1 #-5.15
#lambda <- 1
X_update_d_alphaco_i <- function(id){ # id
if(id[1]==FALSE){ # Pilot runs
id <- 1:BBB2$locnum
locnum2 <- BBB2$locnum
}else{ # Main Loop
locnum2 <- length(id)
}
##
sub_freq_locus <- BBB2$freq_locus[id,]
sub_freq_pop <- lapply(BBB2$freq_pop,function(yy){return(yy[id,])})
old_alpha <- BBB2$d_alpha[id]
sub_sample_size <- BBB2$sample_size[id,]
sub_GROUP <- BBB2$GROUP[id]
##
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)
## modification ------------------------------------
oldL <- tapply(oldL,sub_GROUP,sum)
##--------------------------------------------------
# e <- rnorm(locnum2,0,sqrt(BBB2$var_prop_alpha))
# new_alpha <- old_alpha + e
### modification -----------------------------------------------
if(!BBB2$PILOT){
whatgroup <- unique(sub_GROUP)
xxx <- BASIX.match(whatgroup,BBB2$GROUP)
e_group <- rnorm(length(whatgroup),0,sqrt(BBB2$var_prop_alpha[xxx]))
#
BBB2$new_alpha <- old_alpha
BBB2$jjj <- 1
sapply(whatgroup,function(x){
what <- x==sub_GROUP
BBB2$new_alpha[what] <- old_alpha[what] + e_group[BBB2$jjj]
BBB2$jjj <- BBB2$jjj + 1
})
new_alpha <- BBB2$new_alpha
#
}else{
e_group <- rnorm(BBB2$N.REGIONS,0,sqrt(BBB2$var_prop_alpha))
new_alpha <- old_alpha + e_group[BBB2$GROUP]
whatgroup <- 1:BBB2$N.REGIONS
}
#
##--------------------------------------------------------------
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)
## modification ----------------------------
newL <- tapply(newL,sub_GROUP,sum)
old_alpha_group <- tapply(old_alpha,sub_GROUP,function(x){return(x[1])})
new_alpha_group <- tapply(new_alpha,sub_GROUP,function(x){return(x[1])})
# print(old_alpha)
# print(new_alpha)
##-----------------------------------------
#A <- -oldL + newL + (old_alpha^2-new_alpha^2)/(2*BBB2$sd_prior_alpha^2)
## mod ------------------------
A <- -oldL + newL + (old_alpha_group^2-new_alpha_group^2)/(2*BBB2$sd_prior_alpha^2)
##-----------------------------
# r <- runif(locnum2)
## mod ---------------------------
r <- runif(length(whatgroup))
## mod ---------------------------
check <- log(r) < A
# BBB2$d_alpha[id[check]] <- new_alpha[check]
# BBB2$acc_alpha[id[check]] <- BBB2$acc_alpha[id[check]] + 1
# modification -------------------------------------
# update groups
group_idX <- which(check)
group_id <- whatgroup[group_idX]
# make this also in C
if(length(group_idX)>0){
#res <- .Call("update_group2",group_id,BBB2$GROUP,BBB2$d_alpha,new_alpha,BBB2$acc_alpha)
#BBB2$d_alpha <- res[[1]]
#BBB2$acc_alpha <- res[[2]]
mm <- sapply(group_id,function(xx){xx==BBB2$GROUP})
BBB2$iii <- 1
apply(mm,2,function(x){
BBB2$d_alpha[x] <- new_alpha_group[group_idX[BBB2$iii]] # new_alpha[x]
BBB2$acc_alpha[x] <- BBB2$acc_alpha[x] + 1
BBB2$iii <- BBB2$iii + 1
})
}
# ---------------------------------------------------
}# End of Function
update_d_betaco <- function(){
# Init
#val <- rep(NA,2)
val <- BBB2$GLOBAL_INIT2
# ---------------
old_beta <- BBB2$d_beta
old_theta <- outer(BBB2$d_alpha,old_beta,"+")
old_theta <- exp(-(old_theta))
#cat("old_theta", old_theta,"\n")
old1 <- lgamma(old_theta)-lgamma(BBB2$sample_size+old_theta)
#cat("old1: ",old1,"\n")
for(xx in 1:BBB2$popnum){
xyz <- old_theta[,xx]*BBB2$freq_locus
val[xx] <- sum(lgamma(BBB2$freq_pop[[xx]]+xyz)-lgamma(xyz),na.rm=TRUE)
}
#cat("val:",val,"\n")
old1 <- colSums(old1)
oldL <- old1 + val
#--------- calculate new likelihood
e <- rnorm(BBB2$popnum,0,BBB2$var_prop_beta)
BBB2$d_beta <- BBB2$d_beta + e
new_theta <- outer(BBB2$d_alpha,BBB2$d_beta,"+")
new_theta <- exp(-(new_theta))
new1 <- lgamma(new_theta)-lgamma(BBB2$sample_size+new_theta)
for(xx in 1:BBB2$popnum){
xyz <- new_theta[,xx]*BBB2$freq_locus
val[xx] <- sum(lgamma(BBB2$freq_pop[[xx]]+xyz)-lgamma(xyz),na.rm=TRUE)
}
new1 <- colSums(new1)
newL <- new1 + val
A <- -oldL + newL
#+ (old_beta -BBB2$d_beta)*(old_beta+BBB2$d_beta-2*BBB2$mean_prior_beta)/(2*BBB2$sd_prior_beta^2)
#+ (old_beta^2-BBB2$d_beta^2)/(2*BBB2$sd_prior_beta^2)
#+ (old_alpha_group^2-new_alpha_group^2)/(2*BBB2$sd_prior_alpha^2)
r <- runif(BBB2$popnum)
check <- log(r)>A
#cat("-oldL: ",-oldL, "\n")
#cat("newL: ",newL, "\n")
#print(A)
#print(old_beta)
#print(check)
#print(r)
#if(any(is.na(check))){
#print("yes")
#check[is.na(check)] <- TRUE
#}
BBB2$d_beta[check] <- old_beta[check]
BBB2$acc_beta[!check] <- BBB2$acc_beta[!check] + 1
}# End of Function
update_freq_codominant <- function(){
lambda <- 1
# choose randomly two haplotype of each locus
# samp <- sapply(BBB2$hapcount,function(x){return(sample(1:x,size=2)) })
theta <- outer(BBB2$d_alpha,BBB2$d_beta,"+")
theta <- exp(-theta)
# ok
#BBB2$iii <- 1
# frequencies of haplotypes of the two random haps
#samp_freq <- apply(samp,2,function(x){
# r <- BBB2$freq_locus[BBB2$iii,x]
# BBB2$iii <- BBB2$iii+1
# return(r)
# })
#ok
## mod
samp_freq <- t(BBB2$freq_locus)
## end mod
#print(dim(samp_freq))
#print(dim(BBB2$freq_locus))
# haplotype distribution for each pop of the two random haps
# BBB2$iii <- 1
samp_pop_freq <- vector("list",BBB2$popnum)
for(yy in 1:BBB2$popnum){
samp_pop_freq[[yy]] <- t(BBB2$freq_pop[[yy]])
#samp_pop <- apply(samp,2,function(x){
# r <- BBB2$freq_pop[[yy]][BBB2$iii,x]
# BBB2$iii <- BBB2$iii+1
# return(r)
# })
#samp_pop_freq[[yy]] <- samp_pop
#BBB2$iii <- 1
}
#ok
#print(dim(samp_pop_freq[[1]]))
old_L <- BBB2$GLOBAL_INIT3
# 1 <- m
# 2 <- n
for(xx in 1:BBB2$popnum){
old_L <- old_L + lgamma(samp_pop_freq[[xx]][1,] + theta[,xx] * samp_freq[1,]) -
lfactorial (samp_pop_freq[[xx]][1,]) - lgamma(theta[,xx]*samp_freq[1,])+
lgamma(samp_pop_freq[[xx]][2,]+theta[,xx]*samp_freq[2,]) -
lfactorial(samp_pop_freq[[xx]][2,]) - lgamma(theta[,xx]*samp_freq[2,])
}
old_m <- samp_freq[1,]
old_n <- samp_freq[2,]
#### BwPrp
old_m_n <- old_m + old_n
old_m_ancestral <- old_m + BBB2$e_ancestral
old_m_ancestral2 <- old_m - BBB2$e_ancestral
checkmin <- cbind(old_m_n,old_m_ancestral)
checkmax <- cbind(BBB2$GLOBAL_INIT3,old_m_ancestral2)
minn <- apply(checkmin,1,min)
maxx <- apply(checkmax,1,max)
BwPrp <- minn - maxx
###################
# change frequencies
u <- runif(BBB2$locnum,0,BwPrp) + maxx
# samp_freq_new
samp_freq_new <- samp_freq
samp_freq_new[1,] <- u
samp_freq_new[2,] <- old_n + old_m - u
### FwPrp
ss1 <- colSums(samp_freq_new)
ss2 <- samp_freq_new[1,] + BBB2$e_ancestral
ss3 <- samp_freq_new[1,] - BBB2$e_ancestral
minn <- apply(cbind(ss1,ss2),1,min)
maxx <- apply(cbind(BBB2$GLOBAL_INIT3,ss3),1,max)
FwPrp <- minn - maxx
###################
new_L <- BBB2$GLOBAL_INIT3
for(xx in 1:BBB2$popnum){
new_L <- new_L + lgamma(samp_pop_freq[[xx]][1,] + theta[,xx] * samp_freq_new[1,]) -
lfactorial (samp_pop_freq[[xx]][1,]) - lgamma(theta[,xx]*samp_freq_new[1,])+
lgamma(samp_pop_freq[[xx]][2,]+theta[,xx]*samp_freq_new[2,]) -
lfactorial(samp_pop_freq[[xx]][2,]) - lgamma(theta[,xx]*samp_freq_new[2,])
}
A <- -old_L + new_L + log(BwPrp) - log(FwPrp) +
(lambda-1)*(log(samp_freq_new[1,]) + log(samp_freq_new[2,]) - log(old_m) - log(old_n) )
r <- runif(BBB2$locnum)
check <- log(r) < A
#for(xx in which(check)){
# BBB2$freq_locus[xx,samp[1,xx]] <- samp_freq_new[1,xx]
# BBB2$freq_locus[xx,samp[2,xx]] <- samp_freq_new[2,xx]
# BBB2$acc_freq_ancestral[xx] <- BBB2$acc_freq_ancestral[xx] + 1
#}
for(xx in which(check)){
BBB2$freq_locus[xx,1] <- samp_freq_new[1,xx]
BBB2$freq_locus[xx,2] <- samp_freq_new[2,xx]
BBB2$acc_freq_ancestral[xx] <- BBB2$acc_freq_ancestral[xx] + 1
}
}# End of Function
X_jump_model_codominant <- function(){
old_alpha <- BBB2$d_alpha
excluded <- which(!BBB2$alpha_included)
BBB2$d_alpha[BBB2$alpha_included] <- 0 #rnorm(1,0,1)
#cat(1)
if(length(excluded)!=0){
### mod -------------------------------------------
whatgroups <- unique(BBB2$GROUP[excluded])
group_id <- match(whatgroups,BBB2$GROUP)
e_group <- rnorm(length(whatgroups),mean=BBB2$mean_alpha[group_id],sd=sqrt(BBB2$var_alpha[group_id]))
#e_group <- rnorm(length(whatgroups),mean=BBB2$d_alpha[group_id],1)
#e_group <- rlaplace(length(whatgroups),BBB2$mean_alpha[group_id],BBB2$var_alpha[group_id])
.Call("update_group_Jump",
BBB2$d_alpha,
as.integer(BBB2$GROUP),
as.integer(whatgroups),
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)
#+ tapply(oldL,BBB2$GROUP,min)
#- log(1-exp(tapply(oldL,BBB2$GROUP,max)))
## 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)
#+ tapply(newL,BBB2$GROUP,min)
#- log(1-exp(tapply(newL,BBB2$GROUP,max)))
A_group <- -oldL + newL
A <- A_group[BBB2$GROUP]
included <- BBB2$alpha_included #==TRUE
excluded <- !included
A_included <- 0
A_excluded <- 0
# include in the selection model
if(any(included)){
A_included <- - log(BBB2$tempting)
- log(dnorm(BBB2$d_alpha[included],BBB2$mean_alpha[included],sqrt(BBB2$var_alpha[included])))
+ log(dnorm(BBB2$d_alpha[included],0,1))
}
# exclude from the selection model
if(any(excluded)){
A_excluded <- log(BBB2$tempting)
+ log(dnorm(BBB2$d_alpha[excluded],BBB2$mean_alpha[excluded],sqrt(BBB2$var_alpha[excluded])))
- log(dnorm(BBB2$d_alpha[excluded],0,1))
}
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)
# tempting
# BBB2$acc.ratio <- BBB2$acc.ratio + length(which(!check))/length(A)
#
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
X_update_d_alphaco_i_new <- 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,sub_GROUP,sum)
#+ tapply(oldL,sub_GROUP,min)
#- log(1-exp(tapply(oldL,sub_GROUP,max)))
e_group <- rnorm(BBB2$N.REGIONS,0,tapply(BBB2$var_prop_alpha,BBB2$GROUP,unique)) # BUG FIXED
#e_group <- rnorm(BBB2$N.REGIONS,0,1)
new_alpha <- old_alpha + e_group[BBB2$GROUP]
whatgroup <- 1:BBB2$N.REGIONS
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,sub_GROUP,sum)
#+ tapply(newL,sub_GROUP,min)
#- log(1-exp(tapply(newL,sub_GROUP,max)))
old_alpha_group <- tapply(old_alpha,sub_GROUP,function(x){return(x[1])})
new_alpha_group <- tapply(new_alpha,sub_GROUP,function(x){return(x[1])})
A <- -oldL + newL + (old_alpha_group^2-new_alpha_group^2)/(2*BBB2$sd_prior_alpha^2)
r <- runif(length(whatgroup))
check <- log(r) < A
group_idX <- which(check)
group_id <- whatgroup[group_idX]
if(length(group_idX)>0){
#CCode
.Call("update_group2",
as.integer(group_id),
as.integer(BBB2$GROUP),
BBB2$d_alpha,
new_alpha,
BBB2$acc_alpha
)
}
}# End of function
X_update_d_alphaco_i_new_id <- function(id){
sub_freq_locus <- BBB2$freq_locus[id,]
sub_freq_pop <- lapply(BBB2$freq_pop,function(yy){return(yy[id,])})
old_alpha <- BBB2$d_alpha[id]
new_alpha <- old_alpha
sub_sample_size <- BBB2$sample_size[id,]
sub_GROUP <- BBB2$GROUP[id]
locnum2 <- length(id)
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,sub_GROUP,sum)
# create new alpha
group_id <- unique(sub_GROUP) # what groups have to be updated
xxx <- BASIX.match(group_id,BBB2$GROUP)
e_group <- rnorm(length(group_id),0,sqrt(BBB2$var_prop_alpha[xxx]))
# update the whole groups
.Call("update_groupJumpX",
new_alpha,
as.integer(sub_GROUP),
as.integer(group_id),
e_group
)
#whatgroup <- 1:BBB2$N.REGIONS
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,sub_GROUP,sum)
old_alpha_group <- tapply(old_alpha,sub_GROUP,function(x){return(x[1])})
new_alpha_group <- tapply(new_alpha,sub_GROUP,function(x){return(x[1])})
A <- -oldL + newL + (old_alpha_group^2-new_alpha_group^2)/(2*BBB2$sd_prior_alpha^2)
r <- runif(length(group_id))
check <- log(r) < A
group_idX <- which(check)
group_id <- group_id[group_idX]
if(length(group_idX)>0){
#CCode
.Call("update_group2",
as.integer(group_id),
as.integer(BBB2$GROUP),
BBB2$d_alpha,
new_alpha,
BBB2$acc_alpha
)
}
}# End of function update alpha new id
pievos101/BlockFeST documentation built on March 16, 2021, 2:18 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions X_update_d_alphaco_i_new_id X_update_d_alphaco_i_new X_jump_model_codominant update_freq_codominant update_d_betaco X_update_d_alphaco_i
#source("math.R")
#constants
#e_f <- 0.05
#sd_prior_beta <- 1 #3
#mean_prior_beta <- -1 #-5.15
#lambda <- 1
X_update_d_alphaco_i <- function(id){ # id
if(id[1]==FALSE){ # Pilot runs
id <- 1:BBB2$locnum
locnum2 <- BBB2$locnum
}else{ # Main Loop
locnum2 <- length(id)
}
##
sub_freq_locus <- BBB2$freq_locus[id,]
sub_freq_pop <- lapply(BBB2$freq_pop,function(yy){return(yy[id,])})
old_alpha <- BBB2$d_alpha[id]
sub_sample_size <- BBB2$sample_size[id,]
sub_GROUP <- BBB2$GROUP[id]
##
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)
## modification ------------------------------------
oldL <- tapply(oldL,sub_GROUP,sum)
##--------------------------------------------------
# e <- rnorm(locnum2,0,sqrt(BBB2$var_prop_alpha))
# new_alpha <- old_alpha + e
### modification -----------------------------------------------
if(!BBB2$PILOT){
whatgroup <- unique(sub_GROUP)
xxx <- BASIX.match(whatgroup,BBB2$GROUP)
e_group <- rnorm(length(whatgroup),0,sqrt(BBB2$var_prop_alpha[xxx]))
#
BBB2$new_alpha <- old_alpha
BBB2$jjj <- 1
sapply(whatgroup,function(x){
what <- x==sub_GROUP
BBB2$new_alpha[what] <- old_alpha[what] + e_group[BBB2$jjj]
BBB2$jjj <- BBB2$jjj + 1
})
new_alpha <- BBB2$new_alpha
#
}else{
e_group <- rnorm(BBB2$N.REGIONS,0,sqrt(BBB2$var_prop_alpha))
new_alpha <- old_alpha + e_group[BBB2$GROUP]
whatgroup <- 1:BBB2$N.REGIONS
}
#
##--------------------------------------------------------------
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)
## modification ----------------------------
newL <- tapply(newL,sub_GROUP,sum)
old_alpha_group <- tapply(old_alpha,sub_GROUP,function(x){return(x[1])})
new_alpha_group <- tapply(new_alpha,sub_GROUP,function(x){return(x[1])})
# print(old_alpha)
# print(new_alpha)
##-----------------------------------------
#A <- -oldL + newL + (old_alpha^2-new_alpha^2)/(2*BBB2$sd_prior_alpha^2)
## mod ------------------------
A <- -oldL + newL + (old_alpha_group^2-new_alpha_group^2)/(2*BBB2$sd_prior_alpha^2)
##-----------------------------
# r <- runif(locnum2)
## mod ---------------------------
r <- runif(length(whatgroup))
## mod ---------------------------
check <- log(r) < A
# BBB2$d_alpha[id[check]] <- new_alpha[check]
# BBB2$acc_alpha[id[check]] <- BBB2$acc_alpha[id[check]] + 1
# modification -------------------------------------
# update groups
group_idX <- which(check)
group_id <- whatgroup[group_idX]
# make this also in C
if(length(group_idX)>0){
#res <- .Call("update_group2",group_id,BBB2$GROUP,BBB2$d_alpha,new_alpha,BBB2$acc_alpha)
#BBB2$d_alpha <- res[[1]]
#BBB2$acc_alpha <- res[[2]]
mm <- sapply(group_id,function(xx){xx==BBB2$GROUP})
BBB2$iii <- 1
apply(mm,2,function(x){
BBB2$d_alpha[x] <- new_alpha_group[group_idX[BBB2$iii]] # new_alpha[x]
BBB2$acc_alpha[x] <- BBB2$acc_alpha[x] + 1
BBB2$iii <- BBB2$iii + 1
})
}
# ---------------------------------------------------
}# End of Function
update_d_betaco <- function(){
# Init
#val <- rep(NA,2)
val <- BBB2$GLOBAL_INIT2
# ---------------
old_beta <- BBB2$d_beta
old_theta <- outer(BBB2$d_alpha,old_beta,"+")
old_theta <- exp(-(old_theta))
#cat("old_theta", old_theta,"\n")
old1 <- lgamma(old_theta)-lgamma(BBB2$sample_size+old_theta)
#cat("old1: ",old1,"\n")
for(xx in 1:BBB2$popnum){
xyz <- old_theta[,xx]*BBB2$freq_locus
val[xx] <- sum(lgamma(BBB2$freq_pop[[xx]]+xyz)-lgamma(xyz),na.rm=TRUE)
}
#cat("val:",val,"\n")
old1 <- colSums(old1)
oldL <- old1 + val
#--------- calculate new likelihood
e <- rnorm(BBB2$popnum,0,BBB2$var_prop_beta)
BBB2$d_beta <- BBB2$d_beta + e
new_theta <- outer(BBB2$d_alpha,BBB2$d_beta,"+")
new_theta <- exp(-(new_theta))
new1 <- lgamma(new_theta)-lgamma(BBB2$sample_size+new_theta)
for(xx in 1:BBB2$popnum){
xyz <- new_theta[,xx]*BBB2$freq_locus
val[xx] <- sum(lgamma(BBB2$freq_pop[[xx]]+xyz)-lgamma(xyz),na.rm=TRUE)
}
new1 <- colSums(new1)
newL <- new1 + val
A <- -oldL + newL
#+ (old_beta -BBB2$d_beta)*(old_beta+BBB2$d_beta-2*BBB2$mean_prior_beta)/(2*BBB2$sd_prior_beta^2)
#+ (old_beta^2-BBB2$d_beta^2)/(2*BBB2$sd_prior_beta^2)
#+ (old_alpha_group^2-new_alpha_group^2)/(2*BBB2$sd_prior_alpha^2)
r <- runif(BBB2$popnum)
check <- log(r)>A
#cat("-oldL: ",-oldL, "\n")
#cat("newL: ",newL, "\n")
#print(A)
#print(old_beta)
#print(check)
#print(r)
#if(any(is.na(check))){
#print("yes")
#check[is.na(check)] <- TRUE
#}
BBB2$d_beta[check] <- old_beta[check]
BBB2$acc_beta[!check] <- BBB2$acc_beta[!check] + 1
}# End of Function
update_freq_codominant <- function(){
lambda <- 1
# choose randomly two haplotype of each locus
# samp <- sapply(BBB2$hapcount,function(x){return(sample(1:x,size=2)) })
theta <- outer(BBB2$d_alpha,BBB2$d_beta,"+")
theta <- exp(-theta)
# ok
#BBB2$iii <- 1
# frequencies of haplotypes of the two random haps
#samp_freq <- apply(samp,2,function(x){
# r <- BBB2$freq_locus[BBB2$iii,x]
# BBB2$iii <- BBB2$iii+1
# return(r)
# })
#ok
## mod
samp_freq <- t(BBB2$freq_locus)
## end mod
#print(dim(samp_freq))
#print(dim(BBB2$freq_locus))
# haplotype distribution for each pop of the two random haps
# BBB2$iii <- 1
samp_pop_freq <- vector("list",BBB2$popnum)
for(yy in 1:BBB2$popnum){
samp_pop_freq[[yy]] <- t(BBB2$freq_pop[[yy]])
#samp_pop <- apply(samp,2,function(x){
# r <- BBB2$freq_pop[[yy]][BBB2$iii,x]
# BBB2$iii <- BBB2$iii+1
# return(r)
# })
#samp_pop_freq[[yy]] <- samp_pop
#BBB2$iii <- 1
}
#ok
#print(dim(samp_pop_freq[[1]]))
old_L <- BBB2$GLOBAL_INIT3
# 1 <- m
# 2 <- n
for(xx in 1:BBB2$popnum){
old_L <- old_L + lgamma(samp_pop_freq[[xx]][1,] + theta[,xx] * samp_freq[1,]) -
lfactorial (samp_pop_freq[[xx]][1,]) - lgamma(theta[,xx]*samp_freq[1,])+
lgamma(samp_pop_freq[[xx]][2,]+theta[,xx]*samp_freq[2,]) -
lfactorial(samp_pop_freq[[xx]][2,]) - lgamma(theta[,xx]*samp_freq[2,])
}
old_m <- samp_freq[1,]
old_n <- samp_freq[2,]
#### BwPrp
old_m_n <- old_m + old_n
old_m_ancestral <- old_m + BBB2$e_ancestral
old_m_ancestral2 <- old_m - BBB2$e_ancestral
checkmin <- cbind(old_m_n,old_m_ancestral)
checkmax <- cbind(BBB2$GLOBAL_INIT3,old_m_ancestral2)
minn <- apply(checkmin,1,min)
maxx <- apply(checkmax,1,max)
BwPrp <- minn - maxx
###################
# change frequencies
u <- runif(BBB2$locnum,0,BwPrp) + maxx
# samp_freq_new
samp_freq_new <- samp_freq
samp_freq_new[1,] <- u
samp_freq_new[2,] <- old_n + old_m - u
### FwPrp
ss1 <- colSums(samp_freq_new)
ss2 <- samp_freq_new[1,] + BBB2$e_ancestral
ss3 <- samp_freq_new[1,] - BBB2$e_ancestral
minn <- apply(cbind(ss1,ss2),1,min)
maxx <- apply(cbind(BBB2$GLOBAL_INIT3,ss3),1,max)
FwPrp <- minn - maxx
###################
new_L <- BBB2$GLOBAL_INIT3
for(xx in 1:BBB2$popnum){
new_L <- new_L + lgamma(samp_pop_freq[[xx]][1,] + theta[,xx] * samp_freq_new[1,]) -
lfactorial (samp_pop_freq[[xx]][1,]) - lgamma(theta[,xx]*samp_freq_new[1,])+
lgamma(samp_pop_freq[[xx]][2,]+theta[,xx]*samp_freq_new[2,]) -
lfactorial(samp_pop_freq[[xx]][2,]) - lgamma(theta[,xx]*samp_freq_new[2,])
}
A <- -old_L + new_L + log(BwPrp) - log(FwPrp) +
(lambda-1)*(log(samp_freq_new[1,]) + log(samp_freq_new[2,]) - log(old_m) - log(old_n) )
r <- runif(BBB2$locnum)
check <- log(r) < A
#for(xx in which(check)){
# BBB2$freq_locus[xx,samp[1,xx]] <- samp_freq_new[1,xx]
# BBB2$freq_locus[xx,samp[2,xx]] <- samp_freq_new[2,xx]
# BBB2$acc_freq_ancestral[xx] <- BBB2$acc_freq_ancestral[xx] + 1
#}
for(xx in which(check)){
BBB2$freq_locus[xx,1] <- samp_freq_new[1,xx]
BBB2$freq_locus[xx,2] <- samp_freq_new[2,xx]
BBB2$acc_freq_ancestral[xx] <- BBB2$acc_freq_ancestral[xx] + 1
}
}# End of Function
X_jump_model_codominant <- function(){
old_alpha <- BBB2$d_alpha
excluded <- which(!BBB2$alpha_included)
BBB2$d_alpha[BBB2$alpha_included] <- 0 #rnorm(1,0,1)
#cat(1)
if(length(excluded)!=0){
### mod -------------------------------------------
whatgroups <- unique(BBB2$GROUP[excluded])
group_id <- match(whatgroups,BBB2$GROUP)
e_group <- rnorm(length(whatgroups),mean=BBB2$mean_alpha[group_id],sd=sqrt(BBB2$var_alpha[group_id]))
#e_group <- rnorm(length(whatgroups),mean=BBB2$d_alpha[group_id],1)
#e_group <- rlaplace(length(whatgroups),BBB2$mean_alpha[group_id],BBB2$var_alpha[group_id])
.Call("update_group_Jump",
BBB2$d_alpha,
as.integer(BBB2$GROUP),
as.integer(whatgroups),
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)
#+ tapply(oldL,BBB2$GROUP,min)
#- log(1-exp(tapply(oldL,BBB2$GROUP,max)))
## 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)
#+ tapply(newL,BBB2$GROUP,min)
#- log(1-exp(tapply(newL,BBB2$GROUP,max)))
A_group <- -oldL + newL
A <- A_group[BBB2$GROUP]
included <- BBB2$alpha_included #==TRUE
excluded <- !included
A_included <- 0
A_excluded <- 0
# include in the selection model
if(any(included)){
A_included <- - log(BBB2$tempting)
- log(dnorm(BBB2$d_alpha[included],BBB2$mean_alpha[included],sqrt(BBB2$var_alpha[included])))
+ log(dnorm(BBB2$d_alpha[included],0,1))
}
# exclude from the selection model
if(any(excluded)){
A_excluded <- log(BBB2$tempting)
+ log(dnorm(BBB2$d_alpha[excluded],BBB2$mean_alpha[excluded],sqrt(BBB2$var_alpha[excluded])))
- log(dnorm(BBB2$d_alpha[excluded],0,1))
}
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)
# tempting
# BBB2$acc.ratio <- BBB2$acc.ratio + length(which(!check))/length(A)
#
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
X_update_d_alphaco_i_new <- 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,sub_GROUP,sum)
#+ tapply(oldL,sub_GROUP,min)
#- log(1-exp(tapply(oldL,sub_GROUP,max)))
e_group <- rnorm(BBB2$N.REGIONS,0,tapply(BBB2$var_prop_alpha,BBB2$GROUP,unique)) # BUG FIXED
#e_group <- rnorm(BBB2$N.REGIONS,0,1)
new_alpha <- old_alpha + e_group[BBB2$GROUP]
whatgroup <- 1:BBB2$N.REGIONS
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,sub_GROUP,sum)
#+ tapply(newL,sub_GROUP,min)
#- log(1-exp(tapply(newL,sub_GROUP,max)))
old_alpha_group <- tapply(old_alpha,sub_GROUP,function(x){return(x[1])})
new_alpha_group <- tapply(new_alpha,sub_GROUP,function(x){return(x[1])})
A <- -oldL + newL + (old_alpha_group^2-new_alpha_group^2)/(2*BBB2$sd_prior_alpha^2)
r <- runif(length(whatgroup))
check <- log(r) < A
group_idX <- which(check)
group_id <- whatgroup[group_idX]
if(length(group_idX)>0){
#CCode
.Call("update_group2",
as.integer(group_id),
as.integer(BBB2$GROUP),
BBB2$d_alpha,
new_alpha,
BBB2$acc_alpha
)
}
}# End of function
X_update_d_alphaco_i_new_id <- function(id){
sub_freq_locus <- BBB2$freq_locus[id,]
sub_freq_pop <- lapply(BBB2$freq_pop,function(yy){return(yy[id,])})
old_alpha <- BBB2$d_alpha[id]
new_alpha <- old_alpha
sub_sample_size <- BBB2$sample_size[id,]
sub_GROUP <- BBB2$GROUP[id]
locnum2 <- length(id)
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,sub_GROUP,sum)
# create new alpha
group_id <- unique(sub_GROUP) # what groups have to be updated
xxx <- BASIX.match(group_id,BBB2$GROUP)
e_group <- rnorm(length(group_id),0,sqrt(BBB2$var_prop_alpha[xxx]))
# update the whole groups
.Call("update_groupJumpX",
new_alpha,
as.integer(sub_GROUP),
as.integer(group_id),
e_group
)
#whatgroup <- 1:BBB2$N.REGIONS
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,sub_GROUP,sum)
old_alpha_group <- tapply(old_alpha,sub_GROUP,function(x){return(x[1])})
new_alpha_group <- tapply(new_alpha,sub_GROUP,function(x){return(x[1])})
A <- -oldL + newL + (old_alpha_group^2-new_alpha_group^2)/(2*BBB2$sd_prior_alpha^2)
r <- runif(length(group_id))
check <- log(r) < A
group_idX <- which(check)
group_id <- group_id[group_idX]
if(length(group_idX)>0){
#CCode
.Call("update_group2",
as.integer(group_id),
as.integer(BBB2$GROUP),
BBB2$d_alpha,
new_alpha,
BBB2$acc_alpha
)
}
}# End of function update alpha new id
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.