232 lines
6.3 KiB
R
232 lines
6.3 KiB
R
|
## Packages
|
||
|
|
|
||
|
|
library(ggplot2)
|
||
|
|
library(ggthemes)
|
||
|
|
|
||
|
|
## Build the prior
|
||
|
|
|
||
|
|
min_95_ci = 1/5000
|
||
|
|
max_95_ci = 1/100
|
||
|
|
|
||
|
|
magic_constant = 1.6448536269514722
|
||
|
|
mean_lognormal = (log(min_95_ci) + log(max_95_ci))/2
|
||
|
|
sigma_lognormal = (log(max_95_ci) - log(min_95_ci))/
|
||
|
|
( 2 * magic_constant)
|
||
|
|
|
||
|
|
num_points = 50000 ## 100000
|
||
|
|
|
||
|
|
xs = c(0:num_points)/num_points
|
||
|
|
ys = dlnorm(xs, meanlog = mean_lognormal, sdlog = sigma_lognormal)
|
||
|
|
ys = ys / sum(ys)
|
||
|
|
ys[1] = 2/3 ## so that it ends as 40% after normalization: 2/3 / (2/3 + 1) = 0.4
|
||
|
|
ys = ys / sum(ys)
|
||
|
|
cat("Checking that sum(ys) = 1")
|
||
|
|
cat(c("sum(ys) =", sum(ys)))
|
||
|
|
|
||
|
|
chicken_v_human = list()
|
||
|
|
chicken_v_human$xs = xs
|
||
|
|
chicken_v_human$ys = ys
|
||
|
|
chicken_v_human$color = c("red", rep("blue", num_points ))
|
||
|
|
chicken_v_human = as.data.frame(chicken_v_human)
|
||
|
|
|
||
|
|
## Plot the prior
|
||
|
|
title_text = "Prior over human vs chicken relative values"
|
||
|
|
subtitle_text=""
|
||
|
|
label_x_axis = "xs"
|
||
|
|
label_y_axis = "P(x)"
|
||
|
|
|
||
|
|
### Plot the prior over the whole x domain
|
||
|
|
ggplot(data=chicken_v_human, aes(x=xs, y=ys, color=color))+
|
||
|
|
geom_point(size=0.5, aes(colour=color), show.legend = FALSE)+
|
||
|
|
labs(
|
||
|
|
title=title_text,
|
||
|
|
subtitle=subtitle_text,
|
||
|
|
x=label_x_axis,
|
||
|
|
y=label_y_axis
|
||
|
|
) +
|
||
|
|
theme_tufte() +
|
||
|
|
theme(
|
||
|
|
legend.title = element_blank(),
|
||
|
|
plot.title = element_text(hjust = 0.5),
|
||
|
|
plot.subtitle = element_text(hjust = 0.5),
|
||
|
|
legend.position="bottom",
|
||
|
|
legend.box="vertical",
|
||
|
|
axis.text.x=element_text(angle=60, hjust=1),
|
||
|
|
plot.background=element_rect(fill = "white",colour = NA)
|
||
|
|
)+
|
||
|
|
scale_colour_manual(values = c("navyblue", "red"))
|
||
|
|
|
||
|
|
getwd() ## Working directory on which the file will be saved. Can be changed with setwd("/your/directory")
|
||
|
|
height = 5
|
||
|
|
width = floor(height*(1+sqrt(5))/2)
|
||
|
|
ggsave("prior.png", width=width, height=height)
|
||
|
|
|
||
|
|
### Plot the prior over only part of the x domain
|
||
|
|
subtitle_text="(zoomed in)"
|
||
|
|
ggplot(data=chicken_v_human, aes(x=xs, y=ys))+
|
||
|
|
geom_point(color="navyblue", size=0.05)+
|
||
|
|
labs(
|
||
|
|
title=title_text,
|
||
|
|
subtitle=subtitle_text,
|
||
|
|
x=label_x_axis,
|
||
|
|
y=label_y_axis
|
||
|
|
) +
|
||
|
|
theme_tufte() +
|
||
|
|
theme(
|
||
|
|
legend.title = element_blank(),
|
||
|
|
plot.title = element_text(hjust = 0.5),
|
||
|
|
plot.subtitle = element_text(hjust = 0.5),
|
||
|
|
legend.position="bottom",
|
||
|
|
legend.box="vertical",
|
||
|
|
axis.text.x=element_text(angle=60, hjust=1),
|
||
|
|
plot.background=element_rect(fill = "white",colour = NA)
|
||
|
|
)+
|
||
|
|
scale_x_continuous(limits = c(0, 0.05))+
|
||
|
|
scale_y_continuous(limits = c(0, 0.005))
|
||
|
|
|
||
|
|
ggsave("prior-zoomed-in.png", width=width, height=height)
|
||
|
|
|
||
|
|
## Construct p(h|x)
|
||
|
|
|
||
|
|
p_w = 0.5
|
||
|
|
rp_estimate = 0.332
|
||
|
|
|
||
|
|
### Construct p(h|xW)
|
||
|
|
|
||
|
|
chicken_v_human$p_h_cond_x_W = rep(1/num_points, num_points+1)
|
||
|
|
|
||
|
|
### Construct p(h|x(not W))
|
||
|
|
are_within_one_order_of_magnitude = function(p1, p2){
|
||
|
|
return(abs(log(p1/p2)/log(10)) <= 1)
|
||
|
|
}
|
||
|
|
are_within_one_order_of_magnitude(1,10)
|
||
|
|
are_within_one_order_of_magnitude(0.1,10)
|
||
|
|
are_within_one_order_of_magnitude(0.1,1)
|
||
|
|
|
||
|
|
start.time <- Sys.time()
|
||
|
|
count = 0
|
||
|
|
cache = TRUE
|
||
|
|
if(!cache){
|
||
|
|
p_h_cond_x_not_W = c()
|
||
|
|
for(x in (0:num_points)/num_points){
|
||
|
|
print(count)
|
||
|
|
count = count + 1
|
||
|
|
is_within_one_oom_from_x = function(y){
|
||
|
|
are_within_one_order_of_magnitude(x,y)
|
||
|
|
}
|
||
|
|
if(is_within_one_oom_from_x(rp_estimate)){
|
||
|
|
num_close_to_x = sum(sapply(xs, is_within_one_oom_from_x))
|
||
|
|
p_h_cond_x_not_W <- c(p_h_cond_x_not_W, 1/num_close_to_x)
|
||
|
|
}else{
|
||
|
|
p_h_cond_x_not_W <- c(p_h_cond_x_not_W, 0)
|
||
|
|
}
|
||
|
|
}
|
||
|
|
end.time <- Sys.time()
|
||
|
|
time.taken <- end.time - start.time
|
||
|
|
time.taken
|
||
|
|
}
|
||
|
|
|
||
|
|
chicken_v_human$p_h_cond_x_not_W = p_h_cond_x_not_W
|
||
|
|
|
||
|
|
### Construct p(h|x) from p(h|xW) and p(h|x(not W))
|
||
|
|
|
||
|
|
chicken_v_human$p_h_cond_x = p_w * chicken_v_human$p_h_cond_x_W + (1-p_w)* chicken_v_human$p_h_cond_x_not_W
|
||
|
|
|
||
|
|
### plot p(h|x)
|
||
|
|
title_text = "P(h|x) update"
|
||
|
|
subtitle_text=""
|
||
|
|
label_x_axis = "xs"
|
||
|
|
label_y_axis = "P(h|x)"
|
||
|
|
ggplot(data=chicken_v_human, aes(x=xs, y=p_h_cond_x, color=color))+
|
||
|
|
geom_point(size=0.05, aes(colour=color), show.legend = FALSE)+
|
||
|
|
labs(
|
||
|
|
title=title_text,
|
||
|
|
subtitle=subtitle_text,
|
||
|
|
x=label_x_axis,
|
||
|
|
y=label_y_axis
|
||
|
|
) +
|
||
|
|
theme_tufte() +
|
||
|
|
theme(
|
||
|
|
legend.title = element_blank(),
|
||
|
|
plot.title = element_text(hjust = 0.5),
|
||
|
|
plot.subtitle = element_text(hjust = 0.5),
|
||
|
|
legend.position="bottom",
|
||
|
|
legend.box="vertical",
|
||
|
|
axis.text.x=element_text(angle=60, hjust=1),
|
||
|
|
plot.background=element_rect(fill = "white",colour = NA)
|
||
|
|
)+
|
||
|
|
scale_colour_manual(values = c("navyblue", "red"))
|
||
|
|
|
||
|
|
ggsave("p_h_x_update.png", width=width, height=height)
|
||
|
|
|
||
|
|
## Calculate p_h
|
||
|
|
|
||
|
|
p_h = sum(chicken_v_human$ys * chicken_v_human$p_h_cond_x)
|
||
|
|
cat(p_h)
|
||
|
|
|
||
|
|
## Calculate p_x_given_h
|
||
|
|
|
||
|
|
chicken_v_human$p_x_given_h = chicken_v_human$ys * chicken_v_human$p_h_cond_x / p_h
|
||
|
|
|
||
|
|
## Plot p_x_given_h
|
||
|
|
|
||
|
|
title_text = "Posterior, P(x|h)"
|
||
|
|
subtitle_text=""
|
||
|
|
label_x_axis = "xs"
|
||
|
|
label_y_axis = "P(x|h)"
|
||
|
|
ggplot(data=chicken_v_human, aes(x=xs, y=p_x_given_h, color=color))+
|
||
|
|
geom_point(size=0.5, aes(colour=color), show.legend = FALSE)+
|
||
|
|
labs(
|
||
|
|
title=title_text,
|
||
|
|
subtitle=subtitle_text,
|
||
|
|
x=label_x_axis,
|
||
|
|
y=label_y_axis
|
||
|
|
) +
|
||
|
|
theme_tufte() +
|
||
|
|
theme(
|
||
|
|
legend.title = element_blank(),
|
||
|
|
plot.title = element_text(hjust = 0.5),
|
||
|
|
plot.subtitle = element_text(hjust = 0.5),
|
||
|
|
legend.position="bottom",
|
||
|
|
legend.box="vertical",
|
||
|
|
axis.text.x=element_text(angle=60, hjust=1),
|
||
|
|
plot.background=element_rect(fill = "white",colour = NA)
|
||
|
|
)+
|
||
|
|
scale_colour_manual(values = c("navyblue", "red"))
|
||
|
|
|
||
|
|
ggsave("posterior.png", width=width, height=height)
|
||
|
|
|
||
|
|
### Plot the posterior over only part of the x domain
|
||
|
|
subtitle_text="(zoomed in)"
|
||
|
|
ggplot(data=chicken_v_human, aes(x=xs, y=p_x_given_h))+
|
||
|
|
geom_point(color="navyblue", size=0.05)+
|
||
|
|
labs(
|
||
|
|
title=title_text,
|
||
|
|
subtitle=subtitle_text,
|
||
|
|
x=label_x_axis,
|
||
|
|
y=label_y_axis
|
||
|
|
) +
|
||
|
|
theme_tufte() +
|
||
|
|
theme(
|
||
|
|
legend.title = element_blank(),
|
||
|
|
plot.title = element_text(hjust = 0.5),
|
||
|
|
plot.subtitle = element_text(hjust = 0.5),
|
||
|
|
legend.position="bottom",
|
||
|
|
legend.box="vertical",
|
||
|
|
axis.text.x=element_text(angle=60, hjust=1),
|
||
|
|
plot.background=element_rect(fill = "white",colour = NA)
|
||
|
|
)+
|
||
|
|
scale_x_continuous(limits = c(0, 0.1))+
|
||
|
|
scale_y_continuous(limits = c(0, 0.005))
|
||
|
|
|
||
|
|
ggsave("posterior-zoomed-in.png", width=width, height=height)
|
||
|
|
|
||
|
|
## Show some indicators
|
||
|
|
|
||
|
|
chicken_v_human$p_x_given_h[1]
|
||
|
|
|
||
|
|
prior_average_value = sum(chicken_v_human$xs * chicken_v_human$ys)
|
||
|
|
posterior_average_value = sum(chicken_v_human$xs * chicken_v_human$p_x_given_h)
|
||
|
|
prior_average_value
|
||
|
|
posterior_average_value
|