# binomial differential privacy
obs <- rbinom(100, 1, 0.5)
mean(obs)
#> [1] 0.51
rlaplace <- function(n, sensitivity, epsilon) {
u <- runif(n, -0.5, 0.5)
b <- sensitivity / epsilon
sign(u) * b * log(1 - 2 * abs(u))
}
var(rlaplace(10000, 1, 0.25))
#> [1] 31.70462
dp_mean <- function(obs, sensitivity, epsilon) {
(sum(obs) + rlaplace(1, sensitivity, epsilon))/length(obs)
}
dp_mean(obs, 1, 0.1)
#> [1] 0.377109
N <- 1000
sensitivity <- 1
epsilon <- 0.1
means <- sapply(1:10000, \(x) {
obs <- rbinom(N, 1, 0.5)
c(true_mean = mean(obs), dp_mean = dp_mean(obs, sensitivity, epsilon))
})
quantile(means[1,], c(0.025, 0.975))
#> 2.5% 97.5%
#> 0.469 0.531
quantile(means[2,], c(0.025, 0.975))
#> 2.5% 97.5%
#> 0.4575225 0.5414932
# var(X + Y) = var(X) + var(Y) + 2 cov(X, Y)
# var(X) = p*(1-p) / n
var(means[1,])
#> [1] 0.0002520712
0.5^2 / N
#> [1] 0.00025
var(means[2,])
#> [1] 0.0004457994
(0.5^2 / N + (2*(sensitivity/epsilon)^2 / N) / N)
#> [1] 0.00045
rolling_means <- sapply(1:100, \(i) {
obs <- rbinom(i*100, 1, 0.5)
c(true_mean = mean(obs), dp_mean = dp_mean(obs, sensitivity, epsilon))
})
plot(rolling_means[1,], type = "l")
lines(rolling_means[2,], col = "darkgreen")
Created on 2023-12-17 with reprex v2.0.2