tm: Log-likelihood Estimation using the Time Machine

In TimeMachine: Time Machine

Description

Estimates the log-likelihood for the given parameter values using the Time Machine.

Usage

    tm(transitions, pi=NULL, population, n, mu, samples, threads=NULL)

    ## S3 method for class 'tm'
density(x, ...)

    ## S3 method for class 'tm'
mean(x, ...)

    ## S3 method for class 'tm'
print(x, ...)

    ## S3 method for class 'tm'
quantile(x, ...)

    ## S3 method for class 'tm'
summary(object, ..., digits=max(3, getOption("digits")-3))

Arguments

transitions	Transition matrix between types
pi	Stationary distribution associated with the transition matrix, or NULL to compute it
population	Vector representing the initial population
n	Target population size
mu	Mutation rate
samples	Number of samples to simulate
threads	Number of threads used for simulations, or NULL to use the default OpenMP setting
x,object	An object of class tm as returned by a call to tm
...	Further arguments passed to the corresponding generic function
digits	Minimum number of significant digits for number formatting

Value

An object of class tm representing the result of the simulations.

population	Vector representing the initial population
n	Target population size
mu	Mutation rate
logliks	Vector of simulated log-likelihoods
corrections	Vector of correction terms
coalescent.events	Matrix of column vectors (one per simulation) containing SENs when coalescent events were simulated
final.populations	Matrix of columns vectors (one per simulation) containing final population sizes for each type
simulation.times	Computation time for each simulation (in seconds)
total.time	Total computation time (in seconds)

Examples

# Load example dataset
data(pdm)

transitions <- full.transitions(pdm$unitary.transitions, pdm$loci)
pi <- stationary.dist(transitions)
n <- 10
samples <- 10

# Estimate log-likelihood for a fixed value of the mutation rate
mu <- 1
est.res <- tm(transitions, pi, pdm$population, n, mu, samples)
print(est.res)

# Estimate log-likelihood for different values of the mutation rate
mus <- seq(0.1, 10, 0.1)
estimates <- sapply(mus, function(mu) {
    tm(transitions, pi, pdm$population, n, mu, samples)
}, simplify=FALSE)

# Compute mean log-likelihood for each value of the mutation rate mu
mean.logliks <- do.call(rbind, lapply(estimates, function(x) {
    c(x$mu, mean(x))
}))

# Plot mean log-likelihoods
plot(mean.logliks, pch=20, xlab=expression(mu), ylab="Mean log-likelihood",
     main=paste("Mean log-likelihoods - Sample size:", samples))

# Estimate log-likelihood for different target population sizes and fixed
# mutation rate
ns <- 2:50
mu <- 1
estimates <- sapply(ns, function(n) {
    tm(transitions, pi, pdm$population, n, mu, samples)
}, simplify=FALSE)

# Compute mean correction term/log-likelihood ratios for each target population size
mean.corrections <- do.call(rbind, lapply(estimates, function(x) {
    c(x$n, mean(x$correction / x$logliks))
}))

# Plot mean correction term/log-likelihood ratios
plot(mean.corrections[,1], rev(mean.corrections[,2]), pch=20,
     xlab="Target population size", ylab="Mean correction term/log-likelihood",
     axes=FALSE)
axis(1, at=mean.corrections[,1], labels=rev(mean.corrections[,1]))
axis(2)

TimeMachine documentation built on Sept. 11, 2018, 5:03 p.m.