An Introduction to the Private Language

In this post, I am going to provide a tutorial style introduction to the Private programming language. Private is a privacy preserving probabilistic scripting language. This will be the first of a series of three tutorials that cover different aspects of the language. I’ll link these in here when they are complete.

The tutorial assumes you have a basic understanding of Bayes theorem and Bayesian data analysis. If you don’t, it would be good to take a look at some of these resources before attempting the tutorial:

• Intro to Conditional Probability https://www.youtube.com/watch?v=ibINrxJLvlM Duration: 6:14
• You know I’m all about that Bayes: Crash Course Statistics #24 https://www.youtube.com/watch?v=9TDjifpGj-k Duration: 12:02
• Bayes’ Theorem – The Simplest Case https://www.youtube.com/watch?v=XQoLVl31ZfQ Duration: 5:30

If you want to try the examples out for yourself, you can access the Private language at https://private.mall-lab.com.

The interface has three main areas. You enter your Private code in the editor window and execute it by pressing the Run button. There is a list of what keys do what at https://github.com/ajaxorg/ace/wiki/Default-Keyboard-Shortcuts.

Unlike other programming languages, Private allows you to update your code at any time – you don’t have to interrupt execution. If your latest code does not change the value of a variable it will not be recomputed. If, however, a variable depends on a variable that you change it will automatically be recalculated.

On the right hand side, you will see a list of variables and their values. Some of these variables are the ones you define, while others exist before you start coding. You can use the search box at the top of the list to restrict which variables you see and you can click on individual variables to see them in more detail. Click on the Update Variables button to see the latest values.

The line above the editor window allows you to keep track of how computation is progressing. Private distributes jobs over a cluster of computers in order to speed up execution. On the left hand side, you can see the total number of jobs that are currently executing. These can be broken into the compute jobs which are calculating the values of deterministic variables, the sampler jobs which are computing the values of probabilistic variables and manifold privacy jobs which are evaluating variables to test if they are too sensitive to release. To update these values, click the Show Jobs button.

Estimating the bias of a coin

We are going to begin our exploration of the Private language by considering how one would go about estimating the degree of bias in a coin. If a coin is fair, it has a 0.5 probability of landing on heads and a 0.5 probability of landing on tails. However, it is also possible that we have been handed a ‘trick coin’ which is weighted to favour either heads or tails. In this exercise, we are going to be estimating the probability (or rate) of our weighted coin landing on heads.

We’ll start by entering our data. Let’s assume that we flipped our coin 10 times and saw the following pattern:

Head, Tail, Head, Tail, Tail, Tail, Tail, Head, Tail, Tail

We are going to code the data using a one for the heads and a zero for the tails. To enter the data into Private we type and click Run:

flips = [1, 0, 1, 0, 0, 0, 0, 1, 0, 0]

So, we have the data (3 heads and 7 tails) that we are going to model. Now we want to define the Bayesian model that we will use to estimate the rate parameter. 

To define a Bayesian model, we need to specify the likelihood of the data given the parameters and the priors of the parameters. We are going to assume that the model that generated the data is a Bernoulli distribution. The Bernoulli distribution describes a random variable that has two possible outcomes given a parameter that determines the probability of each. So, we can go ahead and define this as follows:

flips ~ Bernoulli(r)

where r is the unknown rate variable. Note we have used the “~” character here rather than the “=” because we are providing the probabilistic definition of the variable now – rather than setting its value. We are saying “we believe flips is a Bernoulli random variable”.

Before we can estimate r, we also need to specify the prior. We know that r is a probability, so it has to lie between 0 and 1, but we don’t know anything else about it, so we will define it as a uniform distribution as follows:

r ~ Uniform(0,1)

Notice that the uniform distribution takes the upper and lower bounds as input parameters, in this case 0 and 1.

If you look down the variable list you will see: 

flips    [1, 0, 1, 0, 0, 0, ...]
r        Computing

Note r is now ‘Computing’ because all definitions have been provided, but Private is still generating posterior samples of the random variable r. This means that the computer is generating many estimates of the true rate, given what it knows from the data, the likelihood, and the prior. There is some randomness in the generation of these estimates, but more likely estimates will be generated more often.

While we wait for the computer to finish, we can define some other variables that will be useful. Let’s compute the mean and the standard deviation of r as follows:

meanr = mean(r)
stdr = std(r)

Now click Update Variables. If Private has finished estimating r you’ll see the following:

flips   [1, 0, 1, 0, 0, 0, …]
meanr   0.325    
stdr    0.132    
flips   [1, 0, 1, 0, 0, 0, …]    
r       [0.339 … 0.446]        

Inside the brackets next to r are all of the samples (i.e. estimates) that have been generated. It is truncated due to space requirements, so only the first and last values are displayed. If you want to see all of the samples click on the name of the variable (r):

