python数值计算 ECE 302 Probabilistic Methods in Electrical and Computer Engineering

ECE 302: Probabilistic Methods in Electrical and Computer Engineering Fall 2018
Instructor: Prof. Stanley H. Chan

Project 1: Numerical Integration by Monte Carlo Sampling

Fall 2018 (Due: Sep 28, 2018)

Project is due at 4:30pm. Please put your project report in the dropbox located at MSEE 330. No late report will be accepted.

1 Introduction

The Monte Carlo method is a very popular simulation method that relies on repeated random sampling to obtain numerical results. Monte Carlo has a wide range of applications in physics, optimization, and probability. You will see Monte Carlo very often if you study further in machine learning, operation research, electro-magnetic simulation, etc. In this project, we will study some basic ideas of Monte Carlo.

To begin with, let us consider a two-dimensional function g(x) defined on some set Ω = {x | h(x) ≥ 0}. If we like to compute the integral

􏰋

Ω

􏰋

I = =

g(x)dx
g(x)dx, (1)

h(x)≥0
we will soon encounter some difficulties because g(x) and h(x) may not allow us to derive a closed form

expression. For example, suppose x = (x1,x2) and g ( x ) = 1 e − x 21 + x 2 2

2

2π
h(x)=x31 +x2 +2x1x2 −sinx1 +cos(2πx1x2)−1, (2)

how would you compute I then? (See Figure 1 for a pictorial illustration.)

Figure 1: Integrating a high-dimensional function g(x) over some nonlinear set Ω = {x | h(x) ≥ 0} can be a challenging problem.

© 2018 Stanley Chan. All Rights Reserved. 1

2 Monte Carlo Method
2.1 Classical Numerical Integration

To understand Monte Carlo, let us first think about what we would do in the classical numerical integration. Suppose we have a one-dimensional integration

􏰋b
I= g(x)dx (3)

a

over some interval [a, b]. To integrate g, we can divide the interval [a, b] into N equally spaced sub-intervals. The width of each sub-interval is b−a . Then the integration is approximated by a summation

N 􏰋 b

g(x)dx ≈ a

􏰊N 􏰃 b − a 􏰄

N g(xn) n=1

| b − a | 􏰊N

g(xn)I(xn ∈[a,b]),
where xn is a point in the nth sub-interval so that g(xn) is the height of that sub-interval. The function

= N I(xn ∈ [a, b]) is the indicator function, defined as

􏰇1, if xn ∈[a,b], I{xn ∈ [a,b]} = 0, if xn ̸∈ [a,b].

Figure 2 illustrates the idea: each rectangle has an area of 􏰁b−a􏰂g(xn). The total area is the sum of all N

these rectangles.

Figure 2: In classical numerical integration, we divide the interval [a, b] into N equally spaced sub-intervals, and calculate the height in each subinterval.

2.2 Monte Carlo Method

The Monte Carlo method is a simple modification of the numerical integration. All it does is to replace the equally spaced intervals by some randomly spaced intervals. That is, if we draw N i.i.d. random variables X(1),…,X(N) uniformly distributed over the sample space S (e.g., S = [0,1]2), then the integration can be approximated by the finite sum

􏰋 | S | 􏰊N

g(x)dx ≈ N Ω

g(X(n))I{X(n) ∈ Ω}, (4) n=1

n=1

© 2018 Stanley Chan. All Rights Reserved.

2

where I{X(n) ∈ Ω} is an indicator function with the property that
􏰇1, if X(n) ∈ Ω,

I{X(n) ∈Ω}= 0, if X(n) ̸∈Ω.

The magic of Monte Carlo is about the random samples X(1), . . . , X(N). To see how these random samples can lead to the desired solution, we can take the expectation on the right hand side and show that

􏰅 | S | 􏰊N

E N
n=1

n=1 Ω The equalities above hold because of the following reasons:

• 􏰎  (a): By linearity of expectation: E[X + Y ] = E[X] + E[Y ]
• 􏰎  (b): X (1) , . . . , X (N ) are i.i.d. random variables. So we only need

􏰆 ( a ) | S | 􏰊N 􏰌 􏰍 g(X(n))I{X(n) ∈ Ω} = N E g(X(n))I{X(n) ∈ Ω}

n=1
( b ) | S | 􏰊N 􏰋 􏰃 1 􏰄

= N n=1 S |S| g(x)I{x ∈ Ω}dx ( c ) | S | 􏰊N 􏰋 􏰃 1 􏰄

=Nn=1Ω |S|g(x)dx

X |S| 􏰉 g(x)f (n) (x)dx for any function g.

