waterman_function(ps,n, d, k)
{
#this program compute the probability of k-interrupt-run, P(Lk>=d), 
# in a iid Bernoulli seq with P(1)=ps, 
# where n is the length of seq., 
# d is the length of Lk, and k is the num. of mismached.
# output 3 values: est. of the prob, and its lower, upper bounds
	s <- d - k
	EEW <- (n * (s/d - ps) * dbinom(s, d, ps)) * (-1)
	EEW <- exp(EEW)
	erro <- 7 * d * dbinom(s, d, ps) + (1 - pbinom(s, d, ps))
	upb <- EEW + erro
	lob <- EEW - erro
	bd <- c(EEW, upb, lob)
	bd <- 1 - bd
    if (bd>1) bd_1
    bd
}