|
1130 | 1130 |
|
1131 | 1131 | $$ r = \sqrt{4y} $$ |
1132 | 1132 |
|
1133 | | -Which means the inverse of our CDF is defined as |
| 1133 | +Which means the inverse of our CDF (which we'll call $ICD(x)$) is defined as |
1134 | 1134 |
|
1135 | | - $$ P^{-1}(r) = \sqrt{4y} $$ |
| 1135 | + $$ P^{-1}(r) = \operatorname{ICD}(r) = \sqrt{4y} $$ |
1136 | 1136 |
|
1137 | 1137 | Thus our random number generator with density $p(r)$ can be created with: |
1138 | 1138 |
|
1139 | | - $$ f(d) = \sqrt{4 \cdot \operatorname{random\_double}()} $$ |
| 1139 | + $$ \operatorname{ICD}(d) = \sqrt{4 \cdot \operatorname{random\_double}()} $$ |
1140 | 1140 |
|
1141 | 1141 | Note that this ranges from 0 to 2 as we hoped, and if we check our work, we replace |
1142 | 1142 | `random_double()` with $1/4$ to get 1, and also replace with $1/2$ to get $\sqrt{2}$, just as |
|
1155 | 1155 | The last time that we tried to solve for the integral we used a Monte Carlo approach, uniformly |
1156 | 1156 | sampling from the interval $[0, 2]$. We didn't know it at the time, but we were implicitly using a |
1157 | 1157 | uniform PDF between 0 and 2. This means that we're using a PDF = $1/2$ over the range $[0,2]$, which |
1158 | | -means the CDF is $P(x) = x/2$, so $f(d) = 2d$. Knowing this, we can make this uniform PDF explicit: |
| 1158 | +means the CDF is $P(x) = x/2$, so $\operatorname{ICD}(d) = 2d$. Knowing this, we can make this |
| 1159 | +uniform PDF explicit: |
1159 | 1160 |
|
1160 | 1161 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ |
1161 | 1162 | #include "rtweekend.h" |
|
1165 | 1166 |
|
1166 | 1167 |
|
1167 | 1168 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight |
1168 | | - double f(double d) { |
| 1169 | + double icd(double d) { |
1169 | 1170 | return 2.0 * d; |
1170 | 1171 | } |
1171 | 1172 |
|
|
1184 | 1185 |
|
1185 | 1186 | for (int i = 0; i < N; i++) { |
1186 | 1187 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight |
1187 | | - auto x = f(random_double()); |
| 1188 | + auto x = icd(random_double()); |
1188 | 1189 | sum += x*x / pdf(x); |
1189 | 1190 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ |
1190 | 1191 | } |
|
1199 | 1200 | </div> |
1200 | 1201 |
|
1201 | 1202 | There are a couple of important things to emphasize. Every value of $x$ represents one sample of the |
1202 | | -function $x^2$ within the distribution $[0, 2]$. We use a function $f$ to randomly select samples |
1203 | | -from within this distribution. We were previously multiplying the average over the interval |
1204 | | -(`sum / N`) times the length of the interval (`b - a`) to arrive at the final answer. Here, we |
1205 | | -don't need to multiply by the interval length--that is, we no longer need to multiply the average |
| 1203 | +function $x^2$ within the distribution $[0, 2]$. We use a function $\operatorname{ICD}$ to randomly |
| 1204 | +select samples from within this distribution. We were previously multiplying the average over the |
| 1205 | +interval (`sum / N`) times the length of the interval (`b - a`) to arrive at the final answer. Here, |
| 1206 | +we don't need to multiply by the interval length--that is, we no longer need to multiply the average |
1206 | 1207 | by $2$. |
1207 | 1208 |
|
1208 | 1209 | We need to account for the nonuniformity of the PDF of $x$. Failing to account for this |
1209 | 1210 | nonuniformity will introduce bias in our scene. Indeed, this bias is the source of our inaccurately |
1210 | | -bright image--if we account for nonuniformity, we will get accurate results. The PDF will "steer" |
| 1211 | +bright image. Accounting for the nonuniformity will yield accurate results. The PDF will "steer" |
1211 | 1212 | samples toward specific parts of the distribution, which will cause us to converge faster, but at |
1212 | 1213 | the cost of introducing bias. To remove this bias, we need to down-weight where we sample more |
1213 | 1214 | frequently, and to up-weight where we sample less frequently. For our new nonuniform random number |
1214 | | -generator, the PDF defines how much or how little we sample a specific portion. |
1215 | | -So the weighting function should be proportional to $1/\mathit{pdf}$. |
1216 | | -In fact it is _exactly_ $1/\mathit{pdf}$. |
1217 | | -This is why we divide `x*x` by `pdf(x)`. |
| 1215 | +generator, the PDF defines how much or how little we sample a specific portion. So the weighting |
| 1216 | +function should be proportional to $1/\mathit{pdf}$. In fact it is _exactly_ $1/\mathit{pdf}$. This |
| 1217 | +is why we divide `x*x` by `pdf(x)`. |
1218 | 1218 |
|
1219 | | -We can try to solve for the integral using the linear PDF $p(r) = \frac{r}{2}$, for which we were |
1220 | | -able to solve for the CDF and its inverse. To do that, all we need to do is replace the functions |
1221 | | -$f = \sqrt{4d}$ and $pdf = x/2$. |
| 1219 | +We can try to solve for the integral using the linear PDF, $p(r) = \frac{r}{2}$, for which we were |
| 1220 | +able to solve for the CDF and its inverse, ICD. To do that, all we need to do is replace the |
| 1221 | +functions $\operatorname{ICD}(d) = \sqrt{4d}$ and $p(x) = x/2$. |
1222 | 1222 |
|
1223 | 1223 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ |
1224 | | - double f(double d) { |
| 1224 | + double icd(double d) { |
1225 | 1225 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight |
1226 | 1226 | return std::sqrt(4.0 * d); |
1227 | 1227 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ |
|
1243 | 1243 | if (z == 0.0) // Ignore zero to avoid NaNs |
1244 | 1244 | continue; |
1245 | 1245 |
|
1246 | | - auto x = f(z); |
| 1246 | + auto x = icd(z); |
1247 | 1247 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ |
1248 | 1248 | sum += x*x / pdf(x); |
1249 | 1249 | } |
|
1290 | 1290 |
|
1291 | 1291 | and |
1292 | 1292 |
|
1293 | | - $$ P^{-1}(x) = f(d) = 8d^\frac{1}{3} $$ |
| 1293 | + $$ P^{-1}(x) = \operatorname{ICD}(d) = 8d^\frac{1}{3} $$ |
1294 | 1294 |
|
1295 | 1295 | <div class='together'> |
1296 | 1296 | For just one sample we get: |
1297 | 1297 |
|
1298 | 1298 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ |
1299 | | - double f(double d) { |
| 1299 | + double icd(double d) { |
1300 | 1300 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight |
1301 | 1301 | return 8.0 * std::pow(d, 1.0/3.0); |
1302 | 1302 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ |
|
1319 | 1319 | if (z == 0.0) // Ignore zero to avoid NaNs |
1320 | 1320 | continue; |
1321 | 1321 |
|
1322 | | - auto x = f(z); |
| 1322 | + auto x = icd(z); |
1323 | 1323 | sum += x*x / pdf(x); |
1324 | 1324 | } |
1325 | 1325 | std::cout << std::fixed << std::setprecision(12); |
|
1342 | 1342 | nonuniform PDF is usually called _importance sampling_. |
1343 | 1343 |
|
1344 | 1344 | In all of the examples given, we always converged to the correct answer of $8/3$. We got the same |
1345 | | -answer when we used both a uniform PDF and the "correct" PDF ($i.e. f(d)=8d^{\frac{1}{3}}$). While |
1346 | | -they both converged to the same answer, the uniform PDF took much longer. After all, we only needed |
1347 | | -a single sample from the PDF that perfectly matched the integral. This should make sense, as we were |
1348 | | -choosing to sample the important parts of the distribution more often, whereas the uniform PDF just |
1349 | | -sampled the whole distribution equally, without taking importance into account. |
| 1345 | +answer when we used both a uniform PDF and the "correct" PDF (that is, $\operatorname{ICD}(d) = |
| 1346 | +8d^{\frac{1}{3}}$). While they both converged to the same answer, the uniform PDF took much longer. |
| 1347 | +After all, we only needed a single sample from the PDF that perfectly matched the integral. This |
| 1348 | +should make sense, as we were choosing to sample the important parts of the distribution more often, |
| 1349 | +whereas the uniform PDF just sampled the whole distribution equally, without taking importance into |
| 1350 | +account. |
1350 | 1351 |
|
1351 | 1352 | Indeed, this is the case for any PDF that you create--they will all converge eventually. This is |
1352 | 1353 | just another part of the power of the Monte Carlo algorithm. Even the naive PDF where we solved for |
1353 | 1354 | the 50% value and split the distribution into two halves: $[0, \sqrt{2}]$ and $[\sqrt{2}, 2]$. That |
1354 | 1355 | PDF will converge. Hopefully you should have an intuition as to why that PDF will converge faster |
1355 | | -than a pure uniform PDF, but slower than the linear PDF ($i.e. f(d) = \sqrt{4d}$). |
| 1356 | +than a pure uniform PDF, but slower than the linear PDF (that is, $\operatorname{ICD}(d) = |
| 1357 | +\sqrt{4d}$). |
1356 | 1358 |
|
1357 | 1359 | The perfect importance sampling is only possible when we already know the answer (we got $P$ by |
1358 | 1360 | integrating $p$ analytically), but it’s a good exercise to make sure our code works. |
|
0 commit comments