似然函数和最大似然估计


拿到一组样本时，在大致判断了它服从什么样的分布以后，最重要的步骤就是求得分布的最优参数，以下简称参数为$\theta$。这里回顾以下最大似然估计法。

1 最大似然估计

最大似然估计，以我的理解，就是找到一个最优参数$\theta$，使得观测到已知样本的可能性最大，是一种反推方法。
对于随机变量$X$，观测到一组相互独立的n维样本$x_1,x_2...x_n$。若$X$在某一取值$x_i$处的概率密度为$f_X(x_i|\theta)$，那么随机变量$X$取值为这一组样本的概率密度为
$\begin{aligned}
f_X(x|\theta)=f_X(x_1,x_2...x_n|\theta)=\prod_{i=1}^{n}f(xi|\theta)
\end{aligned}$
似然估计函数：$L(\theta|x)=f_X(x|\theta)$。在得到了$L(\theta|x)$后，我们想找出的最优参数(即最大似然参数) $\hat\theta_{MV}$就是那个能让$L(\theta|x)$取到最大值的参数，为什么是最大值我们稍后通过例子可以推断出。
$\begin{aligned}
\hat\theta_{MV}=\mathop{\arg\max}_{\theta}L(\theta|x)
\end{aligned}$
如果一个函数一阶可导，二阶倒数<0，那么这一定是一个凸函数，函数的最大值就是在一阶倒数=0时候的取值。

