|
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