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

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16``` ``` 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

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47``` ```# 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 May 30, 2017, 1:36 a.m.