ourMELONS/R/computePopulationLogml.R

36 lines
1.1 KiB
R
Raw Normal View History

2021-02-15 09:49:21 +01:00
computePopulationLogml <- function(pops, adjprior, priorTerm) {
2020-10-19 14:08:25 +02:00
# Palauttaa length(pops)*1 taulukon, jossa on laskettu korikohtaiset
2021-02-01 09:22:58 +01:00
# ======================================================== #
# Limiting COUNTS size #
# ======================================================== #
COUNTS <- COUNTS[seq_len(nrow(adjprior)), seq_len(ncol(adjprior)), pops, drop=FALSE]
2020-10-19 14:08:25 +02:00
x <- size(COUNTS, 1)
y <- size(COUNTS, 2)
z <- length(pops)
2021-02-01 09:22:58 +01:00
# ======================================================== #
# Computation #
# ======================================================== #
2021-02-01 10:01:02 +01:00
isarray <- length(dim(repmat(adjprior, c(1, 1, length(pops))))) > 2
2021-06-29 10:37:42 +02:00
# FIXME: 3rd dimension of COUNTS getting dropped
2021-02-01 09:22:58 +01:00
term1 <- squeeze(
2020-10-19 14:08:25 +02:00
sum(
sum(
reshape(
lgamma(
repmat(adjprior, c(1, 1, length(pops))) +
2021-02-01 10:01:02 +01:00
COUNTS[seq_len(nrow(adjprior)), seq_len(ncol(adjprior)), pops, drop=!isarray]
2020-10-19 14:08:25 +02:00
),
c(x, y, z)
),
1
),
2
)
2021-02-01 09:22:58 +01:00
)
2021-02-15 09:49:21 +01:00
if (is.null(priorTerm)) priorTerm <- 0
2021-02-01 09:22:58 +01:00
popLogml <- term1 - sum(lgamma(1 + SUMCOUNTS[pops, ]), 2) - priorTerm
2020-10-19 14:08:25 +02:00
return(popLogml)
}