多维正态分布中的边际分布、条件分布及Bayes公式

这个是概率论的一个期末小论文,大概删掉了开头结尾一些废话放在这里。

这篇文章主要是想要证明几个关于多维正态分布的很有用的定理,一是已知一个多维正态分布的联合分布,求边际分布以及条件分布的公式。二是已知一个条件分布和先验分布求后验分布的Bayes公式。这几个公式在进行统计的时候会比较有用。

预警:下面会有大量的公式,如果有公式恐惧症的话看每一节开头的结论就好了。

(more…)

Read More

三维空间中的旋转:旋转矩阵、欧拉角

考虑这样一个问题:如何计算三维空间中一个点绕着某一条向量旋转一个特定角度之后的坐标?旋转矩阵、欧拉角和四元数都是用来解决这个问题的方法。

接下来我们来讨论一下旋转矩阵和欧拉角这两个方法,并且我们选取右手坐标系作为我们的坐标系。

旋转矩阵

首先,对于一个三维空间的点 P(x, y, z),要将其绕 z 轴旋转 \theta 角度是可以很简单地用旋转矩阵来表示的

 \begin{equation}
R_z(\theta) =
\begin{bmatrix}
\cos \theta & -\sin \theta & 0 \\
\sin \theta & \cos \theta & 0 \\
0 & 0 & 1 \\
\end{bmatrix}
\end{equation}

类似地,绕另外两个坐标轴旋转的矩阵可以表示如下(注意 R_y 有些不同)

 \begin{equation}
R_x(\theta) =
\begin{bmatrix}
1 & 0 & 0 \\
0 & \cos \theta & -\sin \theta \\
0 & \sin \theta & \cos \theta \\
\end{bmatrix}
\end{equation}

\begin{equation}
R_y(\theta) =
\begin{bmatrix}
\cos \theta & 0 & \sin \theta \\
0 & 1 & 0 \\
-\sin \theta & 0 & \cos \theta \\
\end{bmatrix}
\end{equation}

对于这三个特殊的旋转轴我们已经有了解决方案了,那么对于任意的轴 p 要怎么办呢?我们可以将这个旋转分解:

  • 将整个坐标轴旋转,使得旋转轴 pz 轴重合
  • 再将点 Pz 轴旋转 \theta
  • 再将整个坐标轴旋转回原位

3d-angle-rotation-matrix

