diff --git a/.nojekyll b/.nojekyll index ed82e10..c021562 100644 --- a/.nojekyll +++ b/.nojekyll @@ -1 +1 @@ -8679afee \ No newline at end of file +acd55f95 \ No newline at end of file diff --git a/labs/lab0-setup.pdf b/labs/lab0-setup.pdf index b3a8278..158d461 100644 Binary files a/labs/lab0-setup.pdf and b/labs/lab0-setup.pdf differ diff --git a/labs/lab1-submission.pdf b/labs/lab1-submission.pdf index 2499ae8..4f4fa78 100644 Binary files a/labs/lab1-submission.pdf and b/labs/lab1-submission.pdf differ diff --git a/labs/lab2-testing.pdf b/labs/lab2-testing.pdf index 1035973..4babcbb 100644 Binary files a/labs/lab2-testing.pdf and b/labs/lab2-testing.pdf differ diff --git a/labs/lab3-debugging.pdf b/labs/lab3-debugging.pdf index b43e696..b6c1e8d 100644 Binary files a/labs/lab3-debugging.pdf and b/labs/lab3-debugging.pdf differ diff --git a/labs/lab4-reproducibility.pdf b/labs/lab4-reproducibility.pdf index b61f32d..f0639ab 100644 Binary files a/labs/lab4-reproducibility.pdf and b/labs/lab4-reproducibility.pdf differ diff --git a/labs/lab5-codereview.pdf b/labs/lab5-codereview.pdf index fc31e92..56d21bd 100644 Binary files a/labs/lab5-codereview.pdf and b/labs/lab5-codereview.pdf differ diff --git a/labs/lab6-scfparallel.pdf b/labs/lab6-scfparallel.pdf index d11bef8..12071b0 100644 Binary files a/labs/lab6-scfparallel.pdf and b/labs/lab6-scfparallel.pdf differ diff --git a/labs/lab7-pythonvsr.pdf b/labs/lab7-pythonvsr.pdf index d9fda18..3956813 100644 Binary files a/labs/lab7-pythonvsr.pdf and b/labs/lab7-pythonvsr.pdf differ diff --git a/labs/lab8-simulation.pdf b/labs/lab8-simulation.pdf index 4bae4c1..3bae04e 100644 Binary files a/labs/lab8-simulation.pdf and b/labs/lab8-simulation.pdf differ diff --git a/ps/ps1.pdf b/ps/ps1.pdf index d642629..6ec99e8 100644 Binary files a/ps/ps1.pdf and b/ps/ps1.pdf differ diff --git a/ps/ps2.pdf b/ps/ps2.pdf index db327ea..b021dc6 100644 Binary files a/ps/ps2.pdf and b/ps/ps2.pdf differ diff --git a/ps/ps3.pdf b/ps/ps3.pdf index a9c7c49..177d8c8 100644 Binary files a/ps/ps3.pdf and b/ps/ps3.pdf differ diff --git a/ps/ps4.pdf b/ps/ps4.pdf index ba281b7..337b8f2 100644 Binary files a/ps/ps4.pdf and b/ps/ps4.pdf differ diff --git a/ps/ps5.pdf b/ps/ps5.pdf index 7aea9e7..a5f2675 100644 Binary files a/ps/ps5.pdf and b/ps/ps5.pdf differ diff --git a/ps/ps6.pdf b/ps/ps6.pdf index a8f0996..27f56ae 100644 Binary files a/ps/ps6.pdf and b/ps/ps6.pdf differ diff --git a/search.json b/search.json index ce0bc4f..11f7cba 100644 --- a/search.json +++ b/search.json @@ -2013,7 +2013,7 @@ "href": "units/unit9-sim.html#generating-random-uniforms-on-a-computer", "title": "Simulation", "section": "Generating random uniforms on a computer", - "text": "Generating random uniforms on a computer\nGenerating a sequence of random standard uniforms is the basis for all generation of random variables, since random uniforms (either a single one or more than one) can be used to generate values from other distributions. Most random numbers on a computer are pseudo-random. The numbers are chosen from a deterministic stream of numbers that behave like random numbers but are actually a finite sequence (recall that both integers and real numbers on a computer are actually discrete and there are finitely many distinct values), so it’s actually possible to get repeats. The seed of a RNG is the place within that sequence where you start to use the pseudo-random numbers.\n\nSequential congruential generators\nMany RNG methods are sequential congruential methods. The basic idea is that the next value is \\[u_{k}=f(u_{k-1},\\ldots,u_{k-j})\\mbox{mod}\\,m\\] for some function, \\(f\\), and some positive integer \\(m\\) . Often \\(j=1\\). mod just means to take the remainder after dividing by \\(m\\). One then generates the random standard uniform value as \\(u_{k}/m\\), which by construction is in \\([0,1]\\). For our discussion below, it is important to distinguish the state (\\(u\\)) from the output of the RNG.\nGiven the construction, such sequences are periodic if the subsequence ever reappears, which is of course guaranteed because there is a finite number of possible subsequence values given that all the \\(u_{k}\\) values are remainders of divisions by a fixed number . One key to a good random number generator (RNG) is to have a very long period.\nAn example of a sequential congruential method is a basic linear congruential generator: \\[u_{k}=(au_{k-1}+c)\\mbox{mod}\\,m\\] with integer \\(a\\), \\(m\\), \\(c\\), and \\(u_{k}\\) values. (Note that in some cases \\(c=0\\), in which case the periodicity can’t exceed \\(m-1\\) as the method is then set up so that we never get \\(u_{k}=0\\) as this causes the algorithm to break.) The seed is the initial state, \\(u_{0}\\) - i.e., the point in the sequence at which we start. By setting the seed you guarantee reproducibility since given a starting value, the sequence is deterministic. In general \\(a\\), \\(c\\) and \\(m\\) are chosen to be ‘large’. The standard values of \\(m\\) are Mersenne primes, which have the form \\(2^{p}-1\\) (but these are not prime for all \\(p\\)). Here’s an example of a linear congruential sampler (with \\(c=0\\)):\n\nn = 100\na = 171\nm = 30269\n\nu = np.empty(n)\nu[0] = 7306\n\nfor i in range(1, n):\n u[i] = (a * u[i-1]) % m\n\nu = u / m\nuFromNP = np.random.uniform(size = n)\n\nplt.figure(figsize=(10, 8))\n\nplt.subplot(2, 2, 1)\nplt.plot(range(1, n+1), u)\nplt.title(\"manual\")\nplt.xlabel(\"Index\"); plt.ylabel(\"Value\")\n\nplt.subplot(2, 2, 2)\nplt.plot(range(1, n+1), uFromNP)\nplt.title(\"numpy\")\nplt.xlabel(\"Index\"); plt.ylabel(\"Value\")\n\nplt.subplot(2, 2, 3)\nplt.hist(u, bins=25)\nplt.xlabel(\"Value\"); plt.ylabel(\"Frequency\")\n\nplt.subplot(2, 2, 4)\nplt.hist(uFromNP, bins=25)\nplt.xlabel(\"Value\"); plt.ylabel(\"Frequency\")\n\nplt.tight_layout()\nplt.show()\n\n\n\n\n\n\n\n\nA wide variety of different RNG have been proposed. Many have turned out to have substantial defects based on tests designed to assess if the behavior of the RNG mimics true randomness. Some of the behavior we want to ensure is uniformity of each individual random deviate, independence of sequences of deviates, and multivariate uniformity of subsequences. One test of a RNG that many RNGs don’t perform well on is to assess the properties of \\(k\\)-tuples - subsequences of length \\(k\\), which should be independently distributed in the \\(k\\)-dimensional unit hypercube. Unfortunately, linear congruential methods produce values that lie on a simple lattice in \\(k\\)-space, i.e., the points are not selected from \\(q^{k}\\) uniformly spaced points, where \\(q\\) is the the number of unique values. Instead, points often lie on parallel lines in the hypercube.\nCombining generators can yield better generators. The Wichmann-Hill is an option in R and is a combination of three linear congruential generators with \\(a=\\{171,172,170\\}\\), \\(m=\\{30269,30307,30323\\}\\), and \\(u_{i}=(x_{i}/30269+y_{i}/30307+z_{i}/30323)\\mbox{mod}\\,1\\) where \\(x\\), \\(y\\), and \\(z\\) are generated from the three individual generators. Let’s mimic the Wichmann-Hill manually:\n\nRNGkind(\"Wichmann-Hill\")\nset.seed(1)\nsaveSeed <- .Random.seed\nuFromR <- runif(10)\na <- c(171, 172, 170)\nm <- c(30269, 30307, 30323)\nxyz <- matrix(NA, nr = 10, nc = 3)\nxyz[1, ] <- (a * saveSeed[2:4]) %% m\nfor( i in 2:10)\n xyz[i, ] <- (a * xyz[i-1, ]) %% m\nfor(i in 1:10)\n print(c(uFromR[i],sum(xyz[i, ]/m)%%1))\n\n[1] 0.1297134 0.1297134\n[1] 0.9822407 0.9822407\n[1] 0.8267184 0.8267184\n[1] 0.242355 0.242355\n[1] 0.8568853 0.8568853\n[1] 0.8408788 0.8408788\n[1] 0.3421633 0.3421633\n[1] 0.7062672 0.7062672\n[1] 0.6212432 0.6212432\n[1] 0.6537663 0.6537663\n\n## we should be able to recover the current value of the seed\nxyz[10, ]\n\n[1] 24279 14851 10966\n\n.Random.seed[2:4]\n\n[1] 24279 14851 10966\n\n\n\n\nPCG generators\nSomewhat recently O’Neal (2014) proposed a new approach to using the linear congruential generator in a way that gives much better performance than the basic versions of such generators described above. This approach is now the default random number generator in numpy (see numpy.random.default_rng()), called the PCG-64 generator. ‘PCG’ stands for permutation congruential generator and encompasses a family of such generators.\nThe idea of the PCG approach goes like this:\n\nLinear congruential generators (LCG) are simple and fast, but for small values of \\(m\\) don’t perform all that well statistically, in particular having values on a lattice as discussed above.\nUsing a large value of \\(m\\) can actually give good statistical performance.\nApplying a technique called permutation functions to the state of the LCG in order to produce the output at each step (the random value returned to the user) can improve the statistical performance even further.\n\nInstead of using relatively small values of \\(m\\) seen above, in the PCG approach one uses \\(m=2^k\\), for ‘large enough’ \\(k\\), usually 64 or 128. It turns out that if \\(m=2^k\\) then the period of the \\(b\\)th bit of the state is \\(2^b\\) where \\(b=1\\) is the right-most bit. Small periods are of course bad for RNG, so the bits with small period cause the LCG to not perform well. Thankfully, one simple fix is simply to discard some number of the right-most bits (this is one form of bit shift). Note that if one does this, the output of the RNG is based on a subset of the bits, which means that the number of unique values that can be generated is smaller than the period. This is not a problem given we start with a state with a large number of bits (64 or 128 as mentioned above).\nO’Neal then goes further; instead of simply discarding bits, she proposes to either shift bits by a random amount or rotate bits by a random amount, where the random amount is determined by a small number of the initial bits. This improves the statistical performance of the generator. The choice of how to do this gives the various members of the PCG family of generators. The details are fairly complicated (the PCG paper is 50-odd pages) and not important for our purposes here.\n\n\nMersenne Twister\nA commonly used generator (including in both R and Python) is the Mersenne Twister. It’s the default in R and the old default in numpy (see next section for what I mean by “old default”).\nThe Mersenne Twister has some theoretical support, has performed reasonably on standard tests of pseudorandom numbers and has been used without evidence of serious failure. (But note that O’Neal criticizes it in her technical report.) Plus it’s fast (because bitwise operations are fast). The particular Mersenne twister used has a periodicity of \\(2^{19937}-1\\approx10^{6000}\\). Practically speaking this means that if we generated one random uniform per nanosecond for 10 billion years, then we would generate \\(10^{25}\\) numbers, well short of the period. So we don’t need to worry about the periodicity! The state (sometimes also called the seed) for the Mersenne twister is a set of 624 32-bit integers plus a position in the set, where the position is .Random.seed[2] in R and (I think) np.random.get_state()[2] in Python.\nThe Mersenne twister is in the class of generalized feedback shift registers (GFSR). The basic idea of a GFSR is to come up with a deterministic generator of bits (i.e., a way to generate sequences of 0s and 1s), \\(B_{i}\\), \\(i=1,2,3,\\ldots\\). The pseudo-random numbers are then determined as sequential subsequences of length \\(L\\) from \\(\\{B_{i}\\}\\), considered as a base-2 number and dividing by \\(2^{L}\\) to get a number in \\((0,1)\\). In general the sequence of bits is generated by taking \\(B_{i}\\) to be the exclusive or [i.e., 0+0 = 0; 0 + 1 = 1; 1 + 0 = 1; 1 + 1 = 0] summation of two previous bits further back in the sequence where the lengths of the lags are carefully chosen.\nnumpy provides access to the Mersenne Twister via the MT19937 generator; more on this below. It looks like PCG-64 only became available as of numpy version 1.17.\n\n\nThe period versus the number of unique values generated\nThe output of the PCG-64 is 64 bits while for the Mersenne Twister the output is 32 bits. The result is that the generators generate fewer unique values than their periods. This means you could get duplicated values in long runs, but this does not violate the comment about the periodicity of PCG-64 and Mersenne-Twister being longer than \\(2^{64}\\) and \\(2^{32}\\). Why not? Because the two values after the two duplicated numbers will not be duplicates of each other – as noted previously, there is a distinction between the output presented to the user and the state of the RNG algorithm.\n\n\nThe seed and the state\nSetting the seed picks a position in the periodic sequence of the RNG, i.e., in the state of the RNG. The state can be a single number or something much more complicated. As mentioned above, the state for the Mersenne Twister is a set of 624 32-bit integers plus a position in the set. For the PCG-64 in numpy, the state is two numbers – the actual state and the increment (c above). This means that when the user passes a single number as the seed, there needs to be a procedure that deterministically sets the state based on that single number seed. The details of this are not usually well-documented or viewable by the user.\nIdeally, nearby seeds generally should not correspond to getting sequences from the RNG stream that are closer to each other than far away seeds. According to Gentle (CS, p. 327) the input to set.seed() in R should be an integer, \\(i\\in\\{0,\\ldots,1023\\}\\) , and each of these 1024 values produces positions in the RNG sequence that are “far away” from each other. I don’t see any mention of this in the R documentation for set.seed() and furthermore, you can pass integers larger than 1023 to set.seed(), so I’m not sure how much to trust Gentle’s claim. More on generating parallel streams of random numbers below.\nWhen one invokes a RNG without a seed, RNG implementations generally have a method for choosing a seed (often based on the system clock). The numpy documentation says that it “mixes sources of entropy in a reproducible way” to do this.\nGenerators should give you the same sequence of random numbers, starting at a given seed, whether you ask for a bunch of numbers at once, or sequentially ask for individual numbers.\n\nAdditional notes\nThere have been some attempts to generate truly random numbers based on physical randomness. One that is based on quantum physics is http://www.idquantique.com/true-random-number-generator/quantis-usb-pcie-pci.html. Another approach is based on lava lamps!", + "text": "Generating random uniforms on a computer\nGenerating a sequence of random standard uniforms is the basis for all generation of random variables, since random uniforms (either a single one or more than one) can be used to generate values from other distributions. Most random numbers on a computer are pseudo-random. The numbers are chosen from a deterministic stream of numbers that behave like random numbers but are actually a finite sequence (recall that both integers and real numbers on a computer are actually discrete and there are finitely many distinct values), so it’s actually possible to get repeats. The seed of a RNG is the place within that sequence where you start to use the pseudo-random numbers.\n\nSequential congruential generators\nMany RNG methods are sequential congruential methods. The basic idea is that the next value is \\[x_{k}=f(x_{k-1},\\ldots,x_{k-j})\\mbox{mod}\\,m\\] for some function, \\(f\\), and some positive integer \\(m\\) . Often \\(j=1\\). mod just means to take the remainder after dividing by \\(m\\). One then generates the random standard uniform value as \\(u_k = x_{k}/m\\), which by construction is in \\([0,1]\\). For our discussion below, it is important to distinguish the state (\\(u\\)) from the output of the RNG.\nGiven the construction, such sequences are periodic if the subsequence ever reappears, which is of course guaranteed because there is a finite number of possible subsequence values given that all the \\(u_{k}\\) values are remainders of divisions by a fixed number . One key to a good random number generator (RNG) is to have a very long period.\nAn example of a sequential congruential method is a basic linear congruential generator: \\[x_{k}=(ax_{k-1}+c)\\mbox{mod}\\,m\\] with integer \\(a\\), \\(m\\), \\(c\\), and \\(u_{k}\\) values. (Note that in some cases \\(c=0\\), in which case the periodicity can’t exceed \\(m-1\\) as the method is then set up so that we never get \\(x_{k}=0\\) as this causes the algorithm to break.) The seed is the initial state, \\(x_{0}\\) - i.e., the point in the sequence at which we start. By setting the seed you guarantee reproducibility since given a starting value, the sequence is deterministic. In general \\(a\\), \\(c\\) and \\(m\\) are chosen to be ‘large’. The standard values of \\(m\\) are Mersenne primes, which have the form \\(2^{p}-1\\) (but these are not prime for all \\(p\\)). Here’s an example of a linear congruential sampler (with \\(c=0\\)):\n\nn = 100\na = 171\nm = 30269\n\nx = np.empty(n)\nx[0] = 7306\n\nfor i in range(1, n):\n x[i] = (a * x[i-1]) % m\n\nu = x / m\n\nrng = np.random.default_rng(seed=1)\nuFromNP = rng.uniform(size=n)\n\nplt.figure(figsize=(10, 8))\n\nplt.subplot(2, 2, 1)\nplt.plot(range(1, n+1), u)\nplt.title(\"manual\")\nplt.xlabel(\"Index\"); plt.ylabel(\"Value\")\n\nplt.subplot(2, 2, 2)\nplt.plot(range(1, n+1), uFromNP)\nplt.title(\"numpy\")\nplt.xlabel(\"Index\"); plt.ylabel(\"Value\")\n\nplt.subplot(2, 2, 3)\nplt.hist(u, bins=25)\nplt.xlabel(\"Value\"); plt.ylabel(\"Frequency\")\n\nplt.subplot(2, 2, 4)\nplt.hist(uFromNP, bins=25)\nplt.xlabel(\"Value\"); plt.ylabel(\"Frequency\")\n\nplt.tight_layout()\nplt.show()\n\n\n\n\n\n\n\n\nA wide variety of different RNG have been proposed. Many have turned out to have substantial defects based on tests designed to assess if the behavior of the RNG mimics true randomness. Some of the behavior we want to ensure is uniformity of each individual random deviate, independence of sequences of deviates, and multivariate uniformity of subsequences. One test of a RNG that many RNGs don’t perform well on is to assess the properties of \\(k\\)-tuples - subsequences of length \\(k\\), which should be independently distributed in the \\(k\\)-dimensional unit hypercube. Unfortunately, linear congruential methods produce values that lie on a simple lattice in \\(k\\)-space, i.e., the points are not selected from \\(q^{k}\\) uniformly spaced points, where \\(q\\) is the the number of unique values. Instead, points often lie on parallel lines in the hypercube.\nCombining generators can yield better generators. The Wichmann-Hill is an option in R and is a combination of three linear congruential generators with \\(a=\\{171,172,170\\}\\), \\(m=\\{30269,30307,30323\\}\\), and \\(u_{i}=(x_{i}/30269+y_{i}/30307+z_{i}/30323)\\mbox{mod}\\,1\\) where \\(x\\), \\(y\\), and \\(z\\) are generated from the three individual generators. Let’s mimic the Wichmann-Hill manually:\n\nRNGkind(\"Wichmann-Hill\")\nset.seed(1)\nsaveSeed <- .Random.seed\nuFromR <- runif(10)\na <- c(171, 172, 170)\nm <- c(30269, 30307, 30323)\nu <- rep(0, 10)\nxyz <- matrix(NA, nr = 10, nc = 3)\nxyz[1, ] <- (a * saveSeed[2:4]) %% m\nu[1] <- sum(xyz[1, ]/m) %% 1\nfor(i in 2:10) {\n xyz[i, ] <- (a * xyz[i-1, ]) %% m\n u[i] <- sum(xyz[i, ]/m) %% 1\n}\nfor(i in 1:10)\n print(c(uFromR[i], u[i]))\n\n[1] 0.1297134 0.1297134\n[1] 0.9822407 0.9822407\n[1] 0.8267184 0.8267184\n[1] 0.242355 0.242355\n[1] 0.8568853 0.8568853\n[1] 0.8408788 0.8408788\n[1] 0.3421633 0.3421633\n[1] 0.7062672 0.7062672\n[1] 0.6212432 0.6212432\n[1] 0.6537663 0.6537663\n\n## we should be able to recover the current value of the seed\nxyz[10, ]\n\n[1] 24279 14851 10966\n\n.Random.seed[2:4]\n\n[1] 24279 14851 10966\n\n\n\n\nPCG generators\nSomewhat recently O’Neal (2014) proposed a new approach to using the linear congruential generator in a way that gives much better performance than the basic versions of such generators described above. This approach is now the default random number generator in numpy (see numpy.random.default_rng()), called the PCG-64 generator. ‘PCG’ stands for permutation congruential generator and encompasses a family of such generators.\nThe idea of the PCG approach goes like this:\n\nLinear congruential generators (LCG) are simple and fast, but for small values of \\(m\\) don’t perform all that well statistically, in particular having values on a lattice as discussed above.\nUsing a large value of \\(m\\) can actually give good statistical performance.\nApplying a technique called permutation functions to the state of the LCG in order to produce the output at each step (the random value returned to the user) can improve the statistical performance even further.\n\nIn the PCG approach, the state is usually a 64 or 128 bit integer. Instead of using relatively small values of \\(m\\) seen above, in the PCG approach one uses \\(m=2^k\\), where \\(k\\) is either 64 or 128. It turns out that if \\(m=2^k\\) then the period of the \\(b\\)th bit of the state is \\(2^b\\) where \\(b=1\\) is the right-most bit. Small periods are of course bad for RNG, so the bits with small period cause the LCG to not perform well. Thankfully, one simple fix is simply to discard some number of the right-most bits (this is one form of bit shift). Note that if one does this, the output of the RNG is based on a subset of the bits, which means that the number of unique values that can be generated is smaller than the period. This is not a problem given we start with a state with a large number of bits (64 or 128 as mentioned above).\nO’Neal then goes further; instead of simply discarding bits, she proposes to either shift bits by a random amount or rotate bits by a random amount, where the random amount is determined by a small number of the initial bits. This improves the statistical performance of the generator. The choice of how to do this gives the various members of the PCG family of generators. The details are fairly complicated (the PCG paper is 50-odd pages) and not important for our purposes here.\n\n\nMersenne Twister\nA commonly used generator (including in both R and Python) is the Mersenne Twister. It’s the default in R and the old default in numpy (see next section for what I mean by “old default”).\nThe Mersenne Twister has some theoretical support, has performed reasonably on standard tests of pseudorandom numbers and has been used without evidence of serious failure. (But note that O’Neal criticizes it in her technical report.) Plus it’s fast (because bitwise operations are fast). The particular Mersenne twister used has a periodicity of \\(2^{19937}-1\\approx10^{6000}\\). Practically speaking this means that if we generated one random uniform per nanosecond for 10 billion years, then we would generate \\(10^{25}\\) numbers, well short of the period. So we don’t need to worry about the periodicity! The state (sometimes also called the seed) for the Mersenne twister is a set of 624 32-bit integers plus a position in the set, where the position is .Random.seed[2] in R and (I think) np.random.get_state()[2] in Python.\nThe Mersenne twister is in the class of generalized feedback shift registers (GFSR). The basic idea of a GFSR is to come up with a deterministic generator of bits (i.e., a way to generate sequences of 0s and 1s), \\(B_{i}\\), \\(i=1,2,3,\\ldots\\). The pseudo-random numbers are then determined as sequential subsequences of length \\(L\\) from \\(\\{B_{i}\\}\\), considered as a base-2 number and dividing by \\(2^{L}\\) to get a number in \\((0,1)\\). In general the sequence of bits is generated by taking \\(B_{i}\\) to be the exclusive or [i.e., 0+0 = 0; 0 + 1 = 1; 1 + 0 = 1; 1 + 1 = 0] summation of two previous bits further back in the sequence where the lengths of the lags are carefully chosen.\nnumpy provides access to the Mersenne Twister via the MT19937 generator; more on this below. It looks like PCG-64 only became available as of numpy version 1.17.\n\n\nThe period versus the number of unique values generated\nThe output of the PCG-64 is 64 bits while for the Mersenne Twister the output is 32 bits. The result is that the generators generate fewer unique values than their periods. This means you could get duplicated values in long runs, but this does not violate the comment about the periodicity of PCG-64 and Mersenne-Twister being longer than \\(2^{64}\\) and \\(2^{32}\\). Why not? Because the two values after the two duplicated numbers will not be duplicates of each other – as noted previously, there is a distinction between the output presented to the user and the state of the RNG algorithm.\n\n\nThe seed and the state\nSetting the seed picks a position in the periodic sequence of the RNG, i.e., in the state of the RNG. The state can be a single number or something much more complicated. As mentioned above, the state for the Mersenne Twister is a set of 624 32-bit integers plus a position in the set. For the PCG-64 in numpy, the state is two numbers – the actual state and the increment (c above). This means that when the user passes a single number as the seed, there needs to be a procedure that deterministically sets the state based on that single number seed. The details of this are not usually well-documented or viewable by the user.\nIdeally, nearby seeds generally should not correspond to getting sequences from the RNG stream that are closer to each other than far away seeds. According to Gentle (CS, p. 327) the input to set.seed() in R should be an integer, \\(i\\in\\{0,\\ldots,1023\\}\\) , and each of these 1024 values produces positions in the RNG sequence that are “far away” from each other. I don’t see any mention of this in the R documentation for set.seed() and furthermore, you can pass integers larger than 1023 to set.seed(), so I’m not sure how much to trust Gentle’s claim. More on generating parallel streams of random numbers below.\nWhen one invokes a RNG without a seed, RNG implementations generally have a method for choosing a seed (often based on the system clock). The numpy documentation says that it “mixes sources of entropy in a reproducible way” to do this.\nGenerators should give you the same sequence of random numbers, starting at a given seed, whether you ask for a bunch of numbers at once, or sequentially ask for individual numbers.\n\nAdditional notes\nThere have been some attempts to generate truly random numbers based on physical randomness. One that is based on quantum physics is http://www.idquantique.com/true-random-number-generator/quantis-usb-pcie-pci.html. Another approach is based on lava lamps!", "crumbs": [ "Units", "Unit 9 (Simulation)" diff --git a/sitemap.xml b/sitemap.xml index b04588d..16334e1 100644 --- a/sitemap.xml +++ b/sitemap.xml @@ -94,7 +94,7 @@ https://stat243.berkeley.edu/fall-2024/units/unit9-sim.html - 2024-10-27T17:45:05.802Z + 2024-10-30T15:07:16.923Z https://stat243.berkeley.edu/fall-2024/units/regex.html diff --git a/units/regex.pdf b/units/regex.pdf index 4b75827..ffb2bfe 100644 Binary files a/units/regex.pdf and b/units/regex.pdf differ diff --git a/units/unit1-intro.pdf b/units/unit1-intro.pdf index 31889fc..54b26ff 100644 Binary files a/units/unit1-intro.pdf and b/units/unit1-intro.pdf differ diff --git a/units/unit2-dataTech.pdf b/units/unit2-dataTech.pdf index 6315dc7..042a318 100644 Binary files a/units/unit2-dataTech.pdf and b/units/unit2-dataTech.pdf differ diff --git a/units/unit3-bash.pdf b/units/unit3-bash.pdf index c4e0a15..53622dc 100644 Binary files a/units/unit3-bash.pdf and b/units/unit3-bash.pdf differ diff --git a/units/unit4-goodPractices.pdf b/units/unit4-goodPractices.pdf index fe43db5..f561d4f 100644 Binary files a/units/unit4-goodPractices.pdf and b/units/unit4-goodPractices.pdf differ diff --git a/units/unit5-programming.pdf b/units/unit5-programming.pdf index 9512169..ad3c2cb 100644 Binary files a/units/unit5-programming.pdf and b/units/unit5-programming.pdf differ diff --git a/units/unit6-parallel.pdf b/units/unit6-parallel.pdf index b43e73b..f1f9853 100644 Binary files a/units/unit6-parallel.pdf and b/units/unit6-parallel.pdf differ diff --git a/units/unit7-bigData.pdf b/units/unit7-bigData.pdf index a6a6821..87fdff1 100644 Binary files a/units/unit7-bigData.pdf and b/units/unit7-bigData.pdf differ diff --git a/units/unit8-numbers.pdf b/units/unit8-numbers.pdf index 203e731..f1651e4 100644 Binary files a/units/unit8-numbers.pdf and b/units/unit8-numbers.pdf differ diff --git a/units/unit9-sim.html b/units/unit9-sim.html index f46ae25..b335bad 100644 --- a/units/unit9-sim.html +++ b/units/unit9-sim.html @@ -690,45 +690,47 @@

G

Generating a sequence of random standard uniforms is the basis for all generation of random variables, since random uniforms (either a single one or more than one) can be used to generate values from other distributions. Most random numbers on a computer are pseudo-random. The numbers are chosen from a deterministic stream of numbers that behave like random numbers but are actually a finite sequence (recall that both integers and real numbers on a computer are actually discrete and there are finitely many distinct values), so it’s actually possible to get repeats. The seed of a RNG is the place within that sequence where you start to use the pseudo-random numbers.

Sequential congruential generators

-

Many RNG methods are sequential congruential methods. The basic idea is that the next value is \[u_{k}=f(u_{k-1},\ldots,u_{k-j})\mbox{mod}\,m\] for some function, \(f\), and some positive integer \(m\) . Often \(j=1\). mod just means to take the remainder after dividing by \(m\). One then generates the random standard uniform value as \(u_{k}/m\), which by construction is in \([0,1]\). For our discussion below, it is important to distinguish the state (\(u\)) from the output of the RNG.

+

Many RNG methods are sequential congruential methods. The basic idea is that the next value is \[x_{k}=f(x_{k-1},\ldots,x_{k-j})\mbox{mod}\,m\] for some function, \(f\), and some positive integer \(m\) . Often \(j=1\). mod just means to take the remainder after dividing by \(m\). One then generates the random standard uniform value as \(u_k = x_{k}/m\), which by construction is in \([0,1]\). For our discussion below, it is important to distinguish the state (\(u\)) from the output of the RNG.

Given the construction, such sequences are periodic if the subsequence ever reappears, which is of course guaranteed because there is a finite number of possible subsequence values given that all the \(u_{k}\) values are remainders of divisions by a fixed number . One key to a good random number generator (RNG) is to have a very long period.

-

An example of a sequential congruential method is a basic linear congruential generator: \[u_{k}=(au_{k-1}+c)\mbox{mod}\,m\] with integer \(a\), \(m\), \(c\), and \(u_{k}\) values. (Note that in some cases \(c=0\), in which case the periodicity can’t exceed \(m-1\) as the method is then set up so that we never get \(u_{k}=0\) as this causes the algorithm to break.) The seed is the initial state, \(u_{0}\) - i.e., the point in the sequence at which we start. By setting the seed you guarantee reproducibility since given a starting value, the sequence is deterministic. In general \(a\), \(c\) and \(m\) are chosen to be ‘large’. The standard values of \(m\) are Mersenne primes, which have the form \(2^{p}-1\) (but these are not prime for all \(p\)). Here’s an example of a linear congruential sampler (with \(c=0\)):

+

An example of a sequential congruential method is a basic linear congruential generator: \[x_{k}=(ax_{k-1}+c)\mbox{mod}\,m\] with integer \(a\), \(m\), \(c\), and \(u_{k}\) values. (Note that in some cases \(c=0\), in which case the periodicity can’t exceed \(m-1\) as the method is then set up so that we never get \(x_{k}=0\) as this causes the algorithm to break.) The seed is the initial state, \(x_{0}\) - i.e., the point in the sequence at which we start. By setting the seed you guarantee reproducibility since given a starting value, the sequence is deterministic. In general \(a\), \(c\) and \(m\) are chosen to be ‘large’. The standard values of \(m\) are Mersenne primes, which have the form \(2^{p}-1\) (but these are not prime for all \(p\)). Here’s an example of a linear congruential sampler (with \(c=0\)):

n = 100
 a = 171
 m = 30269
 
-u = np.empty(n)
-u[0] = 7306
+x = np.empty(n)
+x[0] = 7306
 
 for i in range(1, n):
-    u[i] = (a * u[i-1]) % m
+    x[i] = (a * x[i-1]) % m
 
-u = u / m
-uFromNP = np.random.uniform(size = n)
-
-plt.figure(figsize=(10, 8))
+u = x / m
+
+rng = np.random.default_rng(seed=1)
+uFromNP = rng.uniform(size=n)
 
-plt.subplot(2, 2, 1)
-plt.plot(range(1, n+1), u)
-plt.title("manual")
-plt.xlabel("Index"); plt.ylabel("Value")
-
-plt.subplot(2, 2, 2)
-plt.plot(range(1, n+1), uFromNP)
-plt.title("numpy")
-plt.xlabel("Index"); plt.ylabel("Value")
-
-plt.subplot(2, 2, 3)
-plt.hist(u, bins=25)
-plt.xlabel("Value"); plt.ylabel("Frequency")
-
-plt.subplot(2, 2, 4)
-plt.hist(uFromNP, bins=25)
-plt.xlabel("Value"); plt.ylabel("Frequency")
-
-plt.tight_layout()
-plt.show()
+plt.figure(figsize=(10, 8)) + +plt.subplot(2, 2, 1) +plt.plot(range(1, n+1), u) +plt.title("manual") +plt.xlabel("Index"); plt.ylabel("Value") + +plt.subplot(2, 2, 2) +plt.plot(range(1, n+1), uFromNP) +plt.title("numpy") +plt.xlabel("Index"); plt.ylabel("Value") + +plt.subplot(2, 2, 3) +plt.hist(u, bins=25) +plt.xlabel("Value"); plt.ylabel("Frequency") + +plt.subplot(2, 2, 4) +plt.hist(uFromNP, bins=25) +plt.xlabel("Value"); plt.ylabel("Frequency") + +plt.tight_layout() +plt.show()
@@ -746,12 +748,16 @@

Sequent uFromR <- runif(10) a <- c(171, 172, 170) m <- c(30269, 30307, 30323) -xyz <- matrix(NA, nr = 10, nc = 3) -xyz[1, ] <- (a * saveSeed[2:4]) %% m -for( i in 2:10) - xyz[i, ] <- (a * xyz[i-1, ]) %% m -for(i in 1:10) - print(c(uFromR[i],sum(xyz[i, ]/m)%%1))

+u <- rep(0, 10) +xyz <- matrix(NA, nr = 10, nc = 3) +xyz[1, ] <- (a * saveSeed[2:4]) %% m +u[1] <- sum(xyz[1, ]/m) %% 1 +for(i in 2:10) { + xyz[i, ] <- (a * xyz[i-1, ]) %% m + u[i] <- sum(xyz[i, ]/m) %% 1 +} +for(i in 1:10) + print(c(uFromR[i], u[i]))
[1] 0.1297134 0.1297134
 [1] 0.9822407 0.9822407
@@ -784,7 +790,7 @@ 
PCG generators

 
  • Using a large value of \(m\) can actually give good statistical performance.
    • 
 
  • Applying a technique called permutation functions to the state of the LCG in order to produce the output at each step (the random value returned to the user) can improve the statistical performance even further.
    • 
 
-
    Instead of using relatively small values of \(m\) seen above, in the PCG approach one uses \(m=2^k\), for ‘large enough’ \(k\), usually 64 or 128. It turns out that if \(m=2^k\) then the period of the \(b\)th bit of the state is \(2^b\) where \(b=1\) is the right-most bit. Small periods are of course bad for RNG, so the bits with small period cause the LCG to not perform well. Thankfully, one simple fix is simply to discard some number of the right-most bits (this is one form of bit shift). Note that if one does this, the output of the RNG is based on a subset of the bits, which means that the number of unique values that can be generated is smaller than the period. This is not a problem given we start with a state with a large number of bits (64 or 128 as mentioned above).
    
+
    In the PCG approach, the state is usually a 64 or 128 bit integer. Instead of using relatively small values of \(m\) seen above, in the PCG approach one uses \(m=2^k\), where \(k\) is either 64 or 128. It turns out that if \(m=2^k\) then the period of the \(b\)th bit of the state is \(2^b\) where \(b=1\) is the right-most bit. Small periods are of course bad for RNG, so the bits with small period cause the LCG to not perform well. Thankfully, one simple fix is simply to discard some number of the right-most bits (this is one form of bit shift). Note that if one does this, the output of the RNG is based on a subset of the bits, which means that the number of unique values that can be generated is smaller than the period. This is not a problem given we start with a state with a large number of bits (64 or 128 as mentioned above).
    
 
    O’Neal then goes further; instead of simply discarding bits, she proposes to either shift bits by a random amount or rotate bits by a random amount, where the random amount is determined by a small number of the initial bits. This improves the statistical performance of the generator. The choice of how to do this gives the various members of the PCG family of generators. The details are fairly complicated (the PCG paper is 50-odd pages) and not important for our purposes here.
    diff --git a/units/unit9-sim.pdf b/units/unit9-sim.pdf index ef16c18..58975e1 100644 Binary files a/units/unit9-sim.pdf and b/units/unit9-sim.pdf differ diff --git a/units/unit9-sim_files/figure-html/Time series and histograms of pseudo-random numbers-3.png b/units/unit9-sim_files/figure-html/Time series and histograms of pseudo-random numbers-3.png index 0b376ae..ce31645 100644 Binary files a/units/unit9-sim_files/figure-html/Time series and histograms of pseudo-random numbers-3.png and b/units/unit9-sim_files/figure-html/Time series and histograms of pseudo-random numbers-3.png differ diff --git a/units/unit9-sim_files/figure-html/unnamed-chunk-2-1.png b/units/unit9-sim_files/figure-html/unnamed-chunk-2-1.png index 45f5c9e..c48cf40 100644 Binary files a/units/unit9-sim_files/figure-html/unnamed-chunk-2-1.png and b/units/unit9-sim_files/figure-html/unnamed-chunk-2-1.png differ