Probability & Statistics - Random Sampling

Hello everyone, My name is Eduardo Serna. I'm gonna be giving a talk on random sampling in the Wolfram language. This will be hopefully the first part of a series of talks from our team, the Probability and Statistics group in Wolfram. I'm a consultant there. The other half of my team is Gosa Converska who's hopefully around to, will be helping me answer questions today, if there are any. And we'll be giving some of the talks in this series. So the first talk is going to focus on univariate distributions. And this is the outline of the talk. This is a very, very big subject and I've chosen to. The aim of the talk, I guess is to balance between, balance between it being informative. So I will try and hopefully expose the functionality that is available in the Wolfram language. I'm going to try and make it educational so give you an idea of the different methods and techniques that are used in sampling. And I'm also going to try and make it entertaining. And this is obviously very subjective, but from my point of view, I am a physicist by training. That was very interesting probability when I was a physicist. And this is by far my favorite part of my job here at Waltham. The random number generation is very, very fun for me and there's lots of very interesting math and you know, computing that goes into generating these elusive things which are pseudo random numbers. So the outline of the talk is gonna. I'm gonna first talk about uniform RNGs. I then talk about continuous distributions, introducing the different techniques that are available to us. And this is by no means exhaustive at all. And I've left some very big methods here, but those are for follow up talks. So I'll talk about functional transformations, inversion as an example of it acceptance, rejection. Use these methods to, to introduce how one samples from discrete distributions. And then that's going to be me. So randomness is what it means for something to be random is actually kind of quite a profound concept and it's not entirely obvious and loads of mathematicians and physicists and physicists and philosophers have discussed it. But from a physicist point of view, I think a useful way to think about this is to look that we have some complicated process that we're not privy to all the internal details, the internal states. And we're just looking at some cross section of it where we have some, we would call them observables in physics. And we have, from the point of view of those observables we've averaged out, we integrated out all of the internal state and we just get something that we. That looks random. Essentially it's not intrinsically random. If we had access to everything we, it wouldn't be random. We would be able to see exactly where everything comes from. But because we don't have access to all the internal details, it looks random. And computers are largely deterministic. There's some issues about like cosmic rays and cache races and stuff like that, but the most deterministic and so generating a pseudo random variable from a computer is a not entirely obvious, but B we can use essentially a procedure that one could think is inspired by nature and just, just it involves transformations to just garble the inputs in the internal state enough that we lose track of where we started and we start. The building blocks for producing random number generators are uniform random number generators from a pseudo random distribution and the builder box that we use to build more complicated samplers from non uniform distributions, processes, spatial distributions, etc. And we want them to look as random as possible while at the same time balancing other factors such as crucially speed, but also memory and simplicity of the. And this and other useful things such as reproducibility is a key aspect because with computers we kind of want to have our cake and eat it too. And we also want to be able to reproduce the results by controlling the seed of the what's called the seed of the random number generator. And here we have this Random real is a function that generates reals from a uniform distribution inhibiting serum 1. We compute a million and we plot the histogram and we see that it's essentially uniform with small jitters that are plausible and reassuring. These jitters are what this histogram would look essentially the same if this was a true random number generator and not pseudo random number generator. Okay, so in Wolfram language have a family of functions to sample from different uniform different types or different tasks in this case. So we have random integer or real or complex for numbers, we have random choice and random sample for lists, with or without replacement. We have primes for priming numbers and more esoteric stuff like random color, random image, etc. So let's look at two very important aspects of the framework here, which is to be able to specify the seating. As I mentioned earlier, we want to be able to specify the ceiling because we want to be able to get reproducible results. This is important because sometimes we want to. We've used the random number generator in a computation and we want to go back and see run it back by essentially specifying the same seed and make sure that nothing funny is going or we understand the computation. And this is done in the Wolfram language by using seedrandom. Once we specified a seed, all of the following, the random numbers that are computed following that seed will be the same. So if we do the same thing twice with the same seed, we will get the same random real. Another very important thing that we can do is to leverage the arbitrary position arithmetic system in the OR from language to be able to get arbitrary precision, random arbitrary precision numbers that have more random digits that are allowed by the hardware. And this is again very useful for certain computations where things might be ill conditioned or you're close to singularity and there's no algorithm that will save you from that and you just need more digits than those that are available. And machine precision, which here I generate the random real and the precision is machine precision and that's about 16 digits, decimal digits, but I can generate 24. And then we would be said we would be able to use these to compute at 24 digits. This is of course more expensive, slower, but it's sometimes needed and it's very important to have the option. Another nice aspect of the framework is that you can plug in your own custom uniform random number generator. And here's a reference that is very complete and it talks a lot. It has much more information than I'm going to give in my talk. It talks a lot of what I'm talking about. The references are there and in particular it's very detailed on how to, how to build, how to plug in your own generator. I'm not going to go over the details, all of the details. I just wanted to show an example of one generator you can plug in, which is pcg. It's the permuted congruential generator. It's a relatively recent algorithm that's, that appeared recently. It has very strong statistical properties and among other nice things like look ahead and it's very light and simple. And I wanted to bring it up because it has. Well, first of all, it's state of the art. So it's possibly in many ways the best one there is. We should really, we haven't had time to, but we really should add it to the, to, to the family that we support natively. And it has, it's also very simple and you can see some of the aspects that I was discussing earlier where you just have a state and you, you have a state and you garble it a little bit, do some computation with it that mixes it about and the state is bigger than the thing that you're going to be returning. The result in this case has a bit width that is the state is twice as big as the result. So you do some operation on the state and then you do some further operation that is just bit operations that are mixing things about xoring, shift, writing, rotating, right. And then you return the rest, the result and the new state. And this is how the, this is essentially the heart of the code of this generator. And as you can see, it's very, very simple. This, this code is just to, to make sure that we have the correct bit width. And I'm not gonna, you know, show the rest of the boilerplate that involves in plugging it in. But that's the heart of the generator and we can see that if we just specify the method to seed random. This is how you spec. Once you've, once you've defined the code, you specify this method to seed random and we generate, we can generate 10 to the 4. And it looks entirely reassuringly plausible. There's more jitters now because we only have 10 to the 4. The size of the sample is only 10 to the 4. Okay. Another thing I wanted to point out is that there's loads of very like this in this reference. There's loads of different built in samplers that you can use and in particular there's some that can be useful for very specific topics, tasks like the low discrepancy ones that are ones that have like a lower quality of statistical properties. But the trade off is that there's more spread out the numbers that they produce and that could be useful for certain statistical, for certain numerical tasks like integration. Okay, so let's look at now at non uniform random generators and to sample from, in the Wolfram language, to sample from non uniform random nano generators, we can use random variate, random variate, which works for both continuous and discrete distributions. So here we have the normal distribution, which is the most canonical one, the continuous distribution, also called the Gaussian or the bell curve. We sample a million data points and the histogram looks like it should. This blue line is the PDF and the histogram is bell shaped. We can also sample from discrete distribution like the Poisson distribution. And here we're plotting the, we call it the PDF. It's the same function, but it's actually for discrete distributions. It's technically called the probability mass function because it's only defined at discrete integer points. And we see that the histogram agrees with the, with the Poisson be it PDF There's a question here. Okay. Is that a link exists between number theory and probability theory. I believe in regards to prime distribution. I'm not an expert of number in number theory, but okay, yes, yes, there are, there is some, there's some random measures that you can define on these, on these prime distributions. Yes, I will discuss it. I think I'm going to discuss it later if I have the time. Okay. Okay, let's keep going. Okay, so let's talk about functional transformations. So the next step I guess is to generate samples from a non uniform form distribution using uniform random samples. And to do so we going to have to define the cdf, which is the cumulative distribution function. And it's a function that just describes how much probability there is up to a certain number. So it's the CDF is denoted by capital F of a random variable, capital X capital and it's a function of a number X and it just counts how much probability there is when the random variable is equal or smaller to a certain value. And here I show, here I show. This is the PDF, the again of the Gaussian that I showed earlier. And it has the density, the probability is all in close to the center. And the PDF, the CDF, sorry, is a function that actually interpolates between 0 and 1 from minus infinity to infinity. And it just counts how much probability has accumulated as you go, as you go to the right. Technically the CDF is the integral of the, of the PDF, but it's actually much more general because sometimes the PDF is not defined like for instance for discrete distributions. So it's better to work with CDFs. And this formula allows us to construct random variables by applying a function to other random variable. So if we have a function G and we define use it to define a random variable Y, we can try and figure out what's the distribution of Y. And this can be seen by just applying the definition that G of X is smaller or equal to small Y here. And to express this in terms of the known F of X of the CDF of X, we need G to be invertible. Because if it's invertible we can just apply it to both sides of this inequality and we will get X random variable X on this side and G minus 1 of Y on this side. And this is actually the definition here, as we saw, this is the definition of something of X small X of the CDF of X. So there's two possibilities. For G to be invertible, it has to be either strictly increasing or strictly decreasing. There's two possibilities. If it's strictly increasing, then this inequality doesn't switch sides. And we get that the CDF of y is CDF of DG minus 1, the inverse of this G transformation on Y. If it's decreasing, we use the fact that if to one side of this inequality there's one minus the amount of total probability of the other. So we switch so we just can call this one minus the CDF of g minus one of Y. Okay? And we can use this transform distribution to figure out the distribution of these transform random variables. Transform distribution takes an expression and variables that distributed. This symbol means distributed certain distributions, random variables. And we can it can figure out what the output distribution is. So for instance, if we use a simple invertible function like the minus log, we will find that on a uniform distribution we will find that this has transform distribution knows about this one, which is the exponential distribution. And we can see that the CDF is essentially one minus the inverse of this lock that I used here. And this is a valid way to sample from. It's not the one we use, but this is like fast enough so we can generate 10 to the 7 uniform. This here I'm explicitly using uniform distribution. But this could have been random real. 10 to the 7 is exactly the same thing in this case. And we take the log of that data, flip the sign. This is so fast because these are vectorized operations there, which means the machine operates when on the whole data at a time, not one by one. And we see that the data is indeed distributed like an exponential distribution. This is a PDF of an exponential distribution. Okay, so let's look at what's called the inversion formula. Now. So what we actually want is to generate samples from a distribution D having access to samples from the uniform distribution. And we want to figure out so what's the G that we need to apply or was the transformation we need to apply to get there. So, for a uniform random variable, I hope it's not easy. It's not hard to see that the CDF takes the following form. It's just the identity function. And this can be seen by looking at the PDF, which is just a horizontal line because it's the same everywhere. That's what uniform means. And the integral of that. The CDF is essentially the linear function x with slope 1. So once we have that, once we have access to the CDFs of U and the one we're looking for X, we can use the inverse of X, the cdf, the inverse CDF of X as a Transformation function G to construct variables that are distributed as X. So here I just. This is the formula I used before for transforming functions. And we just replace Z here with U g with f minus 1 of x and we substitute, and we get that the CDF of Y is indeed the CDF of X. So X and Y have the same distribution. This means that by applying the inverse of the CDF of X to a uniform sample, we obtain samples from X, which is what we wanted. Okay, so let's construct samples from the normal distribution using its inverse cdf. The CDF of the normal distribution is called the IRFC function. And it's. It just. It's essentially the integral of its. Of its PDF and the inverse distribution. Its inverse function is the inverse earth. And these are not elementary special functions, elementary being no use very often, I guess. And as such, they're much more, much more slower to compute than, say, the exponential. And this is an issue because although we can use this inverse to again apply to uniform random variates and apply the inverse, we will get samples from the normal distribution. It's actually quite slow. This is because the inverse is. Even though it's vectorized, it's relatively slow. So, as we saw in the last slide, evaluating certain inverse CDFs can be actually quite expensive. However, there are all sorts of clever transformations to sample from important distributions. And I'm gonna look at one, which is called the box molar transform. Now, which is important because it's the one that we use and one of the most canonically used to sample from the normal distribution, which is very important. And also I chose it because it has a very clever trick, in my opinion, that is quite fun. So what it does is it uses the observation that the square root of a sum of the square of two random variables that are normally distributed has a CDF that is easy to invert the transform distributions rightly. Distribution with parameter 1. This is not very important. What's important is this is the cdf. This CDF can be inverted using elementary special functions. The elementary special functions are the exponential, the log, the trigonometric ones. Those usually have been optimized by, I guess, the computing community for like decades now and very often have hardware support to. To some extent. So they're like, they're essentially implemented at the CPU level. And this makes it much, much faster. Okay, so we're now going to use basic trigonometry to see how the boxmuller transform works. So we have seen that by the. This is an expression of the Pythagoras theorem, but Pythagoras theorem essentially says that the square of one side of a right angle triangle, one side plus the square of the other, is the square of the hypothenuse. And that's essentially the transformation we were doing here. So we knew that if these two sides were distributed normally, this radius, the circle here would be distributed as the Rayleigh distribution. So now we can go the other way. We can use this to compute first random variables from the Rayleigh distribution, which we, which we found were cheap to compute. And then use trigonometry to project onto the sides of this triangle using the cosine and the sine. And we will get using an auxiliary random variable for the angle which is uniform. If the angle is uniform, the sides of the triangle will be normally distributed. And this can be seen by performing transform distribution on a Rayleigh distribution times the sine or the cosine of an angle that's distributed uniformly. And we do indeed get the the PDF. This is the PDF of the. I'm going a bit slow. The PDF of a normal distribution. Okay, so another reason I wanted to show the box model transform is actually this code is trivial to do annec and run very, very fast at top level as well. So this is the procedure we talked earlier. This is we generate pairs of random variables, uniform random variables. One of them we choose to compute the radius, and then the other one we use to compute the angle. And we project onto X and Y. And this gives us pairs of randomly of normally distributed random variables, as we can see before. And this is computing 2 times 10 to the 7 pairs, as opposed to 10 to the 6 for the inversion formula. So this is over 650 times faster. Okay, so now let's look at a. Let's look at another very general method, which is called acceptance rejection. And in, in its simplest form, it uses two uniform random numbers. For example, suppose the box Muller was using two and giving you two. This is using two and giving you one. However, if the most of the cost of generating these random numbers is in this transformation, so as long as these transformations are cheap enough, this is fine. And these methods can be very fast and efficient. So we use one of the pair of uniform random numbers to propose a possible sample from a hat distribution. And the other pair we used to compare to the ratio of the PDFs. Between a target distribution and this hat distribution, this proposal distribution. And we accept the proposal if it's appropriate. So let's sample from a trapezoidal distribution, for instance, with a uniform Hat and here I'm introducing the trapezes defined here as a piecewise item. Introducing probability distribution which is a container for general purpose user defined probability distributions. And because I haven't bothered with the to normalize this I this F doesn't integrate to one. If I specify method normalize the system will take care of the rest and compute the appropriate coefficient. Okay, so here is what the PDF of the trapezi trapezoidal distribution looks like. And this is the hat that I've chosen for it. And it's important to note that we in real and here this is a toy problem, but in reality we really want the hat to be as tight as possible to the. To the all of this blank space here is not very good. We also want to balance with the fact that we want sampling from the hat has to be easy to do here. It's very easy because it's a uniform. But what you're balancing here is having a hat that is tight as tight as possible is as simple to sample from as possible. And let's write a very simple algorithm that I can show here we generate again 10 to the 6 pairs this we're going to use this dynamic array just to store the accepted proposals. And what we do in the algorithm is we just specify a proposal which normally would be a transformed wouldn't be a uniform but would be a sample from this hat. Then we use this P to compare to this will be a uniform usually in the algorithm. And we compare the PDF of our target distribution the proposal with divided by the PDF of the hat at the proposal which is in this case it's a constant because it's uniform and if it's smaller than if it's bigger than p we append. We do this in a loop and we get the samples that we can see here. The the histogram overlaid with the PDFs and we see that it follows quite nicely. But we also see that we've essentially thrown away 30% of the proposals and this is indeed the ratio of the volume of this PDF here the hat compared sorry the target distribution which here is one because it's normalized to the volume of the hat here. So that's why we need one usually one this hat to be as tight as possible because we don't want to be throwing away 30% of our proposals. Okay, So the probability distribution we saw earlier is a general container for use of functions and it has several sampling methods available, including acceptance rejection and this is true for other built in distributions as well. But the built in distributions has optimized custom codes that is optimized for both speed and correctness. So we wouldn't want to use these, these algorithms, general purpose algorithms for the built in cases. But for if you're defining your own probability distribution, these can. And you for instance, when I use acceptance rejection and you define a suitable hat, this can be a very fast and efficient sampler. So let's look at an example where we have like the product of two sinusoidals divided by 1 over x squared. Essentially the sinc function is sine of x over 1 over x. And this, this complicated distribution that I plot here. The, the orange one has a very natural hat which is just the envelope of the sinusoidals which is the one over x squared. And we can just truncate it here when it reaches close to zero and have a uniform distribution here. So I am gonna, I'm gonna. And here's the piecewise it's constant close to the origin and it's 1 over x further out. So from this hat function we can plug in this distribution to random variants. Select the method acceptance rejection, the sub option hat function. You could also provide a sampler for the hat function which makes it even faster, not needed in this case. Or use other methods like mcmc, which I'll talk about in a different talk. But this is already pretty good. I'm showing the log PDF to show more detail, otherwise it all gets squeezed near the origin a lot. And so this works really nicely for this complicated. It's actually pretty good. Okay, so let's. I'm almost out of time, but I'm almost done. So let's talk a little bit about discrete distributions now and let's. I'm going to show that like a lot of the same algorithms I've been talking about applying, although sometimes they're a little bit more complicated because we're here, we're trying to sample from integers. So we try to transform some uniform random numbers in the real line and truncate them to integers. And that has some subtleties. Like for instance in this algorithm by Newth, which is not one that is used by anyone today, but it, it's. I'm choosing it here because it shows some of the things I want to talk about. So the PDF of the Poisson distribution looks like this. It has a constant factor, there's an exponential and then it has an exponential with base lambda and divided by a factorial. And we've seen the PDF of PMF of the exponential before at the Beginning of the talk. And here we can see that the algorithm by Nouth is very, very simple. We define some constants and then we iterated a loop multiplying in each iteration we specify a uniform random number. We multiply it and when this accumulated product reaches a threshold, we just return the number of the iteration, how many iterations this loop did. And we can see that, we can see that this indeed gives the same as before with the built in sampler, this gives the right result. Actually, there are subtleties here that I'm glossing over. Like for small enough land for specific values of lambda, this algorithm will have numeric problems. But this is good enough for now. It's toy problem and I just wanted to. It's a toy algorithm that I want to be able to show the code. However, if we want to have a look at the real thing, which is what we're using in our system, it's a much more complicated algorithm. It switches between inversion and acceptance rejection and for different ranges of the parameter lambda. The acceptance rejection algorithm is a very clever algorithm by Aren't S and Dieter, which uses a discrete normal hat function. And this hat function is above the Poisson density on one known side of the range. Like on one side of the range here, up to here, it's below, and up to here is above. So you can split the range in two, compute a certain constant to majorize these numbers and then we can use that. And these hat function will actually be very, very close to the actual PDF. So this algorithm will be very efficient and has very few rejections. So I showed a total implementation of the algorithm, but different facets of acceptance rejection. I had a sketch now line of of a more complicated use of acceptance rejection in the Irons Dieter implementation for the Poisson. But acceptance rejection is used in most sampling algorithms and they can get very clever and be very efficient. I showed that it becomes expensive when one of the ways it becomes expensive is when the hat gets very far away from the density when it's not tight. And this unfortunately happens inevitably at higher dimensions due to something called the curse of dimensionality, which is the fact that a lot of them in a high dimensions, a lot of them measure in a shape, although the volume in a shape is near the boundaries. And that's why in high dimensions, very quickly one has to use different algorithms, mostly mcmc, Markov Chain, Monte Carlo, which we'll discuss in another talks, because we mostly use it in other parts of the system, not so much in univariate samplers. Okay, so closing Remarks here. I haven't been looking at the questions because I'm a terrible person. Yes, this is a huge subject which is scratching the surface. I just wanted to like whet the audience's appetite here. I think it's very, very fun and that might be just me. And I deliberately avoided talking about MCMC because this is better addressed in different and subsequent talks that will be given by myself and Gorsia. Okay, thanks so much for listening. Yeah, thanks so much. This is the. I'm done with the talk and now I can hopefully maybe answer some questions. Okay. Okay, so it's Rachel, Andrew, you. Okay, so let me go in in order. So the link between number theory and probability theory. Yeah there's actually quite a lot and again I don't know enough about number theory to really answer this authoritatively but I've read about it. It's just like also like it's a very big topic and I really, I'm not sure I could answer this well enough in the time that we have here. Okay. Andrew Yu. Okay, this is topic for probability. Any plans for this theme of incorporating penalized regression loas we are regression is in the pipeline but it was definitely talked about it. I don't know when it's going to be done though. But we definitely want to revamp progression. There's something that we talked about. Generalized additive models, not so much. When are we going to talk about mcmc? I don't know. I'm guessing there's other people that want to present. I imagine not immediately but I don't know. We've prepared some of these talks already. That will be a question for Charles I guess. Okay, so that seems to be it I think question wise. Okay. Missing something. Okay. Okay, thanks so much. I hope this was made some sense to you guys and also was interesting because that was the point. So yeah, goodbye and thank you so much for listening.