如图,我们可以用两个角 \phi\psi 表示一个旋转轴的位置(这里认为旋转轴是个单位向量

 \begin{equation}
\left \{ \begin{aligned}
x &= \sin \phi \cos \psi \\
y &= \sin \phi \sin \psi \\
z &= \cos \phi
\end{aligned} \right.
\end{equation}

这样整个旋转就可以表示如下(最先进行的旋转它对应的旋转矩阵在最右侧)

R_z(\psi)R_y(\phi)R_z(\theta)R_y(-\phi)R_z(-\psi)

或者利用 R(-\alpha) = R^{-1}(\alpha) = R^T(\alpha),可以改写为

R_z(\psi)R_y(\phi)R_z(\theta)R_y^T(\phi)R_z^T(\psi)

经过化简,就可以得到最终的旋转矩阵

\begin{equation}
\begin{bmatrix}
\cos \theta + x^2(1 - \cos \theta) & xy(1 - \cos \theta) - z\sin \theta & xz(1 - \cos \theta) + y\sin \theta \\
yx(1 - \cos \theta) + z\sin \theta & \cos \theta + y^2(1 - \cos \theta) & yz(1 - \cos \theta) - x\sin \theta \\
zx(1 - \cos \theta) - y\sin \theta & zy(1 - \cos \theta) + x\sin \theta & \cos \theta + z^2(1 - \cos \theta)
\end{bmatrix}
\end{equation}

我们将旋转矩阵左乘需要旋转的向量就可以得到旋转后的结果了!

旋转矩阵一个很方便的地方是它可以沿着任意轴任意角度的旋转,但是,旋转矩阵缺点是它需要有 9 个元素来表示一个旋转,而且矩阵的乘法也比较慢。

欧拉角

euler-angle

欧拉角是用三个旋转角度 \alpha, \beta, \gamma 来标示旋转的。如图,图中蓝色坐标系是起始的坐标系,红色的坐标系是最后旋转完成的坐标系。整个旋转分为三个步骤:

  • 将坐标系绕 z 轴旋转 \alpha
  • 将旋转后坐标系绕自己本身的 x 轴(也就是图中的 N 轴)旋转 \beta
  • 将旋转后坐标系绕自己本身的 z 轴旋转 \gamma

由于绕不同的轴旋转所得到的欧拉角是不同的,所以欧拉角在使用的时候必须要先指明旋转的顺序,这里使用的是“zxz”的顺序。

欧拉角表示的旋转转换成旋转矩阵就是

 \begin{equation} R_z(\alpha)R_x(\beta)R_z(\gamma) \end{equation}

需要注意,这里的后面两次旋转并不是在原本固定的坐标系下的旋转。在旋转矩阵中,每一次旋转的叠加都是在左边乘上对应的旋转矩阵,然而,在这里乘法的顺序是相反的。

这可以这样来理解,假设最原始的固定的坐标系是 C_0

  • 假设有一个和 C_0 重合坐标系 C_3,先将 C_3C_0z 轴旋转 \gamma
  • 假设有一个和 C_0 重合坐标系 C_2,将 C_2 和前一步旋转后的 C_3 一起绕 C_0x 轴旋转 \beta
  • 假设有一个和 C_0 重合坐标系 C_1,将 C_1 和前一步旋转后的 C_2, C_3 一起绕 C_0z 轴旋转 \alpha

然后我们仅看这三个坐标系的关系:C_0 绕自己的 z 轴旋转 \alpha 角就可以和 C_1 重合;C_1 绕自己的 x 轴旋转 \beta 角可以和 C_2 重合,这是因为 C_1C_2 在最后一步是一起旋转的,它们的相对位置不会改变;同样可以知道 C_2 绕自己的 z 轴旋转 \gamma 角就可以和 C_3 重合。

在这里 C_1, C_2 相当于是前面欧拉角旋转的前两步得到的坐标系。

这样的话从 C_0C_3 的变换就相当于之前欧拉角的旋转变换,因此按照这个过程,旋转矩阵就是按照上面的顺序相乘了。

Read More

反演魔术:反演原理及二项式反演

首先我来说说什么是反演(inversion),对于一个数列\{f_n\}来说,如果我们知道另外一个数列\{g_n\},满足如下条件

 g_n = \sum_{i=0}^n a_{ni}f_i

反演过程就是利用 g_0, g_1, \cdots, g_n 来表示出 f_n

 f_n = \sum_{i = 0}^n b_{ni}g_i

本质上来说,反演其实是一个解线性方程组的过程,但是你观察后发现,这个矩阵实际上是一个三角矩阵,那必然存在着快捷的方法。对于使得这两个反演公式成立的这些系数应该满足什么条件呢?我们现在就来探讨一番!

我们首先讨论一下在什么情况下能够比较容易建立反演公式,之后介绍二项式反演以及它的两个应用,其中一个是错位排列问题

(more…)

Read More

特殊多项式在整点上的线性插值方法

首先我来解释一下题目是什么东西,假如给你一个 k 次多项式 P(x),并且已知了 P(0), P(1), \cdots, P(k),要求计算出 P(n), n \in \mathbb N 的值,由于已经知道了 k + 1 个点值,那么要插值出系数是可以在 \mathcal O(k^2) 完成,如果使用 FFT 的话是可以在 \mathcal O(k \log^3 k) 的时间内完成的[1],但是对于这个特别的问题,我们有一个 \mathcal O(k) 的算法!这好像是杜教先提出来的,然后还有这篇文章

(more…)

Read More

Catalan 数列及其应用

卡特兰数(Catalan number)是组合数学中一个重要的计数数列,在很多看似毫不相关地方都能见到它的身影

这篇文章介绍卡特兰数的几个应用以及它的推导过程,从组合推理和生成函数两个方面来推导出 Catalan 数的公式

带限制条件的路径总数

首先我们来看一个问题:

在一个平面直角坐标系中,只能往右或往上走一个单位长度,问有多少种不同的路径可以从左下角 (1, 1) 走到右上角 (n, n),并且要求路径不能经过直线 y = x 上方的点,下图中的路径都是合法的(图片来源 Wikipedia)

450px-Catalan_number_4x4_grid_example.svg

如果没有限制条件,那么从左下角走到右上角一共有 2n 步,有 n 步是往右,另外 n 步是往上,那么路径方案数就是 2n 步中选择 n 步往右,一共有 {2n} \choose {n} (即 C_{2n}^n)种方案

那么我们考虑一下这里面有多少种方案是不合法的

首先对于每一种不合法的方案,它的路径一定与 y = x + 1 有交。我们找到它与 y = x + 1 的第一个交点,然后将这个点后面部分的路径关于 y = x + 1 做一个对称。由于原来路径到达 (n, n),新的对称之后的路径就会到达 (n - 1, n + 1)。这样我们把每一种不合法方案都对应到了一条从 (1, 1)(n - 1, n + 1) 的路径,现在再来看是否每一条这样的路径都能对应到一种不合法方案,如果是,那么这就建立了一个一一映射的关系,也就是它们的方案总数相同。这是肯定的,因为每一条这样的路径必定与 y = x + 1 有交,那么对称回去,就得到一条不合法方案

由于从 (1, 1)(n - 1, n + 1) 的路径有 {n - 1 + n + 1} \choose {n - 1} 条,那么合法的方案就是

 C_n = {{2n} \choose {n}} - {{2n} \choose {n - 1}} = \frac{1}{n + 1} {{2n} \choose {n}}

我们把这个方案数记为 C_n,这就是著名的 Catalan 数

(more…)

Read More

Codeforces 上的一些组合计数问题

Codeforces 15E. Triangles

Problem:给出一个有规律的金字塔(如图),黑色的边是可以走路径,要求求出包含 H 点的简单回路(不经过一个顶点两次)的条数,并且这个回路围住的部分不可以包含灰色的三角形,答案对 10^9+9 取模,其中 n \leq 10^6,并且一定是一个偶数

5945210be972aa5fe947f3b8a3a0378a4cade844

Example:当 n = 2 时,一共有 10 种方案,其中 5 种是这样,另外的是顺时针顺序地走

5E-triangles

Solution:首先第 1 层是最特殊的,因为这里没有灰色的三角形,然后注意到灰色的部分把整个金字塔分成了左部和右部

5E-triangles-level1

路径现在可以分成三个部分

  • 在第一层移动,这一共 10 种(见样例)
  • 进入左部或右部后直接回到 H
  • 进入左部(右部)后经过 E 点进入右部(左部)

对于第二种情况,在第一层的部分走法可以是(只有逆时针部分,实际上有 8 种)

  • A \rightarrow B \rightarrow D \rightarrow \text{LeftPart} \rightarrow E \rightarrow C \rightarrow A
  • A \rightarrow B \rightarrow D \rightarrow \text{LeftPart} \rightarrow E \rightarrow F \rightarrow C \rightarrow A
  • A \rightarrow B \rightarrow E \rightarrow \text{RightPart} \rightarrow F \rightarrow C \rightarrow A
  • A \rightarrow B \rightarrow D \rightarrow E \rightarrow \text{RightPart} \rightarrow F \rightarrow C \rightarrow A

对于第三种情况,在第一层(原来图中的两层现在算一层,也就是下面的 n 是题目中 n\frac{1}{2})的部分走法有 2 种,左部到右部和右部到左部

现在要统计的就是在左部和右部内有多少种方案,因为左部右部对称,所以答案是一样的,我们来考虑左部

S_n 表示在有 n 层的金字塔中,走左部的方案数,因为走进左部后每一层都会遇到凹进去的部分,这部分只要枚举一下走了多远可以统计出第 n 层一共有 2\sum_{i=1}^{n-2} 2^i + 1 = 2^n - 3 这么多的情况

然后要走到第 n 层后就是 \prod_{i=2}^{n} (2^i - 3) 种情况,然后走回去的路径有 4 种选择,因此

 S_n = \sum_{i=2}^{n} 4\prod_{j=2}^i (2^j - 3)

答案是 2S_n^2 + 8S_n + 10 (more…)

Read More

牛顿迭代法在多项式运算的应用

总算是快把 FFT 和生成函数的各种东西补了好多,膜拜策爷的论文和 Picks 的博客 QAQ

这篇文章大概就是说如何用牛顿迭代法来解满足 G(F(z)) \equiv 0 \pmod {z^n}F(z)

然后这东西可以比较方便地计算 \sqrt{A(z)}e^{A(z)},也就是多项式开根、求指数对数之类鬼畜的东西,在生成函数计数中十分有用

顺便一提,这里说的“多项式”实际上你可以直接理解为生成函数或者形式幂级数

(more…)

Read More

多项式的多点求值与快速插值

多项式的多点求值(multi-point evaluation)是给出一个多项式 A(x),和 n 个点 x_0, x_1, \cdots, x_{n-1},要求求出 A(x_0), A(x_1), \cdots, A(x_{n-1})

相反的,多项式的插值(interpolation)是给出 n+1 个点 (x_0, y_0), (x_1, y_1), \cdots, (x_n, y_n),求出一个 n 次多项式,使得这些点都在这个多项式上

这两个问题实际上是在多项式的点值表示(point-value representation)和系数表示(coefficient representation)之间转换的方法,在快速傅里叶变换中由于带入值的特殊性质,可以在 \mathcal O(n\log n) 的之间内将两种东西互相转换,但是,如果是任意给定点要求求值或者插值,就没有比较好的性质可以利用,但是仍然有比较快速的方法来计算它们

(more…)

Read More

多项式除法及求模

问题是这样的:给定一个 n 次多项式 A(x) 和一个 m(m \leq n) 次多项式 B(x),要求求出两个多项式 D(x), R(x),满足

 \begin{equation} \label{div0} A(x) = D(x)B(x) + R(x) \end{equation}

并且 deg D \leq degA - degB = n - mdegR < m

在解决这个问题之前你需要掌握多项式求逆,利用快速傅里叶变换可以在 \mathcal O(n\log n) 的复杂度内求出这个问题的解

多项式除法在很多问题中都有应用,例如多项式的扩展欧几里得算法、线性递推的优化等等

(more…)

Read More