Skip to content

Instantly share code, notes, and snippets.

@thomvolker
Created December 17, 2023 11:18
Show Gist options
  • Save thomvolker/13fe1b3efb1988ab4ecb2a9664534f36 to your computer and use it in GitHub Desktop.
Save thomvolker/13fe1b3efb1988ab4ecb2a9664534f36 to your computer and use it in GitHub Desktop.
Differential privacy of a proportion
# 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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment