6.2. Applied Stochastic Analysis’s home work 2

6.2.1. No.1

omitting..

6.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

6.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\)
\[\begin{split}x &= \cos \theta \sin \phi\\ y &= \sin \theta \sin \phi\\ z &= \cos \phi\end{split}\]
  1. pick \(u=cos\phi\) to be uniformly distributed on \([-1,1]\), and \(\theta\) the same as above.
\[\begin{split}x &= \sqrt{1-u^2} \cos\theta\\ y &= \sqrt{1-u^2} \sin\theta\\ z &= u\end{split}\]
  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.
\[\begin{split}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)\\\end{split}\]

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.