[0.339 0.353 0.302 0.219 0.171 0.456 0.256 0.38 0.407 0.294 0.18 0.41
0.335 0.119 0.163 0.533 0.345 0.432 0.113 0.249 0.653 0.417 0.621 0.077
0.359 0.256 0.439 0.268 0.095 0.459 0.381 0.365 0.321 0.231 0.325 0.447
0.146 0.355 0.277 0.349 0.402 0.492 0.402 0.328 0.488 0.282 0.462 0.127
0.392 0.1 0.522 0.339 0.329 0.282 0.324 0.205 0.388 0.394 0.245 0.164
0.263 0.474 0.451 0.347 0.218 0.203 0.261 0.472 0.294 0.264 0.495 0.178
0.419 0.576 0.487 0.534 0.516 0.665 0.47 0.219 0.494 0.252 0.259 0.519
0.511 0.082 0.604 0.282 0.381 0.223 0.274 0.313 0.159 0.096 0.18 0.326
0.132 0.242 0.445 0.261 0.367 0.359 0.278 0.235 0.585 0.289 0.126 0.192
0.373 0.269 0.316 0.421 0.206 0.373 0.146 0.458 0.312 0.47 0.178 0.16
0.511 0.192 0.161 0.389 0.477 0.347 0.181 0.265 0.353 0.373 0.419 0.306
0.161 0.27 0.213 0.294 0.276 0.316 0.571 0.33 0.31 0.235 0.438 0.472
0.29 0.288 0.244 0.515 0.333 0.394 0.454 0.481 0.189 0.186 0.232 0.288
0.387 0.273 0.117 0.325 0.302 0.126 0.652 0.206 0.104 0.325 0.262 0.482
0.249 0.366 0.271 0.178 0.366 0.238 0.339 0.275 0.114 0.349 0.063 0.276
0.195 0.431 0.528 0.193 0.377 0.357 0.175 0.381 0.349 0.254 0.393 0.125
0.276 0.255 0.116 0.408 0.407 0.23 0.463 0.412 0.406 0.342 0.271 0.64
0.264 0.22 0.528 0.424 0.6 0.232 0.315 0.374 0.47 0.324 0.488 0.329
0.282 0.319 0.351 0.525 0.549 0.327 0.549 0.38 0.569 0.287 0.247 0.431
0.478 0.468 0.535 0.563 0.314 0.523 0.552 0.156 0.164 0.377 0.216 0.394
0.278 0.203 0.188 0.225 0.208 0.483 0.22 0.387 0.121 0.431 0.128 0.213
0.447 0.409 0.118 0.447 0.38 0.41 0.286 0.166 0.334 0.492 0.308 0.36
0.403 0.301 0.445 0.289 0.2 0.259 0.405 0.272 0.367 0.2 0.204 0.384
0.212 0.456 0.445 0.249 0.732 0.281 0.256 0.393 0.224 0.247 0.25 0.219
0.109 0.393 0.356 0.377 0.355 0.325 0.44 0.681 0.393 0.245 0.348 0.346
0.172 0.423 0.271 0.218 0.337 0.474 0.211 0.145 0.597 0.34 0.197 0.153
0.212 0.384 0.38 0.331 0.353 0.324 0.395 0.256 0.093 0.1 0.169 0.264
0.076 0.09 0.107 0.292 0.574 0.699 0.299 0.104 0.14 0.206 0.416 0.32
0.217 0.409 0.333 0.464 0.341 0.158 0.129 0.3 0.329 0.208 0.58 0.263
0.297 0.155 0.239 0.289 0.309 0.371 0.64 0.314 0.33 0.093 0.458 0.31
0.286 0.323 0.137 0.56 0.258 0.399 0.391 0.203 0.131 0.184 0.351 0.33
0.288 0.366 0.322 0.365 0.358 0.611 0.38 0.36 0.471 0.106 0.41 0.532
0.485 0.405 0.386 0.347 0.268 0.202 0.229 0.422 0.292 0.177 0.299 0.283
0.512 0.266 0.381 0.201 0.225 0.475 0.385 0.308 0.587 0.371 0.49 0.402
0.375 0.474 0.492 0.538 0.376 0.145 0.475 0.269 0.223 0.31 0.528 0.273
0.117 0.249 0.422 0.435 0.722 0.212 0.304 0.345 0.225 0.24 0.128 0.251
0.249 0.421 0.1 0.384 0.311 0.426 0.366 0.21 0.186 0.24 0.165 0.198
0.323 0.41 0.477 0.442 0.435 0.279 0.233 0.233 0.408 0.45 0.437 0.255
0.343 0.258 0.227 0.537 0.224 0.477 0.27 0.06 0.198 0.447 0.432 0.141
0.107 0.207 0.081 0.312 0.232 0.287 0.206 0.306 0.178 0.34 0.045 0.646
0.732 0.183 0.275 0.228 0.343 0.359 0.32 0.442 0.275 0.268 0.484 0.406
0.508 0.23 0.29 0.219 0.269 0.277 0.225 0.216 0.129 0.595 0.459 0.276
0.27 0.535 0.097 0.356 0.277 0.299 0.326 0.485 0.312 0.244 0.122 0.361
0.143 0.588 0.334 0.178 0.368 0.071 0.335 0.348 0.064 0.344 0.553 0.432
0.204 0.204 0.259 0.384 0.39 0.298 0.4 0.044 0.422 0.409 0.376 0.473
0.086 0.169 0.139 0.435 0.693 0.433 0.04 0.272 0.475 0.609 0.233 0.203
0.201 0.519 0.619 0.328 0.256 0.426 0.33 0.311 0.154 0.42 0.252 0.365
0.282 0.305 0.222 0.383 0.431 0.48 0.272 0.229 0.259 0.486 0.485 0.182
0.204 0.187 0.592 0.455 0.647 0.083 0.424 0.281 0.216 0.38 0.479 0.45
0.162 0.408 0.455 0.363 0.553 0.133 0.276 0.36 0.208 0.28 0.263 0.405
0.391 0.138 0.408 0.197 0.399 0.349 0.4 0.319 0.315 0.406 0.394 0.272
0.258 0.508 0.193 0.294 0.37 0.21 0.17 0.291 0.321 0.208 0.301 0.505
0.403 0.542 0.344 0.303 0.191 0.351 0.328 0.177 0.42 0.251 0.216 0.214
0.323 0.549 0.224 0.069 0.42 0.209 0.449 0.524 0.455 0.403 0.421 0.228
0.182 0.235 0.309 0.419 0.24 0.592 0.261 0.398 0.309 0.448 0.282 0.656
0.181 0.096 0.23 0.64 0.37 0.451 0.353 0.339 0.253 0.193 0.424 0.199
0.444 0.213 0.28 0.267 0.258 0.317 0.507 0.331 0.457 0.351 0.366 0.68
0.53 0.469 0.318 0.158 0.137 0.409 0.506 0.285 0.24 0.281 0.21 0.38
0.647 0.382 0.414 0.155 0.106 0.391 0.391 0.545 0.414 0.266 0.349 0.512
0.313 0.268 0.602 0.276 0.652 0.33 0.486 0.138 0.538 0.172 0.461 0.131
0.194 0.449 0.586 0.265 0.108 0.18 0.098 0.193 0.244 0.238 0.184 0.303
0.31 0.27 0.47 0.223 0.3 0.312 0.479 0.28 0.32 0.503 0.411 0.462
0.216 0.44 0.45 0.374 0.326 0.283 0.283 0.198 0.407 0.193 0.185 0.512
0.285 0.405 0.412 0.311 0.619 0.207 0.267 0.321 0.223 0.112 0.5 0.178
0.521 0.119 0.596 0.16 0.38 0.307 0.258 0.35 0.336 0.308 0.493 0.364
0.368 0.166 0.283 0.12 0.35 0.253 0.242 0.281 0.4 0.32 0.544 0.521
0.252 0.305 0.272 0.127 0.362 0.447 0.181 0.229 0.317 0.186 0.29 0.502
0.397 0.325 0.305 0.344 0.27 0.259 0.594 0.37 0.395 0.251 0.344 0.316
0.236 0.27 0.576 0.595 0.244 0.481 0.161 0.374 0.268 0.602 0.176 0.127
0.281 0.641 0.281 0.283 0.344 0.624 0.284 0.221 0.436 0.329 0.267 0.397
0.32 0.129 0.393 0.23 0.242 0.697 0.108 0.37 0.346 0.406 0.12 0.192
0.464 0.528 0.318 0.436 0.318 0.56 0.155 0.204 0.388 0.227 0.171 0.373
0.31 0.484 0.489 0.138 0.355 0.3 0.229 0.363 0.472 0.417 0.314 0.314
0.26 0.421 0.38 0.36 0.537 0.335 0.219 0.311 0.38 0.231 0.381 0.589
0.508 0.079 0.141 0.558 0.115 0.477 0.341 0.186 0.37 0.341 0.343 0.335
0.108 0.155 0.23 0.455 0.333 0.368 0.248 0.239 0.35 0.314 0.456 0.356
0.253 0.458 0.431 0.099 0.231 0.154 0.335 0.566 0.413 0.225 0.301 0.703
0.3 0.407 0.667 0.528 0.627 0.061 0.256 0.446 0.386 0.705 0.365 0.198
0.165 0.391 0.37 0.148 0.391 0.228 0.34 0.188 0.512 0.285 0.451 0.189
0.348 0.329 0.26 0.302 0.458 0.139 0.165 0.335 0.288 0.273 0.402 0.403
0.404 0.394 0.329 0.447 0.358 0.193 0.365 0.253 0.541 0.165 0.369 0.11
0.352 0.601 0.145 0.392 0.149 0.355 0.348 0.388 0.194 0.175 0.346 0.13
0.498 0.327 0.515 0.287 0.323 0.38 0.356 0.349 0.421 0.159 0.455 0.18
0.307 0.426 0.485 0.283 0.251 0.337 0.275 0.564 0.271 0.275 0.411 0.449
0.176 0.118 0.527 0.422 0.37 0.456 0.175 0.425 0.431 0.256 0.438 0.322
0.258 0.233 0.277 0.184 0.181 0.559 0.334 0.244 0.214 0.228 0.213 0.335
0.186 0.524 0.25 0.531 0.113 0.391 0.26 0.331 0.24 0.222 0.486 0.612
0.278 0.34 0.232 0.128 0.322 0.402 0.235 0.209 0.274 0.328 0.228 0.288
0.506 0.126 0.213 0.234 0.539 0.383 0.347 0.457 0.404 0.309 0.301 0.362
0.434 0.278 0.43 0.158 0.069 0.353 0.356 0.381 0.313 0.394 0.116 0.309
0.379 0.404 0.155 0.235 0.233 0.238 0.273 0.406 0.232 0.472 0.31 0.261
0.482 0.539 0.214 0.488 0.204 0.2 0.333 0.642 0.399 0.438 0.223 0.253
0.272 0.447 0.451 0.396 0.27 0.237 0.3 0.371 0.231 0.536 0.204 0.149
0.287 0.472 0.287 0.482 0.125 0.523 0.239 0.539 0.372 0.15 0.341 0.175
0.287 0.303 0.408 0.384 0.49 0.043 0.196 0.626 0.286 0.521 0.306 0.311
0.329 0.168 0.444 0.489 0.135 0.281 0.36 0.378 0.078 0.382 0.154 0.403
0.509 0.362 0.37 0.262 0.171 0.297 0.507 0.449 0.225 0.07 0.604 0.206
0.105 0.672 0.383 0.472 0.271 0.315 0.402 0.2 0.291 0.264 0.554 0.463
0.367 0.168 0.639 0.355 0.271 0.609 0.398 0.304 0.28 0.277 0.064 0.275
0.488 0.163 0.312 0.356 0.213 0.392 0.537 0.346 0.479 0.282 0.413 0.385
0.314 0.218 0.119 0.377 0.258 0.189 0.371 0.134 0.365 0.359 0.34 0.436
0.253 0.454 0.265 0.427 0.191 0.367 0.281 0.401 0.068 0.219 0.211 0.112
0.396 0.578 0.394 0.198 0.288 0.58 0.242 0.545 0.252 0.453 0.173 0.236
0.181 0.326 0.438 0.451 0.118 0.358 0.041 0.109 0.188 0.254 0.288 0.207
0.406 0.14 0.477 0.384 0.05 0.154 0.251 0.479 0.234 0.202 0.27 0.318
0.224 0.455 0.338 0.101 0.328 0.4 0.284 0.335 0.373 0.219 0.247 0.201
0.456 0.214 0.325 0.434 0.453 0.438 0.274 0.552 0.115 0.312 0.189 0.364
0.25 0.204 0.109 0.305 0.398 0.281 0.221 0.215 0.269 0.219 0.201 0.246
0.726 0.251 0.318 0.417 0.333 0.437 0.375 0.296 0.127 0.118 0.398 0.126
0.443 0.189 0.488 0.259 0.365 0.399 0.267 0.368 0.352 0.338 0.155 0.235
0.462 0.138 0.203 0.445 0.383 0.143 0.606 0.456 0.175 0.211 0.437 0.452
0.315 0.181 0.17 0.196 0.207 0.251 0.142 0.472 0.212 0.298 0.256 0.168
0.108 0.504 0.425 0.332 0.286 0.356 0.493 0.304 0.272 0.162 0.244 0.311
0.416 0.364 0.335 0.308 0.403 0.464 0.191 0.346 0.321 0.305 0.105 0.344
0.2 0.447 0.205 0.257 0.091 0.098 0.456 0.38 0.428 0.413 0.243 0.287
0.594 0.425 0.351 0.337 0.577 0.27 0.365 0.321 0.417 0.041 0.298 0.286
0.279 0.274 0.062 0.264 0.134 0.349 0.226 0.396 0.599 0.36 0.309 0.208
0.308 0.182 0.249 0.323 0.365 0.459 0.459 0.386 0.314 0.22 0.162 0.225
0.45 0.351 0.287 0.248 0.595 0.213 0.475 0.462 0.199 0.22 0.13 0.292
0.346 0.417 0.273 0.4 0.388 0.282 0.646 0.566 0.215 0.363 0.347 0.242
0.252 0.217 0.284 0.351 0.34 0.687 0.169 0.543 0.299 0.311 0.48 0.141
0.399 0.512 0.157 0.442 0.154 0.39 0.321 0.375 0.072 0.42 0.349 0.234
0.459 0.096 0.319 0.423 0.213 0.48 0.469 0.09 0.646 0.12 0.188 0.129
0.36 0.408 0.387 0.497 0.305 0.344 0.341 0.407 0.316 0.279 0.356 0.376
0.254 0.525 0.216 0.619 0.189 0.281 0.321 0.121 0.33 0.377 0.196 0.451
0.359 0.187 0.463 0.459 0.336 0.261 0.4 0.224 0.23 0.43 0.52 0.336
0.316 0.149 0.402 0.334 0.308 0.416 0.331 0.237 0.272 0.28 0.521 0.301
0.186 0.396 0.217 0.519 0.292 0.32 0.294 0.405 0.33 0.04 0.278 0.556
0.151 0.082 0.125 0.485 0.283 0.42 0.362 0.433 0.138 0.518 0.266 0.4
0.313 0.316 0.39 0.31 0.288 0.43 0.216 0.385 0.104 0.282 0.36 0.254
0.417 0.555 0.33 0.358 0.122 0.242 0.346 0.207 0.423 0.365 0.208 0.138
0.386 0.418 0.377 0.451 0.253 0.211 0.295 0.357 0.476 0.326 0.584 0.474
0.316 0.17 0.134 0.278 0.25 0.31 0.232 0.225 0.648 0.482 0.394 0.215
0.253 0.411 0.117 0.461 0.328 0.095 0.472 0.156 0.248 0.331 0.432 0.071
0.654 0.119 0.363 0.317 0.256 0.383 0.469 0.153 0.611 0.361 0.475 0.366
0.415 0.477 0.384 0.161 0.555 0.452 0.359 0.392 0.186 0.253 0.191 0.521
0.238 0.282 0.127 0.383 0.349 0.218 0.096 0.365 0.491 0.213 0.176 0.288
0.478 0.205 0.469 0.32 0.518 0.328 0.32 0.35 0.458 0.149 0.243 0.409
0.48 0.111 0.107 0.281 0.168 0.378 0.168 0.158 0.238 0.268 0.362 0.157
0.127 0.564 0.356 0.494 0.15 0.056 0.439 0.317 0.439 0.262 0.301 0.343
0.165 0.269 0.479 0.301 0.298 0.252 0.311 0.198 0.27 0.303 0.447 0.287
0.49 0.25 0.37 0.311 0.495 0.361 0.344 0.268 0.193 0.331 0.235 0.073
0.276 0.257 0.162 0.142 0.19 0.464 0.225 0.454 0.409 0.397 0.229 0.311
0.316 0.71 0.218 0.115 0.192 0.328 0.457 0.261 0.451 0.383 0.331 0.481
0.286 0.447 0.267 0.142 0.145 0.275 0.292 0.173 0.44 0.545 0.194 0.228
0.388 0.182 0.378 0.187 0.251 0.441 0.413 0.235 0.318 0.491 0.483 0.146
0.18 0.435 0.572 0.584 0.128 0.447 0.42 0.33 0.343 0.039 0.04 0.356
0.427 0.475 0.188 0.062 0.156 0.443 0.488 0.254 0.223 0.232 0.143 0.3
0.489 0.467 0.314 0.321 0.276 0.235 0.398 0.521 0.224 0.382 0.129 0.34
0.512 0.715 0.276 0.469 0.442 0.393 0.483 0.277 0.38 0.306 0.232 0.376
0.183 0.358 0.566 0.196 0.33 0.149 0.349 0.466 0.328 0.267 0.321 0.602
0.225 0.358 0.169 0.675 0.333 0.383 0.21 0.24 0.29 0.116 0.254 0.259
0.233 0.418 0.39 0.55 0.44 0.145 0.64 0.146 0.18 0.279 0.365 0.313
0.264 0.395 0.317 0.225 0.221 0.337 0.118 0.51 0.451 0.226 0.502 0.346
0.506 0.517 0.514 0.274 0.314 0.331 0.307 0.181 0.503 0.31 0.196 0.326
0.46 0.216 0.169 0.183 0.371 0.134 0.288 0.318 0.4 0.084 0.187 0.275
0.268 0.18 0.343 0.351 0.421 0.523 0.463 0.27 0.253 0.196 0.405 0.274
0.278 0.282 0.235 0.278 0.235 0.347 0.314 0.344 0.22 0.375 0.234 0.393
0.278 0.47 0.311 0.477 0.448 0.194 0.513 0.215 0.309 0.184 0.475 0.374
0.378 0.572 0.436 0.298 0.405 0.371 0.544 0.441 0.169 0.189 0.407 0.516
0.235 0.288 0.399 0.27 0.465 0.384 0.381 0.425 0.23 0.341 0.407 0.263
0.218 0.287 0.431 0.139 0.374 0.076 0.145 0.608 0.388 0.184 0.476 0.504
0.191 0.475 0.33 0.488 0.399 0.548 0.417 0.595 0.196 0.272 0.136 0.491
0.514 0.348 0.11 0.355 0.49 0.136 0.149 0.352 0.407 0.297 0.58 0.167
0.272 0.605 0.265 0.593 0.372 0.292 0.365 0.275 0.136 0.373 0.49 0.505
0.562 0.646 0.307 0.126 0.208 0.49 0.296 0.372 0.443 0.042 0.379 0.168
0.453 0.262 0.142 0.407 0.558 0.275 0.32 0.446]

