d <- read.table('inp_RF_M.txt',skip=1) vec_q <- sort(d[,1]) nn <- length(vec_q) alp <- 0.4 #Cunnane formula vec_p <- ppoints(nn,alp) #plot(vec_q,vec_p) fx <- function(bb,vec_q){ bbq <- bb*vec_q sum1 <- sum(sqrt(bbq)) sum2 <- sum(bbq*exp(-sqrt(bbq))) sum3 <- sum(exp(-sqrt(bbq))) sum4 <- sum(sqrt(bbq)*exp(-sqrt(bbq))) difa <- (sum1-2*nn)/sum2-nn/(sum3+sum4) return (difa) } eps <-0.0001 b1 <- (2*nn/sum(sqrt(vec_q)))^2+eps b2 <- b1+0.5 bb <- 0.5*(b1+b2) f1 <- fx(b1,vec_q) f2 <- fx(b2,vec_q) ff <- fx(bb,vec_q) while(eps