Compute a tree-sequence expression across a given number of replicates
Source:R/tree-sequences.R
ts_replicate.RdThis function can be used to compute a given number of replicates of a specified tree-sequence statistical expression
Examples
init_env()
p1 <- population("p1", time = 4000, N = 1000)
p2 <- population("p2", time = 3000, N = 1000, parent = p1)
p3 <- population("p3", time = 2000, N = 1000, parent = p2)
model <- compile_model(list(p1, p2, p3), generation_time = 1)
schedule <- schedule_sampling(model, times = 0, list(p1, 10), list(p2, 10), list(p3, 10))
samples <- ts_names(model, split = "pop", schedule = schedule)
# this is how we can compute a single tree-sequence statistic with slendr:
# 1. first we simulate a tree sequence from a model
ts <- msprime(model, sequence_length = 1e6, recombination_rate = 1e-8, schedule = schedule)
# 2. then we run a desired tskit-wrapper function
ts_f3(ts, A = samples["p1"], B = samples["p2"], C = samples["p3"], mode = "branch")
#> # A tibble: 1 × 4
#> A B C f3
#> <chr> <chr> <chr> <dbl>
#> 1 p1 p2 p3 3808.
# this is how we can compute the same function across multiple replicates at once
ts_replicate(
n = 10,
{
ts <- msprime(model, sequence_length = 1e6, recombination_rate = 1e-8, schedule = schedule)
ts_f3(ts, A = samples["p1"], B = samples["p2"], C = samples["p3"], mode = "branch")
}
)
#> # A tibble: 10 × 5
#> A B C f3 rep
#> <chr> <chr> <chr> <dbl> <int>
#> 1 p1 p2 p3 5163. 1
#> 2 p1 p2 p3 4241. 2
#> 3 p1 p2 p3 3995. 3
#> 4 p1 p2 p3 4123. 4
#> 5 p1 p2 p3 3949. 5
#> 6 p1 p2 p3 4539. 6
#> 7 p1 p2 p3 3357. 7
#> 8 p1 p2 p3 2981. 8
#> 9 p1 p2 p3 3884. 9
#> 10 p1 p2 p3 3702. 10