https://github.com/cran/pracma
Raw File
Tip revision: 708a2ad382a163d1eef5af0665e3ae2aad200ced authored by HwB on 21 March 2013, 00:00:00 UTC
version 1.4.5
Tip revision: 708a2ad
agmean.R
##
##  a g m e a n . R  Arithmetic-geometric Mean
##


agmean <- function(a, b) {
    eps <- .Machine$double.eps
    stopifnot(is.numeric(a) || is.complex(a),
              is.numeric(b) || is.complex(b))
    if (is.numeric(a) && any(a < 0) || is.numeric(b) && any(b < 0)) {
        a <- as.complex(a)
        b <- as.complex(b)
    }

    if (length(a) == 1) {
        n <- length(b)
        a <- rep(a, n)
    } else if (length(b) == 1) {
        n <- length(a)
        b <- rep(b, n)
    } else 
        if (length(a) != length(b)) 
            stop("Arguments must have the same length or one has length 1.")

    while ( any(abs(a-b) >= eps) ) {
        a1 <- (a + b) / 2
        b1 <- sqrt(a * b)
        if (max(abs(a-a1)) < eps && max(abs(b-b1)) < eps) break
        a <- a1
        b <- b1
    }

    return( (a+b)/2 )
}
back to top