Graphics Reference
In-Depth Information
Inline Exercise 31.9:
(a) Show that the expected value of the output of this
small program is exactly
A
+
B
.
(b) Modify the program by writing
if (u < 0.3)...
, and adjusting the two
0.5s as needed to get a different estimator for the same value.
(c) Suppose that
A
=
8 and
B
=
12. What fraction
p
would you use to replace
0.5 in line 4 to ensure that the estimates produced by multiple runs of the
program had the smallest possible variance?
(d) What if
A
=
0 and
B
=
12?
Now let's suppose that
B
is in fact a sum of
n
terms:
B
=
B
1
+
...
B
n
. We can
modify the program to handle this:
1
2
3
4
5
6
7
define estimate():
u
= uniform(0, 1)
if (
u
< 0.5):
return
A
/
0.5
else:
i
= randint(1,
n
)
// random integer from
1
to
n
.
return
B
i
/ (0.5
*
(1/
n
))
We're just using a uniform single-sample estimate of
B
now as well as the earlier
random choice to estimate
A
+
B
.
Let's apply these ideas to estimating the solution
x
to
(
I
−
M
)
x
=
b
, which is
x
=
b
+
Mb
+
M
2
b
+
...
(31.90)
=
b
+
M
(
b
+
Mb
+
...
)
.
(31.91)
To keep things simple in what follows, and in analogy with the light-transport
application we'll soon examine, we'll assume that all entries of
M
are non-
negative, the row-sum
r
i
=
j
m
ij
is less than 1 for all
i
, and the eigenvalues
of
M
have magnitude less than 1, so that the series solution converges. We'll also
assume that
M
is an
n
2 matrix.
Let's use the ideas of the short programs above to estimate
x
1
, the first entry
of
x
in Equation 31.90. According to the equation, the value
x
1
is the sum of two
numbers,
b
1
and
(
Mx
)
1
. That makes it easy to find
x
1
using our sum-of-numbers
estimator, but only if we already know
x
!But
Mx
1
is a sum:
m
11
x
1
+
×
n
matrix rather than just a 2
×
+
m
1
n
x
n
.
We can estimate this by selecting among these terms at random, as usual. In each
term, we know the
m
ij
factor, and we need to
estimate
the
x
j
factor. So we can
estimate
x
1
only if we can estimate
x
j
for any
j
. This begs for recursion. And as is
typical in recursion, the recursion gets easier if we broaden the problem. So we'll
write a procedure,
estimate(i)
, that estimates the
i
th entry of
x
, and to find
x
1
we'll invoke it for
i
=
1.
...
1
2
3
4
5
6
7
define estimate(i):
u
= uniform(0, 1)
if (
u
> 0.5):
return
b
i
/ 0.5
else:
k
= randint(1,
n
)
return
m
ik
·
estimate(
k
) / (0.5
*
(
1
/
n
))
The topic's website has actual code that implements this approach, and com-
pares it with the answer obtained by actually solving the system of equations.
