Skip to content

Instantly share code, notes, and snippets.

@goldingn
Created May 12, 2014 10:44
Show Gist options
  • Save goldingn/7e7e3dfffa4d57c8e114 to your computer and use it in GitHub Desktop.
Save goldingn/7e7e3dfffa4d57c8e114 to your computer and use it in GitHub Desktop.
Playing with probability-of-presence from Maxent-type models
# converting maxent/Poisson model output to true probability of presence
# clear workspace
rm(list = ls())
# set rng seed
set.seed(1)
# number of datapoints
n <- 100000
# 'true' model
f <- function(x) pnorm(-3 * sin(x * 0.6))
# covariate
x <- runif(n, -6, 6)
# 'true' probability
p <- f(x)
# observed presence/absence
y <- rbinom(n, 1, p)
# x where at presence locations: \pi = p(x|y=1)
x_y1 <- x[y == 1]
adj <- 0.6
# Density estimate of p(x|y=1)
f1_star <- density(x_y1,
adjust = adj,
from = -4,
to = 4)
# Density estimate of p(x)
f_star <- density(sample(x, length(x_y1), replace = FALSE),
adjust = adj,
from = -4,
to = 4)
# ~~~~~~~
# Attempts to reproduce p(y=1|x)
# Quantities:
# p(y=1) == prevalence in x == mean(p) ~= mean(y)
py1 <- mean(f(f1_star$x))
# Answer using Bayes theorem:
# p(y=1|x) = (p(x|y=1) * p(y=1)) / (p(x))
plot(f1_star$y / f_star$y * py1 ~ f1_star$x,
col = 'blue',
type = 'l',
ylim = c(-0.3, 1.3))
# Answer using Dave's approximation
# Z = mean(p(x|y=1)) / p(y=1)
# p(y=1|x) = p(x|y=1) / Z
Z <- mean(f1_star$y) / py1
lines(f1_star$y / Z ~ f1_star$x,
col = 'red')
# true answer
lines(f(f1_star$x) ~ f1_star$x,
lty = 2)
# add in lines for 0 and 1
abline(h = c(0, 1),
col = 'grey')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment