max/ourMELONS
1
0
Fork 0
Code Issues 1 Pull requests Projects Releases Packages Wiki Activity Actions

Retranslated computePopulationLogml() (#24)

Browse source
This commit is contained in:
Waldir Leoncio 2024-04-10 14:19:20 +02:00
parent 98b5e7a154
commit 55d302ab67
1 changed files with 11 additions and 33 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

44
R/computePopulationLogml.R
View file

@ -1,43 +1,21 @@
computePopulationLogml <- function(pops, adjprior, priorTerm) {
# Palauttaa length(pops)*1 taulukon, jossa on laskettu korikohtaiset
# ======================================================== #
# Limiting COUNTS size #
# ======================================================== #
if (!is.null(adjprior)) {
nr <- seq_len(nrow(adjprior))
nc <- seq_len(ncol(adjprior))
COUNTS <- COUNTS[nr, nc, pops, drop = FALSE]
} else {
COUNTS <- NA
}
nr <- seq_len(nrow(adjprior))
nc <- seq_len(ncol(adjprior))
x <- size(COUNTS, 1)
y <- size(COUNTS, 2)
x <- size(baps.globals$COUNTS, 1)
y <- size(baps.globals$COUNTS, 2)
z <- length(pops)
# ======================================================== #
# Computation #
# ======================================================== #
term1 <- NULL
if (!is.null(adjprior)) {
isarray <- length(dim(repmat(adjprior, c(1, 1, length(pops))))) > 2
term1 <- squeeze(
sum(
sum(
reshape(
lgamma(
repmat(adjprior, c(1, 1, length(pops))) + COUNTS[nr, nc, pops, drop = !isarray]
),
c(x, y, z)
),
1
),
2
)
)
}
if (is.null(priorTerm)) priorTerm <- 0
popLogml <- term1 - sum(lgamma(1 + SUMCOUNTS[pops, ]), 2) - priorTerm
return(popLogml)
rep_adj <- repmat(adjprior, c(1, 1, z))
gamma_rep_counts <- matlab2r::gammaln(rep_adj + baps.globals$COUNTS[, , pops])
gamma_sum_counts <- rowSums(matlab2r::gammaln(1 + baps.globals$SUMCOUNTS[pops, , drop = FALSE]))
gamma_rep_counts_sum <- colSums(colSums(reshape(gamma_rep_counts, c(x, y, z))))
gamma_rep_counts_reshaped <- squeeze(gamma_rep_counts_sum)
popLogml <- gamma_rep_counts_reshaped - gamma_sum_counts - priorTerm
return(popLogml[, , drop = FALSE])
}