Because the samples of r are now available, meanr and stdr have been calculated. Note meanr is 0.325, which makes sense as there were 3 heads and 7 tails in our data, so we would expect the rate to be around 0.3.

To make it easier to try out different data sets, we are going to generate fake data with a known probability of the coin landing on heads, so that we can see how well we are able to recover the correct value. To generate 100 flips with a coin that has a probability of 0.3 remove the existing deterministic definition of flips and type:

realr = 0.3
flips = Bernoulli(realr, 100)

so your code is now:

realr = 0.3
flips = Bernoulli(realr, 100)
flips ~ Bernoulli(r)
r ~ Uniform(0,1)
meanr = mean(r)
stdr = std(r)

We say the variable “flips” contains 100 samples of a Bernoulli variable with a rate parameter of 0.3. 

To view the 100 samples that were created, click on the variable name – in this case flips:

[0 0 0 1 0 1 0 0 0 0 0 1 1 0 0 0 1 0 1 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1
0 0 1 0 0 1 1 1 1 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 1 1 0 0 0 0 1 0 1 0 0 0 0
0 0 0 1 0 0 1 1 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0]

Now if you click Update Variables you’ll see:

flips   [0.000 … 0.000]    
meanr   0.264    
stdr    0.044    
realr   0.300    
flips   [0.000 … 0.000]    
r       [0.274 … 0.245]    

Note that meanr is close to the realr value that we used to generate the data.

One of the nice features of Private is that if we change the value of a variable then the other variables change accordingly. For instance, if we now type:

realr = 0.9

And then click Update Variables, we will get:

flips   [1.000 … 1.000]    
meanr   0.892    
stdr    0.032    
realr   0.900    
flips   [1.000 … 1.000]    
r       [0.906 … 0.874]  

So again meanr is a reasonable approximation to the rate we used to generate the data.

You can now generate a histogram of the samples of r like this:

rplot = distplot(r)

If you show the variables using Update Variables you will see:

rplot   [PNG Image]    

Notice that rplot appears as a variable. Click on rplot and you will also see the histogram of the values:

Most of the values of r are around 0.9. The blue line is called the kernel density estimator (kde) and is a smooth representation of the histogram. 

If you change the data again

realr = 0.5

The values will change and the plot will update automatically. Click on rplot now and you will see:

Exercise 1: Look at the standard deviations of the rate variable when the real rate was 0.3, 0.5 and 0.9. Notice that it is greatest for 0.5 and smallest for 0.9. Why is this the case? (Solutions to exercises are provided at the end of the tutorial).

We have now generated code for a model that estimates the bias of the coin. We have also shown that the estimate of coin bias will change appropriately if we alter the data that we give to the model we have created.

Estimating the bias of multiple coins

Now that we have mastered the estimation of the rate of a single coin, let’s spread our wings and consider the case of estimating the rate of multiple identical coins when we flip them all. We can use the binomial distribution to describe the probability of a given number of heads if we have N coins and the probability of each coin being a head is r. 

Again we are going to cheat (I mean facilitate our understanding) by using fake data for which we know the correct rate parameter. First remove all of the code we used in the first example. Now we can generate binomial values as follows:

NumberOfHeads = Binomial(10, 0.3, 100)

which means flip 10 coins, weighted r=0.3 towards tails, 100 times. Now you will see:

NumberOfHeads  [1.000 … 4.000]

To see all of the values click on NumberOfHeads:

[1 4 3 2 1 3 4 3 3 2 2 1 4 3 5 3 3 2 2 3 2 6 5 2 3 4 6 4 2 3 1 1 3 2 4 5 2
1 4 2 5 3 0 4 2 1 5 4 4 3 4 3 5 3 5 2 1 3 1 2 3 4 3 3 1 5 3 5 2 4 3 2 4 3
1 5 3 4 4 5 2 3 3 3 2 4 2 3 3 1 4 3 4 2 4 4 5 5 5 4]

We can generate the histogram of the samples using distplot:

hplot = distplot(NumberOfHeads, bins=[0,1,2,3,4,5,6,7,8,9], kde=False)

Note how we used “kde=False” to stop distplot from adding the smoothing line. (kde is an optional parameter within the distplot function, and if it is not present, it is assumed that kde is true, and a smoothing line will be automatically created). We have also used the bins parameter to indicate the values on the x axis of the plot.

Most of the values are around 3 with a spread from 0 to 7. 

As in the case of one coin, we now need to define the likelihood and the prior. We know that the data was drawn from a binomial distribution, so we can define the likelihood as:

NumberOfHeads ~ Binomial(10, r)

where 10 is the number of coins we are flipping and r is the probability that each of them will land on heads (r being the parameter we are trying to estimate).

Again, we will assume apriori that the rate is drawn from a uniform distribution:

r ~ Uniform(0,1)

Click Update Variables and you should see:

NumberOfHeads     [1.000 … 4.000]    
hplot             [PNG Image]    
NumberOfHeads     [1.000 … 4.000]    
r                 [0.323 … 0.294]   

To see the mean we can type:

meanr = mean(r)

which returns .309, so again our model has generated a pretty good estimate of the rate parameter 0.3. (Reminder, we know the actual rate parameter is 0.3 because we generated samples based on a rate of 0.3. When performing this kind of analysis in a real life setting, we typically wouldn’t know the true rate that we are trying to estimate.)

Exercise 2: Repeat the Bernoulli experiment but with 1000 data points instead of 100? What happens to the mean r? What happens to the standard deviation of r? Why?

Exercise 3: Change the number of coins in the binomial example to 100. What happens to the mean rate? What happens to the standard deviation of the rate? Why?

