#
 Exp <- function(n, r, N = 256, k = 1, wav = "d8", ...)
{
#
# This S+ function calculates the expectation of
# r:n-th order statistics for a standard normal
# distribution.
# 
# This function needs S+Wavelets [procedure  dwt].
        ret <- list()
        eps <- 1/(2 * N)
        ret$x <- seq(eps, 1 - eps, length = N)
        #-------define lik function------------------------
        fun <- function(x, n, r)
        {
                return(exp((r - 1) * log(x) + (n - r) * log(1 - x) - lgamma(r) - 
                        lgamma(n - r + 1) + lgamma(n + 1)))
        }
#---------------------------------------------------
        ret$g <- fun(ret$x, n, r)
        ret$g <- ret$g/sum(ret$g)
        ret$h <- qnorm(ret$x)^k
        ret$gd <- as.vector(dwt(ret$g, wavelet = wav, ...))
        ret$hd <- as.vector(dwt(ret$h, wavelet = wav, ...))
        #--------------------------------------------------
        ret$res <- print(sum(ret$gd * ret$hd), digits = 18)
        #---------------------------------
        return(ret)
}

#----------------------------------------------------
#
#Examples:
#
#> module(wavelets)
#> a <- Exp( n=5, r=2, N=2^19)
#[1] -0.49501897047558691
#
#> a <- Exp(n=10000000000, r=4000000000, N=2^19)
#[1] -0.25334710323084308
#