( d ) 1 􏰊N 􏰋

( e ) 􏰋
g(x)dx = g(x)dx.

(5)

them, say

= N

Ω

to consider one of

X(n). X(n) is a uniform random variable defined over the sample space S. Therefore, its probability

density function (PDF) is f (n)(x) = 1 . By the definition of expectation, we have E[g(X(n))] =

SX
􏰎 (c): The indicator function changes the integration domain from S to Ω, because for all x ∈ Ω\S, we

have I{X(n) ∈ Ω} = 0.
􏰎 (d): Cancelation between the factor |S| outside the integration and the factor 1 inside the integration.

|S|
􏰎 (e): There are N identical integrals in this summation. So the factor N cancels with the factor 1/N.

Perhaps complicated, the equation in (4) can actually be implemented in a few simple steps: 􏰎 Step 1: Generate random variables X(1),…,X(N);

􏰎 Step 2: Check whether X(n) ∈ Ω.
􏰎 Step 3: Compute the sum 􏰈Nn=1 g(X(n))I{X(n) ∈ Ω}.

You may wonder why do people like about Monte Carlo. It turns out that Monte Carlo can be equipped with some elegant sampling strategies to speed up the numerical integration (e.g., Markov Chain Monte Carlo). These sampling strategies can allow us to use very small N and achieve some results as good as the full numerical integration. We will not touch these sampling strategies here. What we will do is to go through some basic steps of Monte Carlo in this project.

3 Exercise

Your task in this project is to write a Monte Carlo algorithm. Please follow the report instruction posted on the website when preparing your report.

© 2018 Stanley Chan. All Rights Reserved. 3

1. (10 points) First of all, generate N = 1000 random variables X(n) = (X(n),X(n)) representing 1000 12

   random coordinates in the unit square S = [0, 1]2. In your report, state clearly how you did it in Python

2. (10 points) Next, construct the set Ω using the h(x) given in (2). Write a function my_h_function that takes a vector of N coordinates x and outputs a vector output of N values of h(x).

   Similarly, for the function g(x) given in (2), write a function my_g_function that takes a vector of N coordinates x and outputs a vector output of N values of g(x).

   Highlight the region of g(x) where h(x) ≥ 0 by

   Submit your plot.

3. (20 points) Write a program to estimate the integral

   1 􏰊N

   g(X(n))I{X(n) ∈ Ω}. Your function should have the following format

   Submit your code. What is the estimate of 􏰉

4. (20 points) Use your Monte Carlo algorithm to estimate the value of π. (Hint: Set h(x) for a circle;

   Set g(x) = 1.)

x = numpy.random.rand(…) # a vector of N two-dimensional coordinates in the range [0,1]^2.

def my_h_function(x): output = …

return output

def my_g_function(x): output = …

return output

I = N

n=1

Ω

g(x)dx?

© 2018 Stanley Chan. All Rights Reserved. 4

5. (20 points) The experiment from part 1 – part 3 is conducted by drawing N = 1000 random coordinates. Therefore, I itself is a random variable and so it has mean and standard deviation. Repeat the experi- ment for T = 500 times to estimate the mean and standard deviation. Furthermore, repeat the calcula- tion of means and standard deviations for a range of N ’s in numpy.round(numpy.logspace(1,5,100)). Plot the raw estimate I, the mean numpy.mean(I,1), and the standard deviation numpy.std(I,1) as a function of N. Below is a sample code:

Submit your code and the plot. Explain the trend shown in the plot. How large should N be if we want the standard deviation to be within at most 5% of the mean?

6. (20 points) Prove that
􏰅 | S | 􏰊N 􏰆 | S | 􏰋 1 􏰃 􏰋 􏰄 2

g(X(n))I{X(n) ∈ Ω} = N g(x)2dx − N g(x)dx n=1 Ω Ω

Var N
SupposeNislargeenough,e.g.,N=100000.ForS=[0,1]2,estimate􏰉 g(x)2dx.Togetherwiththe

. (6)

result in Part 3, estimate the value
V =|S| g(x)2dx− g(x)dx .

Ω 􏰋 􏰃􏰋 􏰄2

ΩΩ
Plot |S| 􏰉 g(x)2dx− 1 􏰁􏰉 g(x)dx􏰂2 as a function of N. Overlap this curve with the empirical variance

NΩ NΩ
found in Part 5. A sample code is as follows.

Submit your proof (hand write or typed are both acceptable). Submit your code and the plot. Explain the difference between the two curves.

© 2018 Stanley Chan. All Rights Reserved. 5