R/init_coef.R

In PopGenome: An Efficient Swiss Army Knife for Population Genomic Analyses

init_coef <- function (sample_size,p){

#Tajima and Fu coefficients... and Achaz coefficients. Total  <- 23 Most from Simonsen 1995*/
		
       p <- rep(-1,24)
   
		
      #an <- sum(1/1:(n-1))
      #bn <- sum(1/1:((n-1)*(n-1)))
      
	    an <-0
            bn <- 0.0
	    n  <- sample_size
	
	 #schneller !!!!
   if(n>1){
          bn <- 1/((1:(n-1))*(1:(n-1)))
          bn <- sum(bn)
          an <- 1/(1:(n-1))
	  an <- sum(an)
	 }
	#for(x in 1:(n-1)) {
	#	an <- an + 1.0/x
	#	bn <- bn + 1.0/(x*x)
	#}
	
	p[1] <- an
	p[2] <- bn
	
	# vt */
	p[3] <- (2.0*(n*n + n + 3.0)/(9.0*n*(n-1.0))  -(n+2.0)/(an * n) + bn/(an * an)) / (an*an + bn)	
	# ut */
	p[4] <- (( (n+1.0)/(3.0*(n-1.0)) - 1.0/an)/an) - p[3]
	# vd* */
	p[5] <- (bn/(an*an) - 2.0/n *(1.0 + 1.0/an - an + an/n) - 1.0/(n*n))/ (an*an + bn)
	# ud* */
	p[6] <- (((n-1.0)/n - 1.0/an) / an) - p[5]
 
 #/* vf* */
	p[7] <- (((2.0*n*n*n + 110.0 * n*n - 255.0 * n + 153.0)/ (9.0 * n * n * (n-1.0)) + (2.0*(n-1.0) *an)/ (n*n) - (8.0 * bn)/n)) / (an*an + bn)
 #/* uf* */
	p[8] <- (((4.0*n*n + 19.0*n + 3.0 - 12.0 * (n +1.0) * (an + 1.0/n))/ (3.0*n *(n-1.0))) / an) - p[7]	
  
  # vf* */
#	p[7] <- (((2.0*n*n*n + 110.0 * n*n - 255.0 * n + 153.0) / (9.0 * n * n * (n-1.0)) +(2.0*(n-1.0)*an)/ (n*n) - (8.0 * bn)/n)) / (an*an + bn)
	# uf* */
#	p[8] <- (((4.0*n*n + 19.0*n + 3.0 - 12.0 * (n +1.0) * (an + 1.0/n)) / (3.0*n *(n-1.0))) / an) - p[7]
	
  # cn */
	cn <- 2.0 * (n*an-2.0*(n-1.0)) / ((n-1.0)*(n-2.0))
	# vd */
	p[9] <- 1.0 + (an*an/(bn+an*an)) * (cn - ((n+1.0)/(n-1.0)))
	# ud* */
	p[10] <- an -1.0 - p[9]
	# vf */
	p[11] <- (cn + (2.0*(n*n+n+3.0))/(9.0*n*(n-1.0)) - 2.0/(n-1.0)) / (an*an + bn)
	# uf */
	p[12] <- (1.0 + (n+1.0)/(3.0*(n-1.0)) - 4*(n+1.0) /((n-1.0)*(n-1.0)) * (an+(1.0/n) - 2.0*n/(n+1.0)))/ an  -  p[11]
	# Achaz f */
	p[13] <- (n - 2.0)/(n * (an - 1.0))
	# Achaz f* */
	p[14] <- (n - 3.0)/(an * (n - 1.0) - n)
	# alfa_n */
	p[15] <- p[13] * p[13] * (an - 1.0) + p[13] * (an * (4.0*(n+1.0))/((n-1.0)*(n-1.0)) - (2.0 *(n + 1.0)*(n + 2.0))/(n * (n-1.0))) -an * (8.0 * (n+1.0)/(n*(n-1.0)*(n-1.0))) +(n*n*n + n*n + 60.0*n + 12.0)/(3.0*n*n*(n-1.0))
	
	# beta_n */
	p[16] <- p[13] * p[13] * (bn + an * (4.0/((n - 1.0) * (n - 2.0))) - 4.0/(n-2.0)) +p[13] * (-an * (4.0*(n+2.0))/(n*(n-1.0)*(n-2.0)) - (n*n*n - 3.0*n*n - 16.0*n + 20.)/(n*(n-1.0)*(n-2.0))) +an * (8.0/(n*(n-1.0)*(n-2.0))) +2.0 * (n*n*n*n - n*n*n - 17.0*n*n - 42.0*n + 72.0)/(9.0*n*n*(n-1.0)*(n-2.0))
	
	# alfa*_n */
	p[17] <- p[14] * p[14] * (an - n/(n-1.0)) + p[14] * (an * 4.0 * (n+1.0)/((n-1.0)*(n-1.0)) - 2.0 * (n+3.0)/(n-1.0)) -an * (8.0*(n+1.0))/(n*(n-1.0)*(n-1.0)) +(n*n + n + 60.0)/(3.0*n*(n-1.0))
	
	# beta*_n */
	p[18] <- p[14] * p[14] * (bn - (2.0*n-1.0)/((n-1.0)*(n-1.0))) +p[14] * (bn * 8.0/(n-1.0) - an * 4.0/(n*(n-1.0)) - (n*n*n + 12.0*n*n - 35.0*n +18.0)/(n*(n-1.0)*(n-1.0))) -bn * 16.0/(n*(n-1.0)) +an * 8.0/(n*n*(n-1.0)) +2.0 * (n*n*n*n + 110.0*n*n - 255.0*n + 126.0)/(9.0*n*n*(n-1.0)*(n-1.0))
	
	#b1*/
	p[19] <-(n+1.0)/(3.0*(n-1.0))
	#b2*/
	p[20] <-2.0*(n*n + n + 3)/(9*n*(n-1))
	#c1*/
	p[21] <- p[19] - 1.0/p[1]
	#c2*/
	p[22] <-p[20] - (n+2.0)/(p[1]*n) + (p[2]/(p[1]*p[1]))
	#e1*/
	p[23] <-p[21]/p[1]
	#e2*/
	p[24] <-p[22]/(p[1]*p[1] + p[2])


return(p)

}