Temat: Analiza dyskryminacyjna - lda {MASS}

Witam,
interesuje mnie wykorzystanie funkcji lda z pakietu MASS do klasyfikacji obiektów pochodzących z dwóch populacji. Z tego co udało mi się ustalić, funkcja zwraca wektor kanoniczny, a nie postać standardową funkcji klasyfikacyjnej z wyrazem wolnym. Chciałbym dowiedzieć się, jak wygląda procedura klasyfikacji obiektu z wykorzystaniem postaci kanonicznej funkcji oraz jak ustalana jest wartość progowa przynależności obiektu do populacji ? Na podstawie opracowań z sieci starałem się to wyznaczyć, ale wynik klasyfikacji, które uzyskuję nie pokrywają się z tymi wyrzucanymi przez funkcję predict dla obiektu lda.

Pozdrawiam,
BartekTen post został edytowany przez Autora dnia 20.10.15 o godzinie 11:38

Temat: Analiza dyskryminacyjna - lda {MASS}

W takich sytuacjach warto rzucić okiem na postać konkretnej funkcji predict (tu: predict.lda):

> getAnywhere("predict.lda")

A single object matching ‘predict.lda’ was found

It was found in the following places

  registered S3 method for predict from namespace MASS

  namespace:MASS

with value



function (object, newdata, prior = object$prior, dimen, method = c("plug-in", 

    "predictive", "debiased"), ...) 

{

    if (!inherits(object, "lda")) 

        stop("object not of class \"lda\"")

    if (!is.null(Terms <- object$terms)) {

        Terms <- delete.response(Terms)

        if (missing(newdata)) 

            newdata <- model.frame(object)

        else {

            newdata <- model.frame(Terms, newdata, na.action = na.pass, 

                xlev = object$xlevels)

            if (!is.null(cl <- attr(Terms, "dataClasses"))) 

                .checkMFClasses(cl, newdata)

        }

        x <- model.matrix(Terms, newdata, contrasts = object$contrasts)

        xint <- match("(Intercept)", colnames(x), nomatch = 0L)

        if (xint > 0L) 

            x <- x[, -xint, drop = FALSE]

    }

    else {

        if (missing(newdata)) {

            if (!is.null(sub <- object$call$subset)) 

                newdata <- eval.parent(parse(text = paste(deparse(object$call$x, 

                  backtick = TRUE), "[", deparse(sub, backtick = TRUE), 

                  ",]")))

            else newdata <- eval.parent(object$call$x)

            if (!is.null(nas <- object$call$na.action)) 

                newdata <- eval(call(nas, newdata))

        }

        if (is.null(dim(newdata))) 

            dim(newdata) <- c(1L, length(newdata))

        x <- as.matrix(newdata)

    }

    if (ncol(x) != ncol(object$means)) 

        stop("wrong number of variables")

    if (length(colnames(x)) > 0L && any(colnames(x) != dimnames(object$means)[[2L]])) 

        warning("variable names in 'newdata' do not match those in 'object'")

    ng <- length(object$prior)

    if (!missing(prior)) {

        if (any(prior < 0) || round(sum(prior), 5) != 1) 

            stop("invalid 'prior'")

        if (length(prior) != ng) 

            stop("'prior' is of incorrect length")

    }

    means <- colSums(prior * object$means)

    scaling <- object$scaling

    x <- scale(x, center = means, scale = FALSE) %*% scaling

    dm <- scale(object$means, center = means, scale = FALSE) %*% 

        scaling

    method <- match.arg(method)

    dimen <- if (missing(dimen)) 

        length(object$svd)

    else min(dimen, length(object$svd))

    N <- object$N

    if (method == "plug-in") {

        dm <- dm[, 1L:dimen, drop = FALSE]

        dist <- matrix(0.5 * rowSums(dm^2) - log(prior), nrow(x), 

            length(prior), byrow = TRUE) - x[, 1L:dimen, drop = FALSE] %*% 

            t(dm)

        dist <- exp(-(dist - apply(dist, 1L, min, na.rm = TRUE)))

    }

    else if (method == "debiased") {

        dm <- dm[, 1L:dimen, drop = FALSE]

        dist <- matrix(0.5 * rowSums(dm^2), nrow(x), ng, byrow = TRUE) - 

            x[, 1L:dimen, drop = FALSE] %*% t(dm)

        dist <- (N - ng - dimen - 1)/(N - ng) * dist - matrix(log(prior) - 

            dimen/object$counts, nrow(x), ng, byrow = TRUE)

        dist <- exp(-(dist - apply(dist, 1L, min, na.rm = TRUE)))

    }

    else {

        dist <- matrix(0, nrow = nrow(x), ncol = ng)

        p <- ncol(object$means)

        X <- x * sqrt(N/(N - ng))

        for (i in 1L:ng) {

            nk <- object$counts[i]

            dev <- scale(X, center = dm[i, ], scale = FALSE)

            dev <- 1 + rowSums(dev^2) * nk/(N * (nk + 1))

            dist[, i] <- prior[i] * (nk/(nk + 1))^(p/2) * dev^(-(N - 

                ng + 1)/2)

        }

    }

    posterior <- dist/drop(dist %*% rep(1, ng))

    nm <- names(object$prior)

    cl <- factor(nm[max.col(posterior)], levels = object$lev)

    dimnames(posterior) <- list(rownames(x), nm)

    list(class = cl, posterior = posterior, x = x[, 1L:dimen, 

        drop = FALSE])

}

<bytecode: 0x0000000013c92988>

<environment: namespace:MASS>

Temat: Analiza dyskryminacyjna - lda {MASS}

Dziękuję bardzo. Funkcję getAnywhere na pewno zapamiętam :)
Pozdrawiam