$\begin{aligned}
\left\{
             \begin{array}{lr}
             \frac{\partial L(\theta|x)}{\partial \theta} =0 &\\
             &\\
              \frac{\partial^2 L(\theta|x)}{\partial \theta^2}<0
             \end{array}
\right.
\end{aligned}$
2 举例
上面都是理论和结论，为了更好地理解最大似然估计，举一个例子。假设有一枚均匀的硬币，扔到正面和反面的概率都为0.5。现在扔了2次硬币，都是正面。从概率学的角度上讲，扔一枚硬币2次都是正面的概率为0.5*0.5=0.25.




第一次
第二次
概率




正
正
0.25


正
反
0.25


反
正
0.25


反
反
0.25


现在假设我们不知道硬币是否均匀，也不知道扔到正反面的概率分别为多少，令随机变量X为扔到正面的次数，假设扔到正面的概率为$\theta$，那么扔到反面的概率为$1-\theta$, 那么似然函数为：

$\begin{aligned}
L(\theta|x)=f_X(x|\theta)=\prod_{i=1}^{n}f(xi|\theta)=C_n^aP(X=正面)^aP(X=负面)^{n-a}=C_n^a\theta^a(1-\theta)^{n-a}
\end{aligned}$
$\theta$的取值在$0,1$之间。当观测到扔2次硬币均为正面时，似然函数的值为$C_2^2\theta^2(1-\theta)^0=\theta^2$。由于$L(\theta|x)=\theta^2$是一个凹函数，使得似然函数最大的$\hat\theta_{MV}$的值是1。也就是说，在观察到连续两次扔硬币都是正面朝上的情况下，我们认为硬币投掷时正面朝上的概率为1是最合理的。
合理并不代表正确，这仅是在当前观测条件下最好的预测。当样本数量较大时，例如扔1000次硬币，观察正反两面朝上的情况，求出的$\hat\theta_{MV}$会越来越接近真实值1/2。

参考：《统计学习方法》——李航；

摘要

介绍了EM算法的推导过程、步骤及在三硬币模型中的应用等内容。

正文

1. 什么地方要用EM算法

EM算法主要用于含有隐含变量的模型参数估计问题。如果概率模型的变量都是观测变量，那么给定数据，可以直接用极大似然估计或极大后验估计。当概率模型含有隐含变量时，就需要用EM算法来估计参数了。EM算法就是含有隐含变量的概率模型参数的极大似然估计或极大后验估计。

2. EM算法的推导

给定观测数据Y，目标是极大化观测数据（不完全数据）Y关于参数的对数似然函数，即



 表示在模型参数为时，观测数据Y的概率分布。

 

EM算法通过逐步迭代来逐步近似极大化。假设第i次迭代后的估计值为。下一轮的估计值要使。故



利用Jensen不等式得到下界：



令

是

则

是的一个下界。任何可使增大的，都可以使增加。选择能使当前极大的作为新的值。



所以EM算法就是通过迭代不断求Q函数，并将之极大化，直至收敛。下图为EM算法的直观解释，是的一个下界。



图中，在处两个函数的函数值相等，当前迭代中，要找到使极大化。EM算法在点处重新计算Q函数，进行下次迭代。笔者认为这时的图像发生变化，不再是上图。

从图中可以看出，EM算法不能保证全局最优。

3. EM算法的步骤

分为E步和M步，即求Q函数，极大化Q函数。

输入：观测数据Y，隐变量数据Z，联合分布，条件分布；

输出：模型参数。

第一步：选取参数的初始值；（注：模型为）

第二步（E步）：假设第i迭代后，模型代参数为，在第i+1次迭代的E步，计算



是给定观测数据Y和当前的参数估计下隐变量数据Z的条件概率分布；

第三步（M步）：极大化Q函数，求出；

第四步：重复E步，M步直到收敛。

停止迭代的条件： 或。

4. EM算法的收敛性

EM算法的收敛性是指EM算法得到的估计序列是否收敛，主要由两个定理来证明 。

定理一：为观测数据的似然函数，是EM算法得到的参数估计序列，是对应的似然函数序列，则是单调递增的。

定理二：是观测数据的对数似然函数，是EM算法得到的参数估计序列，是对应的对数似然函数序列。如果有上界，则收敛到某一值；在函数与满足一定条件下，EM算法得到的参数估计序列的收敛值是的稳定点。

这里不再展开讨论。

5. 三硬币模型的EM算法

问题：有三枚硬币（ABC）正面向上的概率分别为，p和q。进行如下试验——先掷A，如果A正面向上则掷B，如果A反面向上则掷C。如此独立地重复做n次试验；记录B和C的结果，正面向上记为1，反面向上记为0。观测结果为：1,1,0,0,1,0，...。

若只能观测到结果，不能观测到掷硬币过程，即每一次的观测结果（1或0）由B或C中的哪枚硬币掷出的是未知的。如此情况下估计三枚硬币正面向上的概率，p和q。观测数据表示为，未观测数据表示为。

一次实验得到的概率为



y是观测变量，结果为1或0；

所以三硬币模型可以写作：



故





EM算法先对参数选初始值，再迭代直至收敛；设第i次迭代的估计值为

在第i+1次，E步的计算为，计算数据来自硬币B的概率



M步计算模型参数的新估计值







注1：本博客只作笔者笔记用途，仅供参考用，若有侵权行为请留言联系，笔者会在第一时间回复。

注2：码字不易，若转载请注明出处。

 

 

文章目录
一 实验背景
1.1 针孔相机模型
1.2 两类参数（外参和内参）
1.3 畸变
1.4 相机标定
1.5 张正友法相机标定
二 实验内容
2.1标定步骤
2.2具体实验过程
三 实验总结

一 实验背景
1.1 针孔相机模型



1.2 两类参数（外参和内参）

1.3 畸变
拍照时通常要在相机的镜头前添加透镜，在相机成像的过程中，透镜会对光线的传播产生影响，从而影响相机的成像效果，产生了畸变。

畸变主要分为两种：径向畸变、切向畸变（畸变是相机本身的特性，与内参无关，标定一次后即可）

切向畸变：由于透镜制造上的缺陷使得透镜本身与图像平面不平行而产生的。

径向畸变：透镜质量原因，光线在远离透镜中心的地方比靠近中心的地方更加弯曲。

切向畸变影响较小，通常只考虑径向畸变。

①桶形畸变



②枕形畸变



对畸变进行校正：


1.4 相机标定
同步标定内部参数和外部参数，一般包括两种策略s:

①光学标定: 利用已知的几何信息(如定长棋盘格)实现参数求解。

②自标定: 在静态场景中利用 structure from motion估算参数。
通过世界坐标集(Xi,Yi,Zi)，以 及它们在图像平面上的投影 坐标集(ui,vi)，计算相机投影 矩阵M中的 11个未知参数



线性方法：



非线性方法：



标定工具需要有极高的精度，包括不同平面的角度、特征点的物理距离等。因此制作标定工具十分困难。
1.5 张正友法相机标定
1.5.1 空间平面的三维点与相机平面二维点的映射

设二维相机平面的像素点为：



三维相机平面的像素点为：



齐次坐标：



描述空间坐标到图像坐标的映射 ：



其中：

s: 世界坐标系到图像坐标系的尺度因子

K: 相机内参矩阵

(u0,v0): 像主点坐标

α, β: 焦距与像素横纵比的融合

γ: 径向畸变参数

由于张正友法标定只需要一张棋盘格，所以可设Z=0，

定义旋转矩阵R的第i列为 ri, 则有：



于是空间到图像的映射可改为：



其中H是描述Homographic矩阵,可通过最小二乘，从角点世界坐标到图像坐标的关系求解。（Homgraphic可描述平面到平面之间的映射关系）

令 H 为 H = [h1 h2 h3]

[h1 h2 h3]=λK[r1 r2 t]

通过上述等式的矩阵运算，根据r1和r2正交，以及归一化的约束可以得到如下等式:





1.5.2 内参求解

B是对称阵，其未知量可表示为6D向量b，由于B矩阵是个对称矩阵，所以可以写成一个6维向量形式：



设H中的第i列为 hi



根据b的定义，可以推导出如下公式：





可以推导出



如果有n组观察图像，则V 是 2n x 6 的矩阵



根据最小二乘定义，Vb = 0 的解是 VTV 最小特征值对应的特征向量。

因此, 可以直接估算出b，后续可以通过b求解内参。

因为B中的未知量为6个，所以当观测平面 n ≥ 3 时，可以得到b的唯一解，当 n = 2时, 一般可令畸变参数γ = 0，当 n = 1时, 仅能估算出α 与 β, 此时一般可假定像 主点坐标 u0 与 v0 为0。



B是通过b构造的对称矩阵

内部参数可通过如下公式计算(cholesky分解)：



1.5.3 最大似然估计

上述的推导结果是理想情况下的解，但在实际操作过程中，由于可能存在噪声等一系列的影响，为了增加结果的可靠性，可以使用最大似然估计(Maximum likelihood estimation)来对上述的结果进行优化。

给定 n 张棋盘格图像，每张图像有 m 个角点

最小化下述公式等同于极大似然估计：



上述非线性优化问题可以利用Levenberg- Marquardt算法求解。
二 实验内容
2.1标定步骤
①打印一张棋盘格A4纸张（黑白间距已知），并贴在一个平板上。

②针对棋盘格拍摄若干张图片。

③在图片中检测特征点（Harris角点）。

④根据角点位置信息及图像中的坐标，求解Homographic矩阵。

⑤利用解析解估计方法计算出5个内部参数，以及6个外部参数。

⑥根据最大似然估计策略，设计优化目标并实验参数的refinement。
2.2具体实验过程
（1）打印一张棋盘格，拍摄多张图片（我拍了10张，拍摄工具为iphone7plus）



（2）操作步骤

matlab实现：

打开matlab，在APP中找到Camera Calibrator



进入后，点击Add images将准备好的图片导入



在弹出来的框内输入每个棋盘格的大小



提取出角点



最后点击Calibration，得到标定结果。
python实现：
import _locale
_locale._getdefaultlocale = (lambda *args: ['zh_CN', 'utf8'])
# -*- coding: cp936 -*-
import cv2
import numpy as np
import glob

criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)