In the last two sections, we have investigated the use of the Bernoulli and Binomial distributions. These distributions describe processes that have two outcomes, typically True and False, and are useful for modelling quantities like button responses (e.g. Did the participant press “Happy” during this period?) or the occurrence of events (e.g. Did it rain during this period?). In the next section, we will model continuous quantities using the Normal distribution. 

Estimating the mean and standard deviation of male heights

The average male over the age of 18 is about 178 centimeters tall. There is also, of course, quite a bit of variation. The standard deviation of the distribution of male heights is about 10 centimeters. In this section, we will estimate these values from data. Below is a data set of male heights. Clear your work space and copy and paste this dataset into Private:

heights = [176, 181, 164, 176, 176, 181, 180, 191, 152, 182, 169, 188, 182, 182, 190, 180, 189, 169, 190, 194, 173, 179, 155, 183, 191, 186, 178, 174, 179, 182, 188, 169, 186, 169, 182, 173, 171, 172, 174, 169, 179, 170, 174, 184, 201, 180, 183, 184, 196, 180, 168, 194, 174, 161, 189, 174, 186, 168, 176, 193, 177, 176, 187, 170, 175, 181, 168, 179, 159, 168, 186, 182, 162, 177, 176, 187, 172, 178, 178, 197, 175, 170, 170, 181, 186, 207, 187, 175, 163, 174, 169, 173, 170, 169, 164, 178, 201, 178, 198, 164]

You can check the length of the data using the len function:

lenheights = len(heights)

You can also plot the data:

heightsplot = distplot(heights)

Now we need to specify the likelihood and priors of the model we will use. In this case, we are going to use the Normal distribution:

heights ~ Normal(mu, sigma)

mu is the parameter that will represent the mean and sigma is the parameter that will represent the standard deviation. Now we need to specify the priors. We expect heights to be around 150 centimeters, so we specify a Normal distribution with a mean of 150 and a standard deviation that is large (100) to indicate that we don’t want our preconceptions to influence the result too much. The standard deviation of the heights must be positive. Therefore we will use a prior distribution that can only be positive. There are many to choose from, but a good one is the HalfNormal. The HalfNormal is created by taking the absolute value of values from a Normal distribution with mean 0.  

To see what the HalfNormal distribution looks like we can generate 1000 samples from a HalfNormal distribution with a standard deviation of 100, and then plot them:

hn = HalfNormal(100, 10000)
hnplot = distplot(hn, kde=False)

Specify the priors as follows:

mu ~ Normal(150, 100)
sigma ~ HalfNormal(100)

Exercise 4: Now let’s see how we did. Calculate the mean of mu. How close did we get to the real value? Now calculate the mean of sigma. How close did we get to the real value of the standard deviation?

Exercise 5: Generate a plot of the samples of mu and another of the samples of sigma.

You may be thinking to yourself that we had to do a lot of work just to estimate the bias of a coin or the mean of some heights. In the case of the bias, you could just divide the number of heads by the total number of flips. In the case of the heights, you could just add the heights and divide by 100. There are two main reasons that we go to this trouble. Firstly, by thinking about these simple operations as Bayesian models it is easier to see how we can extend them to more complicated models that would not be so easy to estimate. Secondly, simple statistics like the rate or the mean can be used to reverse engineer people’s private data. In the next section, we will demonstrate the first point by showing how to extend our model of heights to a regression model in which we want to determine the impact of a predictor on height. We will leave the second point, until subsequent tutorials.  

How does gender impact height?

Woman tend to be shorter than men. On average a woman over the age of 18 is about 165 centimeters tall. The standard deviation of the distribution of female heights is slightly smaller than for men – around 9 centimeters. In this section, we are going to define a linear regression model that predicts height based on gender.

Start by clearing your code and entering the data:

heights = [176, 181, 164, 176, 176, 181, 180, 191, 152, 182, 169, 188, 182, 182, 190, 180, 189, 169, 190, 194, 173, 179, 155, 183, 191, 186, 178, 174, 179, 182, 188, 169, 186, 169, 182, 173, 171, 172, 174, 169, 179, 170, 174, 184, 201, 180, 183, 184, 196, 180, 168, 194, 174, 161, 189, 174, 186, 168, 176, 193, 177, 176, 187, 170, 175, 181, 168, 179, 159, 168, 186, 182, 162, 177, 176, 187, 172, 178, 178, 197, 175, 170, 170, 181, 186, 207, 187, 175, 163, 174, 169, 173, 170, 169, 164, 178, 201, 178, 198, 164, 161, 152, 163, 168, 166, 155, 166, 185, 161, 173, 161, 164, 184, 158, 138, 145, 160, 164, 173, 176, 182, 174, 178, 163, 153, 156, 165, 166, 158, 173, 147, 150, 171, 155, 167, 167, 161, 183, 156, 154, 162, 165, 169, 165, 171, 151, 176, 169, 179, 167, 165, 169, 162, 177, 169, 163, 154, 170, 199, 153, 169, 154, 156, 169, 162, 172, 153, 178, 168, 156, 169, 146, 162, 163, 167, 175, 174, 174, 161, 160, 172, 164, 184, 167, 166, 173, 159, 172, 153, 159, 152, 163, 143, 157, 154, 155, 172, 164, 154, 168]

