forked from pierrejacob/couplingsmontecarlo
-
Notifications
You must be signed in to change notification settings - Fork 0
/
2otarrangement.R
64 lines (54 loc) · 1.94 KB
/
2otarrangement.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
rm(list=ls())
library(doParallel)
library(doRNG)
registerDoParallel(cores=10)
library(ggplot2)
theme_set(theme_void())
library(gridExtra)
colors <- c(rgb(1, 0.1, 0.3), rgb(0.3, 0.1, 1))
## two distributions in [0,1]
## density p
densityp <- function(x){
0.25 * dbeta(x, 40,10) +
0.5 * dbeta(x, 20,20) +
0.25 * dbeta(x, 10, 45)
}
## density q
densityq <- function(x) exp(dbeta(x, 4, 4, log=T) + cos(5*pi*x))
constantq <- integrate(densityq, 0, 1)$val
## cdf p
cdfp <- function(x) integrate(densityp, lower=0, upper=x)$val
## cdf q
cdfq <- function(x) integrate(densityq, lower=0, upper=x)$val / constantq
## inverse cdf p
inversecdfp <- function(z){
uniroot(function(x) cdfp(x) - z, interval=c(0, 1))$root
}
## inverse cdf q
inversecdfq <- function(z){
uniroot(function(x) cdfq(x) - z, interval=c(0, 1))$root
}
# nsamples <- 1e4
# unifs <- runif(nsamples)
# xs <- sapply(unifs, inversecdfp)
# ys <- sapply(unifs, inversecdfq)
nsamples <- 1e4
unifs <- runif(nsamples)
xys <- foreach(irep = 1:nsamples) %dorng% c(inversecdfp(unifs[irep]), inversecdfq(unifs[irep]))
xys <- sapply(xys, function(x) x)
xs <- xys[1,]
ys <- xys[2,]
gsegments <- qplot(x=0, xend=1, y=xs, yend=ys, geom="blank") +
geom_segment(alpha=0.05)
gmargleft <- qplot(x=xs, geom="blank") +
geom_histogram(aes(y=..density..), fill=colors[1]) +
stat_function(fun=function(x) densityp(x)) +
coord_flip() +
scale_y_reverse()
gmargright <- qplot(x=ys, geom="blank") +
geom_histogram(aes(y=..density..), fill=colors[2]) +
stat_function(fun=function(x) densityq(x) / constantq) +
coord_flip()
g <- gridExtra::grid.arrange(gmargleft, gsegments, gmargright,
ncol=3, nrow=1, widths=c(1,4,1), heights=4)
# ggsave(filename="otarrangementwithmarginals.pdf", plot=g)