5.2. Applied Stochastic Analysis’s home work 2

5.2.1. No.1

omitting..

5.2.2. No.2

Testify the half order convergence of MC through a numerical example.

Here is the code:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
function ans=hw2_2()
for i = 1:10
    res(i)=estimate();
end
ans = mean(res);

end

function ans=estimate()
syms x;
f = sin(pi*x);
I = int(f,0,1);
N = [4,8,16,32,64,128,256,512]*10;
for i = 1:8
    res(i)=I-mc(N(i),f,x);
end
res=abs(double(res));
ans=polyfit(log(N),log(res),1);
ans = -ans(1);
end

function ans=mc(n,f,x)

ans = mean(subs(f,x,rand(n,1)));
end

5.2.3. No.3

How many ways can you give to construct the uniform distribution on S^2 . Implement them and make a comparison.

Set \theta and \phi be the spherical coordinates. Here I got three mechod:

  1. Let \theta be the uniform distribution of [0,2\pi) and P_{\phi}=\frac{1}{2}sin\phi

x &= \cos \theta \sin \phi\\
y &= \sin \theta \sin \phi\\
z &= \cos \phi

  1. pick u=cos\phi to be uniformly distributed on [-1,1], and \theta the same as above.

x &= \sqrt{1-u^2} \cos\theta\\
y &= \sqrt{1-u^2} \sin\theta\\
z &= u

  1. Marsaglia (1972) derived an elegant method that consists of picking x_1 and x_2 from independent uniform distributions on (-1,1) and rejecting points for which x_1^2+x_2^2\geq 1 . From the remaining points.

x &=2 x_1 \sqrt{1-x_1^2-x_2^2}\\
y &=2 x_2 \sqrt{1-x_1^2-x_2^2}\\
z &=1-2(x_1^2+x_2^2)\\

Here is the code:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
function hw2_3()
    
    num =10000;
    subplot(2,2,1);
    [a,b,c]=sphere();
    mesh(a,b,c);
    hold
    Title('Method a');
    tic
    [x,y,z]=sphere1(num);
    toc
    plot3(x,y,z,'.');


    subplot(2,2,2);
    [a,b,c]=sphere();
    mesh(a,b,c);
    hold
    Title('Method b');
    tic
    [x,y,z]=sphere2(num);
    toc
    plot3(x,y,z,'.');
    
    subplot(2,2,3);
    [a,b,c]=sphere();
    mesh(a,b,c);
    hold
    Title('Method c');
    tic
    [x,y,z]=sphere3(num);
    toc
    plot3(x,y,z,'.');
end

function [x,y,z]=sphere1(n)
    x=zeros(1,n);
    y=zeros(1,n);
    z=zeros(1,n);
    for i=1:n
        theta = 2*pi*rand;
        phi = acos(1-2*rand);
        s = sin(phi);
        x(i) = cos(theta)*s;
        y(i) = sin(theta)*s;
        z(i) = cos(phi);
    end
end

function [x,y,z]=sphere2(n)

    x=zeros(1,n);
    y=zeros(1,n);
    z=zeros(1,n);
    for i=1:n

        theta = 2*pi*rand;
        u = rand*2-1;
        s = sqrt(1-u^2);
        x(i) = s*cos(theta);
        y(i) = s*sin(theta);
        z(i) = u;
    end
    
end

function [x,y,z]=sphere3(n)

    x=zeros(1,n);
    y=zeros(1,n);
    z=zeros(1,n);
    for i=1:n
        
        x1 =2*rand-1;
        x2 =2*rand-1;
        p = x1^2+x2^2;
        while p>=1 
            x1 =2*rand-1;
            x2 =2*rand-1;
            p = x1^2+x2^2;
        end
        s =sqrt(1-p);
        x(i) = 2*x1*s;
        y(i) = 2*x2*s;
        z(i) = 1-2*p;    
    end
end

And here is the result:

alternate text

It’s obvius that the method c is the most effient, but not stable.

a: Elapsed time is 0.026833 seconds.

b: Elapsed time is 0.020886 seconds.

c: Elapsed time is 0.006588 seconds.

PS: I. Another easy way to pick a random point on a sphere is to generate three Gaussian random variables. #. Cook (1957) extended a method of von Neumann (1951) to give a simple method of picking points uniformly distributed on the surface of a unit sphere. This method only need multiply add, sub and divide.