The first 100 values are male heights and the second 100 values are female heights. We need to generate a predictor variable to indicate which data are for females and which data are for males. We can do that as follows:

gender = [1] * 100 + [0] * 100

We are using 1 to indicate a male and a 0 to indicate a female. [1] is a list containing a 1 and when we multiply that by 100 we are indicating we want to create a list that has 100 copies of [1]. We do the same for [0] and then add these lists together to concatenate them into one long list. If you click on gender now you should see:

[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

Now we are ready to specify our regression model. We start by specifying heights as a Normal variable as we did before:

heights ~ Normal(mu, sigma)

we’ll define sigma in the same way also:

sigma ~ HalfNormal(100)

To determine mu, we will assume that we have an intercept called muFemale (which will effectively be an estimate of the female height) and then a parameter “diff” that indicates how many centimeters to increase the mean by to capture male heights. Our equation for mu then is: 

mu ~ muFemale + gender * diff

Now we need to add the priors. As far as we know apriori, muFemale is similar to mu for the male heights so we will use the Normal distribution with mean 150 and standard deviation 100. The difference in height between men and women will be a smaller number, but still positive, so we will use a HalfNormal with a standard deviation of 30:

muFemale ~ Normal(150, 100)
diff ~ HalfNormal(30)

When you click Update Variables you should now see:

heights   [176, 181, 164, 176, 176, 181, …]    
gender    [1, 1, 1, 1, 1, 1, …]    
heights   [176, 181, 164, 176, 176, 181, …]    
sigma     [10.157 … 9.131]    
mu        'Not retained.'    
muFemale  [166.320 … 163.131]    
diff      [13.499 … 15.141]

Private has generated samples of muFemale and diff, but has not kept the samples of mu. In many cases, there can be very large numbers of these internally calculated values and so to avoid storing all that data, Private makes a strategic decision to throw them away. 

Exercise 6: Calculate the mean of muFemale. Does it correspond to the real value? 

Exercise 7: Calculate samples of the mean of the male heights? (Hint: You will need to use muFemale and diff). Do they correspond to reality?

Exercise 8: If we wanted to allow for the possibility that females are taller than males, how would we need to change the code?

Summary

That concludes our introduction to the Private language. You should now understand how to construct models of both binary and continuous data and how to calculate statistics of and plot the samples of the parameters. In the next tutorial, we will talk about the estimation of parameters that might depend on private information.

Solutions to Exercises

Exercise 1: Look at the standard deviations of the rate variable when the real rate was 0.3, 0.5 and 0.9. Notice that it is greatest for 0.5 and smallest for 0.9. Why is this the case?

When the true rate is .5, there are many likely estimates on both sides of 0.5 which the model needs to account for. When the rate is .9 however, the distribution will peak at .9, but because there is an upper bound of 1, there is a much smaller pool of possible values greater than .9 that the estimate can take, and while the distribution has a tail to the left, the smaller amount of variance to the right (greater than .9) will decrease the total variance. A rate of .3 is similar to .9, except it is constrained by the lower bound of 0.

Exercise 2: Repeat the Bernoulli experiment but with 1000 data points instead of 100? What happens to the mean r? What happens to the standard deviation of r? Why?

You will notice that the means stay roughly the same, but the standard deviation decreases when more samples are added, as the posterior distribution is becoming more precise. Visually observing the graphs for both distributions, when you have more information (i.e., high n) the posterior becomes more peaked. This means that you are more certain about what values are plausible, and what values are not.

Exercise 3: Change the number of coins in the binomial example to 100. What happens to the mean rate? What happens to the standard deviation of the rate? Why?

The standard deviation is much smaller when 100 coins are flipped. This is because each data point is providing more information about the rate.

Troubleshooting: To update your code, you will need to change the number of coins in both the samples you are generating and the likelihood distribution. If you only update the number of coins in the samples you generate, the parameters for the likelihood will not match the data and you will get an error message.

Exercise 4: Now let’s see how we did. Calculate the mean of mu. How close did we get to the real value? Now calculate the mean of sigma. How close did we get to the real value of the standard deviation?

Add:

meanmu = mean(mu)

and meanmu should be 178.231 which is pretty close to 178 centimetres and:

meansigma = mean(sigma)

which should give 10.270 which is pretty close to 10 centimetres. So we have recovered the values quite well.

Exercise 5: Generate a plot of the samples of mu and another of the samples of sigma.

Type:

muPlot = distplot(mu)
sigmaPlot = distplot(sigma)

Exercise 6: Calculate the mean of muFemale. Does it correspond to the real value? 

meanMuFemale = mean(muFemale)

The mean is 164.447 which is very close to the true value.

Exercise 7: Calculate samples of the mean of the male heights? (Hint: You will need to use muFemale and diff). Does it correspond to reality?

muMale = muFemale + diff
meanMuMale = mean(muMale)

meanMuMale = 178.269, which corresponds quite well (the mean is 178.253 ) .

Exercise 8: If we wanted to allow for the possibility that females are taller than males, how would we need to change the code?

diff would have to be specified as a distribution that can take on negative values, such as the normal distribution.