objp = np.zeros((6*9,3), np.float32)
objp[:,:2] = np.mgrid[0:9,0:6].T.reshape(-1,2)  # 将世界坐标系建在标定板上，所有点的Z坐标全部为0，所以只需要赋值x和y

obj_points = []  # 存储3D点
img_points = []  # 存储2D点

images = glob.glob(r'D:/qpgs/2.jpg')
i=0;
for fname in images:
    img = cv2.imread(fname)
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    size = gray.shape[::-1]
    ret, corners = cv2.findChessboardCorners(gray,(9,6), None)
    #print(corners)

if ret:
        obj_points.append(objp)
        corners2 = cv2.cornerSubPix(gray, corners, (5, 5), (-1, -1), criteria)  # 在原角点的基础上寻找亚像素角点
        #print(corners2)
        if [corners2]:
            img_points.append(corners2)
        else:
            img_points.append(corners)
        cv2.drawChessboardCorners(img, (9,6), corners, ret)  # 记住，OpenCV的绘制函数一般无返回值
        i+=1;
        cv2.imwrite('conimg'+str(i)+'.jpg', img)
        cv2.waitKey(4000)

print(len(img_points))
cv2.destroyAllWindows()

# 标定
ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(obj_points, img_points, size, None, None)

print("ret:", ret)
print("mtx:\n", mtx) # 内参数矩阵
print("dist:\n", dist)  # 畸变系数   distortion cofficients = (k_1,k_2,p_1,p_2,k_3)
print("rvecs:\n", rvecs)  # 旋转向量  # 外参数
print("tvecs:\n", tvecs ) # 平移向量  # 外参数

img = cv2.imread(images[2])
h, w = img.shape[:2]
newcameramtx, roi = cv2.getOptimalNewCameraMatrix(mtx,dist,(w,h),1,(w,h))#显示更大范围的图片（正常重映射之后会删掉一部分图像）
print (newcameramtx)

dst = cv2.undistort(img,mtx,dist,None,newcameramtx)
x,y,w,h = roi
dst1 = dst[y:y+h,x:x+w]
cv2.imwrite('calibresult3.png', dst1)
print ("dst的大小为:", dst1.shape)

（3）实验结果

①内参矩阵：



②外参





可视化结果如下

以相机为中心：



以棋盘格为中心：


三 实验总结
1.相机标定的目的是获取摄像机的内参和外参矩阵（同时也会得到每一幅标定图像的选择和平移矩阵），内参和外参系数可以对之后相机拍摄的图像就进行矫正，得到畸变相对很小的图像。

2.通过张正友法标定，我们不能得到焦距(f)和像素的物理尺寸(sx,sy)两个参数。因为我们在求解内参阵A时，求解出的是α和β。

α=f/ sx,β= f/ sy 。分别代表焦距长度上，x轴和y轴像素的个数。

虽然没有求得焦距，但这并不影响在三维坐标恢复时，进行三角运算。因为彼时，我们的计算中用到的也是α和β。

3.张正友法标定只需要一个标板，通过在不同方向多次完整拍照，不需要知道定标板的运动方式，就能直接获得相机的内参和畸变系数。该标定方法精度高于自定标法，且不需要高精度的定位仪器。