/
plot_mixt_density.R
75 lines (70 loc) · 3.8 KB
/
plot_mixt_density.R
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
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#' Plot the density of some given sample(s) with mixture distributions.
#'
#' Plot the density of the sample(s) with optional arguments to improve the visualization.
#'
#' @param samples Observed samples (mixture distributions) from which the density will be plotted.
#' @param user.bounds (default to NULL) Bounds to limit the range of x-axis when plotting.
#' @param support Support of the distributions, to know whether density plot or histogram should be displayed.
#' @param main Title for the plot.
#'
#' @return a plot with the densities of the samples provided as inputs.
#'
#' @examples
#' ##### Continuous support:
#' list.comp <- list(f1 = "norm", g1 = "norm",
#' f2 = "norm", g2 = "norm",
#' f3 = "norm", g3 = "norm")
#' list.param <- list(f1 = list(mean = 5, sd = 1), g1 = list(mean = 2, sd = 0.7),
#' f2 = list(mean = 0, sd = 1), g2 = list(mean = -3, sd = 1.1),
#' f3 = list(mean = 9, sd = 1), g3 = list(mean = 6, sd = 2))
#' ## Simulate data:
#' sim1 <- rsimmix(n = 300, unknownComp_weight = 0.8, comp.dist = list(list.comp$f1,list.comp$g1),
#' comp.param = list(list.param$f1, list.param$g1))$mixt.data
#' sim2 <- rsimmix(n= 250, unknownComp_weight = 0.85, comp.dist = list(list.comp$f2,list.comp$g2),
#' comp.param = list(list.param$f2, list.param$g2))$mixt.data
#' sim3 <- rsimmix(n= 400, unknownComp_weight = 0.6, comp.dist = list(list.comp$f3,list.comp$g3),
#' comp.param = list(list.param$f3, list.param$g3))$mixt.data
#' plot_mixt_density(samples = list(sim1,sim2,sim3), user.bounds = NULL, support = "continuous")
#'
#' ####### Countable discrete support:
#' list.comp <- list(f1 = 'pois', g1 = 'pois',
#' f2 = 'pois', g2 = 'pois')
#' list.param <- list(f1 = list(lambda = 7), g1 = list(lambda = 1),
#' f2 = list(lambda = 2), g2 = list(lambda = 15))
#' sim1 <- rsimmix(n=4000, unknownComp_weight=0.5, comp.dist = list(list.comp$f1,list.comp$g1),
#' comp.param=list(list.param$f1,list.param$g1))$mixt.data
#' sim2 <- rsimmix(n=3500, unknownComp_weight=0.3, comp.dist = list(list.comp$f2,list.comp$g2),
#' comp.param=list(list.param$f2,list.param$g2))$mixt.data
#' plot_mixt_density(samples = list(sim1,sim2), user.bounds = NULL, support = "discrete")
#'
#' @author Xavier Milhaud <xavier.milhaud.research@gmail.com>
#' @export
plot_mixt_density <- function(samples, user.bounds = NULL, support = c("continuous","discrete"), main = "")
{
n <- max(sapply(X = samples, FUN = length))
n_pop <- length(samples)
leg <- paste("Sample ", 1:n_pop, sep="")
#leg <- eval(parse(text = paste("latex2exp::TeX('X_", 1:n_pop, "')", sep="")))
densities <- lapply(X = samples, FUN = stats::density)
x.axe.lim.old <- x.axe.lim <- c(min(apply(sapply(densities, "[[", "x"), 2, min)), max(apply(sapply(densities, "[[", "x"), 2, max)))
if (!is.null(user.bounds)) x.axe.lim <- c(user.bounds[1],user.bounds[2])
old_par <- graphics::par()$new
if (support != "discrete") {
y.axe.lim <- max(apply(sapply(densities, "[[", "y"), 2, max))
for (i in 1:n_pop) {
base::plot(densities[[i]], xlim = x.axe.lim, ylim = c(0,y.axe.lim), xlab = "", main = "", col = i, lty = i)
graphics::par(new = TRUE)
}
} else {
## FIXME: pour l'argument 'breaks' defini comme ci-dessous, cela ne marche que si la loi discrete est a support positif!
#?barplot
for (i in 1:n_pop) {
graphics::hist(samples[[i]], freq = TRUE, breaks = 0:ceiling(x.axe.lim.old[2]), xlim = x.axe.lim, ylim = c(0,n/2.3),
xlab = "", main = "", border = i, lty = i)
graphics::par(new = TRUE)
}
}
graphics::title(main = main)
graphics::legend("topright", legend = leg, lty = 1:n_pop, col = 1:n_pop, bty = "n")
on.exit(graphics::par(new = old_par))
}