<a href="https://colab.research.google.com/github/hhk54250/20MA573-HHK/blob/master/hw5/omc_integral%20and%20mcpi.ipynb" target="_parent"></a>


```
import numpy.matlib
import numpy as np

```


```
import random
A=-1
B= 1
C= 17
i= 0
n= 0
for i in range(10000):
  a= random.uniform(A,B)
  b= random.uniform(A,B)
  if a**2+b**2<1:
    n=n+1
print(str(4*n/10000))
```

    3.14


# Example 3
Given i.i.d $\{\alpha_i : i \in 1, 2, ..., N\}$, we use
$$\bar \alpha_N = \frac{1}{N} \sum_{i=1}^N \alpha_i$$
as its estimator of the mean $\mathbb E[\alpha_1]$ and use
$$\beta_N = \frac{1}{N} \sum_{i=1}^N (\alpha_i - \bar \alpha_N)^2$$
as the estimator of $Var(\alpha_1)$. Suppose $\alpha_1 \in L^4$, then 
- Prove $\beta_N$ is biased. 
- Prove that $\beta_N$ is consistent in $L^2$. 
- Can you propose an unbiased estimator? 

__Proof:__
because $\{\alpha_i\}$ is  iid
\begin{equation}
\begin{aligned}
Bias（\beta_N）= &= \mathbb E[\beta_N] - Var(\alpha_1) \\
\end{aligned}
\end{equation}
so 
\begin{equation}
\begin{aligned}
Bias(\beta_N) &= \mathbb E[\beta_N] - Var(\alpha_1) \\
&= \mathbb E[\frac{1}{N} \sum_{i=1}^N (\alpha_i - \bar \alpha_N)^2] - Var(\alpha_1) \\
&= \frac{1}{N} \mathbb E[\sum_{i=1}^N (\alpha_i^2- 2 \alpha_i \bar \alpha_N + \bar \alpha_N^2)] - Var(\alpha_1) \\
&= \frac{1}{N} \mathbb E[\sum_{i=1}^N \alpha_i^2 - 2N \cdot \bar \alpha_N^2 + N \cdot \bar \alpha_N^2] - Var(\alpha_1) \\
&= \frac{1}{N} \sum_{i=1}^N \mathbb E[\alpha_i^2] - \frac{1}{N}\mathbb E[N \cdot \bar \alpha_N^2] - Var(\alpha_1) \\
&= \mathbb E[\alpha_1^2] - \mathbb E[\bar \alpha_N^2]- Var(\alpha_1) \\
&= \mathbb E[\alpha_1]^2 - \mathbb E[\bar \alpha_N^2] .......................(1)\\
\end{aligned}
\end{equation}
because
\begin{equation}
\begin{aligned}
\mathbb E[\bar \alpha_N] &= \mathbb E[\frac{1}{N} \sum_{i=1}^N \alpha_i] \\
&= \frac{1}{N} \sum_{i=1}^N \mathbb E[\alpha_i] \\
&= \frac{1}{N} \cdot N \cdot \mathbb E[\alpha_1] \\
&= \mathbb E[\alpha_1], \\
\end{aligned}
\end{equation}
so
\begin{equation}
\begin{aligned}
(1)=E[\bar \alpha_N]^2-E[\bar \alpha_N^2]=-Var(\bar \alpha_N)=-Var(\frac{1}{N} \sum_{i=1}^N \alpha_i)= -\frac{1}{N^2} \sum_{i=1}^N Var(\alpha_i)=- \frac{1}{N^2} \cdot N \cdot Var(\alpha_1) = -\frac{1}{N} Var(\alpha_1)\\
\end{aligned}
\end{equation}
βN  is biased


 - To show βN is consistent in L^2, we need to show that 
$$\mathbb E[(\beta_N - Var(\alpha_1))^2] \to 0.$$


- find the unbiased estimator:
we already know if $$\beta_N = \frac{1}{N} \sum_{i=1}^N (\alpha_i - \bar \alpha_N)^2$$
\begin{equation}
\begin{aligned}
Bias（\beta_N）= \mathbb E[\beta_N] - Var(\alpha_1) =-\frac{1}{N} Var(\alpha_1)\\
\end{aligned}
\end{equation}
so
\begin{equation}
\begin{aligned}
\mathbb E[\beta_N] = \frac{N-1}{N} Var(\alpha_1)\\
\mathbb E[\beta_2]= \frac{N}{N-1}E[\beta_N] \\
\beta_2= \frac{1}{N-1} \sum_{i=1}^N (\alpha_i - \bar \alpha_N)^2
\end{aligned}
\end{equation}


```
D=0
m=0
pi_list = []
for i in range(100):
  A=-1
  B= 1
  C= 17
  i= 0
  n= 0
  for z in range(10000):
      a= random.uniform(A,B)
      b= random.uniform(A,B)
      if a**2+b**2<1:
       n=n+1
  pi_list.append(4*n/10000)
  D=4*n/10000+D 
d= D/100
for i in range(100):
  m = (pi_list[i] - d)**2+m
print(str(m/100))
```

    0.00029377369600000047



```
import matplotlib.pyplot as plt
y = []
for t in range (6):
  D=0
  m=0
  pi_list = []
  for i in range(100):
    A=-1
    B= 1
    C= 17
    i= 0
    n= 0
    for z in range((2**(5+t))):
        a= random.uniform(A,B)
        b= random.uniform(A,B)
        if a**2+b**2<1:
         n=n+1
    pi_list.append(4*n/(2**(5+t)))
    D=4*n/(2**(5+t))+D 
  d= D/100
  for i in range(100):
    m = (pi_list[i] - d)**2+m
  y.append(m/100)  
x = [n for n in range(5,11)]
print(y)
plt.loglog(x,y)
```


```
import random
A= 0
B= 1
C= 17
i= 0
n= 0
for i in range(100000000):
  a= random.uniform(A,B)
  if a<0.01:
    n=n+100
  if a>0.01:
    n=n+1  
print(str(n/100000000))
```

    1.98990892



```
import matplotlib.pyplot as plt
y = []
for t in range (6):
  D=0
  m=0
  i_list = []
  for i in range(100):
    A= 0
    B= 1
    C= 17
    i= 0
    n= 0
    for z in range((2**(5+t))):
        a= random.uniform(A,B)
        if a<0.01:
          n=n+100
        if a>0.01:
          n=n+1  
    i_list.append(n/(2**(5+t)))
    D=n/(2**(5+t))+D 
  d= D/100
  for i in range(100):
    m = (i_list[i] - d)**2+m
  y.append(m/100)  
x = [n for n in range(5,11)]
print(y)
plt.loglog(x,y)
```

ipynb
02f7c259f68849a06a95ed7db807b3f7d74d53bd
Jupyter Notebook
2f60bda24e8ee77899f3b1a98739f2cfae0009ea
MIT
hhk54250/20MA573-HHK
hw5/omc_integral and mcpi.ipynb
2f60bda24e8ee77899f3b1a98739f2cfae0009ea
MIT
hhk54250/20MA573-HHK
hw5/omc_integral and mcpi.ipynb
2f60bda24e8ee77899f3b1a98739f2cfae0009ea
MIT
hhk54250/20MA573-HHK
hw5/omc_integral and mcpi.ipynb