DNAprofiles (version 0.3.1)

sim.q: Estimate LR exceedance probabilities (with importance sampling)

Description

Estimate LR exceedance probabilities (with importance sampling)

Usage

sim.q(t, dists, N = 1e+05, dists.sample = dists)

Arguments

t

numeric (vector), threshold

dists

list of per-locus probability distributions of a likelihood ratio

N

integer, number of samples

dists.sample

if dists.sample is not equal to dists, then importance sampling is applied. Sampling is done according to dists.sample, while exceedance probabilities are estimated for dists.

Value

• numeric (vector) with estimated probabilities

Details

For a combined likelihood ratio $$LR=LR_1 LR_2 \times LR_m,$$ define $q_{t|H}$ as the probability that the LR exceeds $t$ under hypothesis $H$, i.e.: $$q_{t|H} := P(LR>t|H).$$ The hypothesis $H$ can be $H_p$, $H_d$ or even another hypothesis. The current function estimates $q_{t|H}$ by taking $N$ samples from the distributions specified by the dists parameter and computing the empirical fraction of the product of the samples that exceeds $t$.

Importance sampling can be used by supplying different distributions to sample from and to estimate the exceedance probabilities for. For instance, the exceedance probability for $H_d$ can be estimated by sampling from $H_p$ and an appropriate weighting of the samples. See the paper and examples for details.

Examples

data(freqsNLngm)

# dist of PI for true parent/offspring pairs
hp <- ki.dist(hyp.1="PO",hyp.2="UN",hyp.true="PO",freqs.ki=freqsNLngm)

# dist of PI for unrelated pairs
hd <- ki.dist(hyp.1="PO",hyp.2="UN",hyp.true="UN",freqs.ki=freqsNLngm)

set.seed(100)

# estimate P(PI>1e6) for true PO
sim.q(t=1e6,dists=hp)

# estimate P(PI>1e6) for unrelated pairs
sim.q(t=1e6,dists=hd) # small probability, so no samples exceed t=1e6

# importance sampling can estimate the small probability reliably
# by sampling from H_p and weighting the samples appropriately
sim.q(t=1e6,dists=hd,dists.sample=hp)