exa1901 exa1902

Author: SmirkCao

CH19/README.md

@@ -73,9 +73,27 @@ $p(x)=\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x-\mu)^2}{2\sigma^2}}$
 #### 数学期望估计
 
 #### 积分计算
+**例19.1**
+$$
+\int_0^1e^{-{x^2}/{2}}\mathrm{d}x
+$$
+![fig_e1901](assets/fig_e1901.png)
+实际上这个例子，数据点比较少的时候积分效果不是特别好。十个点，均匀分布出来的也不是特别均匀。
+注意，按照**均匀分布**抽取，和**整个区间均分10等分**是不一样的。
+
+**例19.2**
+10个样本
+![fig_e1902_1](assets/fig_e1902_1.png)
+50个样本
+![fig_e1902](assets/fig_e1902.png)
+
+这两个例子，**注意**一下$f(x)$和$p(x)$的定义。
+
 
 ### 马尔可夫链
+
 ### 马尔可夫链蒙特卡罗法
+
 ### Metropolis-Hastings算法
 ### Gibbs Sampling
 

CH19/assets/fig_e1901.png

@@ -1,5 +1,6 @@
 #! /usr/bin/env python
-#!-*- coding=utf-8 -*-
+
+# -*- coding=utf-8 -*-
 # Project:  Lihang
 # Filename: mcmc
 # Date: 6/10/19

CH19/assets/fig_e1902.png

@@ -16,10 +16,66 @@
 
 class TestMCMCMethod(unittest.TestCase):
     def test_exa_1901(self):
-        pass
+        def f(x):
+            return np.exp(-x**2/2)
+
+        def p(x):
+            return 1
+
+        x = np.random.uniform(0, 1, 10)
+        y = f(x)
+        x_raw = np.arange(0, 1, 0.01)
+        y_raw = f(x_raw)
+        print(x.shape, y.shape)
+        data = np.concatenate((x_raw[:, np.newaxis], y_raw[:, np.newaxis]),
+                              axis=1)
+        print(data)
+        # data = data[data[:, 0].argsort()]
+        rst = np.mean(y)
+
+        plt.rc('text', usetex=True)
+        plt.rc('font', family='serif')
+
+        plt.figure(figsize=(5, 5))
+        plt.scatter(x, y, marker='*')
+        plt.plot(data[:, 0], data[:, 1], alpha=0.3)
+        bbox = dict(boxstyle="round", fc="0.8")
+        comment = r"$\int_0^1e^{-{x^2}/{2}}\mathrm{d}x=$"+str(round(rst, 2))
+        plt.annotate(comment, (0.2, 0.8), bbox=bbox)
+        plt.ylabel(r'$\exp\left(-\frac{1}{2}x^2\right)$')
+        plt.ylim(-0.02, 1.02)
+        plt.xlim(-0.02, 1.02)
+        plt.show()
 
     def test_exa_1902(self):
-        pass
+        def p(x):
+            return np.sqrt(2*np.pi)*np.exp(-x**2/2)
+
+        def f(x):
+            return x
+
+        x = np.random.normal(0, 1, 10)
+        y = f(x)
+        x_raw = np.arange(-2, 2, 0.1)
+        y_raw = f(x_raw)
+        print(x.shape, y.shape)
+        data = np.concatenate((x_raw[:, np.newaxis], y_raw[:, np.newaxis]),
+                              axis=1)
+        print(data)
+        # data = data[data[:, 0].argsort()]
+        rst = np.mean(y)
+
+        plt.rc('text', usetex=True)
+        plt.rc('font', family='serif')
+
+        plt.figure(figsize=(5, 5))
+        plt.scatter(x, y, marker='*', alpha=0.3)
+        plt.plot(data[:, 0], data[:, 1], alpha=0.3)
+        bbox = dict(boxstyle="round", fc="0.8")
+        comment = r"$\int_{-\infty}^{\infty}x\frac{1}{\sqrt{2\pi}}e^{-{x^2}/{2}}\mathrm{d}x=$"+str(round(rst, 2))
+        plt.annotate(comment, (-2, 2), bbox=bbox)
+
+        plt.show()
 
     def test_exa_1903(self):
         pass