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

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

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

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

(more…)

Read More

Vim里一些好用的插件

我主要是想要把一些用过的好用的插件记下来,万一哪天换了系统还可以在这里找回来。

另外这篇文章会不断更新,如果我发现新的有意思的插件也会加进来。

YouCompleteMe

绝对好用的自动补全插件,支持C/C++等各类语言的自动补全

cpp-demo-of-youcompleteme

具体请戳:http://blog.miskcoo.com/2015/12/vim-plugin-youcompleteme

Jedi-Vim

一个很好用的Python的自动补全插件

具体请戳:https://github.com/davidhalter/jedi-vim

Vim-Latex

这是一款帮助编辑Latex的插件,里面集成了很多环境的快捷输入方式,例如输入EFI就会自动变成如下

另外还能够编译并且方便地查看错误信息等。

具体请戳:http://vim-latex.sourceforge.net

Markdown-Preview

一个国人写的插件,支持Markdown的同步预览,感觉挺好用的。

具体请戳:https://github.com/iamcco/markdown-preview.vim

Read More

较短步数复原魔方的算法

这其实是我们上学期徐 BOSS 程设课的一个大作业,我本来很早就想把它记录下来了,可是实在太懒了加上拖延症严重于是就拖到了这个学期。本来作业是要求复原魔方就好了,但是后来发现没有图形界面的话似乎调试很困难,就顺便写了一个 UI,刚好还学了一下 OpenGL。最开始我是写了能够找到最优解的 Krof 算法,但是估计是我写得比较丑再加上状态空间实在太大了(大约在搜索到第 15 步还是 16 步的时候状态数就到百亿级别了,而我 8 个线程一秒只能搜几千万个状态好像),根本没有办法在可以接受的时间内得到解。于是退而求其次,改为去搜索一个比较优的解,最后基本上一秒以内就可以得到解了。

整个工程的代码我放在了我的 Github 上:https://github.com/miskcoo/rubik-cube

魔方的表示

如何表示一个魔方呢?一个很直接的方法就是用六个面一共 54 个色块来表示它,我在最开始写的时候就是用这种方法来表示一个魔方的,但是这样存在着很多的冗余。事实上,我们可以按照魔方的各个棱块和角块的位置和方向来存储。对于魔方的一个角块来说,我们只需要存储这个块在 8 个角的哪一个地方以及它的方向如何(关于方向接下来会进行规定)。对于棱块来说也是相同的,只需要位置和方向就可以决定了。

角块方向的表示

首先对于一个角块来说它可能会有 3 个不同的方向(如下图,暂不考虑这个魔方是否合法,仅关注角块)

corner-0corner-2

想象一下,从左到右相当于单独把面对我们的这个角块(按照从魔方中心指向该角块顶点的向量)旋转了两次。我们可以按照这样来规定角块的方向:首先将对应角块像上图中一样面对我们,并且以绿色或蓝色作为顶面。现在以绿色或蓝色作为参照色(因为一个角块有且仅有这两种颜色之一)。现在来规定角块的方向:

  • 如果参照色位置在如同第一张图一样的位置,就规定方向为 0
  • 如果参照色位置在如同第二张图一样的位置,就规定方向为 1
  • 否则就规定方向为 2

棱块方向的表示

对于棱块方向的规定就稍微有些复杂了,因为没有办法像角块这样找到对立的两个面使得每个棱都包含这两种颜色之一。那么对于棱块的参照色我们规定:如果一个棱存在蓝色或绿色,则该颜色为参照色,否则白色或者黄色为参照色。规定了参照色之后就可以规定棱块的方向了:

  • 如果棱块位于顶层或者底层(也就是绿色或者蓝色面所在的层),那么当参照色在这两个面中的一个时方向为 0,否则为 1
  • 如果棱块位于中间层,那么当参照色在白色或者黄色面时方向为 0,否则为 1

edge-0edge-1

在第一张图中,两个棱块的方向都是 0。这是因为对于红白色棱块,白色为参照色并且其位于顶层绿色面;对于绿白色棱块,绿色为参照色并且位于中层白色面。

在第二张图中,两个棱块的方向都是 1。原因与之前的刚好相反。

Kociemba 算法

首先我们规定一种魔方的中间状态:仅仅通过 U、D、L2、R2、F2、B2 这几种旋转就可以将魔方复原的状态。也就是不允许 90^\circ 地旋转魔方的前后左右这四个面。这种状态有什么特点呢?

首先很明显可以观察到按照我们刚刚规定的方向,旋转顶面和底面都是不会改变棱块和角块的方向的,并且仅按照 180^\circ 来旋转周围四个面也是不会改变棱块和角块的方向。那也就是说,在这种中间状态下,魔方的棱块和角块的方向应该都是 0,这是因为目标状态的魔方如此。除此之外,魔方中间层的四个棱块必然也在中间层。

可以简单地计算出这种中间状态一共只有 8! \times 8! \times 4! = 3.9 \times 10^{10} 种可能,事实证明,在不超过 18 步以内就可以从这种状态复原魔方。我们用 IDA* 来进行搜索,当我们设计启发函数的时候分别针对角块(一共 8! 种可能)、顶层和底层的棱块(一共 8! 种可能)、中间层的棱块(一共 4! 种可能)计算复原它们所需要的最少步数(这可以直接进行搜索),之后的启发函数就是这三者的最大值。

剩下的问题就是如何将魔方旋转到中间状态,我们同样采用 IDA* 来进行搜索。对于启发函数,分别针对角块的方向(一共 3^8 种可能)、棱块的方向(一共 2^8 种可能)、中间层棱块的位置及方向(一共 A^4_{12}\times 2^4 种可能)计算复原它们所需要的最少步数。经过计算,不超过 12 步就可以将一个魔方旋转到中间状态。这样总共在不超过 30 步的情况下就可以复原一个魔方了!

刚刚的状态表示方式在这个算法里显示出了很大的优势,在进行 Hash 的时候可以直接使用这个状态表示,而不必从其它表示方法(例如存储 54 个面的方法)计算出棱和角的方向。事实证明,在通常情况下,这个算法很快就能够复原一个魔方。

Krof 算法

这个算法就是我前面说的搜索最优解的算法,思路和 Kociemba 算法几乎一样,不过 Krof 算法没有划分阶段,启发函数是分别对于角块的方向和位置(一共 8!\times3^7 种可能)、两组每组 6 个梭块的位置和方向(每组都是 A^6_{12}\times2^6 种可能)计算复原它们所需要的最少步数。由于这个启发函数的计算也是十分耗时的,所以可以先计算一次之后记录下来之后直接读取就可以了,而前面 Kociemba 算法的启发函数可以现场计算。

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

上帝的指纹——Mandelbrot集合的绘制

322px-mandel_zoom_00_mandelbrot_set

曼德博集合(Mandelbrot Set)是一个在复平面上的点集。有人认为 Mandelbrot 集合是“人类有史以来做出的最奇异、最瑰丽的几何图形”,曾被称为“上帝的指纹”。

Mandelbrot 集合是一个分形(fractal),将它无限放大都能够有精妙的细节在内,而这瑰丽的图案仅仅由一个简单的公式生成。

Mandelbrot 集合的定义是由法国数学家 Adrien Douady 做出的,而它的命名则是为了纪念被称为“分形学之父”的 Benoit Mandelbrot。

在计算机出现后,对于分形的绘制成为了可能。在这篇文章中,我会介绍如何用 C++ 以及高精度库 MPFR 绘制 Mandelbrot 集合,以及对图像色彩渲染的一些方法。

具体代码可以参考 https://github.com/miskcoo/mandelbrot-render

(more…)

Read More

列表、字符串及字典——Python基本数据结构

Python 的强大之处就在于它有丰富的库以及各种方便的特性。这篇文章将会介绍 Python 中的几个内置类型——列表、元组、字符串及字典。

我是按照 Python3 来介绍,如果和 Python2 有不同的地方会特别指出的。并且,文中的代码如果有输出的部分,都是按照 Python3 语法进行。另外,如果代码中是以 >>> 开头的,就说明这段代码是我直接在交互模式下运行的,并且包含了输出。

List——列表类型

首先,我们来介绍列表类型。什么是列表呢?你可以认为列表是多个变量的集合,并且这个集合是有序的。你在以前遇到的类型可能都只是涉及单个元素,比如说 intfloat 之类的,这都只涉及了一个数字。但是列表不同,它可以用一个变量存储多个元素,如果你了解过其它语言,列表和数组是十分类似的。

一个列表是用一个方括号包围的,列表中的元素使用逗号分隔,例如下面这几个变量都是一个列表:

值得注意的是,一个列表可以不包含任何元素,就像变量 c 存储的列表一样,这样的列表叫做空列表

另外,一个列表内的元素不一定要有相同的类型,甚至列表本身也可以作为元素存储在一个列表中:

比如说上面第一个列表,它存储着三个不同类型的元素,一个是整数,一个是浮点数,另外一个是字符串。

至于第二个列表,它包含四个元素,第一个元素是一个整数 233 ,第二个元素是先前的列表 x ,第三个元素是一个空列表 [] ,最后一个元素是另外一个仅包含一个字符串作为元素的列表。

(more…)

Read More

CodeChef上一道有趣的网络流题目

这是一道2016年二月CodeChef上的一道题目,叫做Call Center Schedule

题目大意是这样的:

一个电话公司有 P 个员工,一个星期有 D 天,每天需要工作 H 个小时。每个员工在每个小时只能做参加会议和与客户通话两件事情中的一件,或者不做。

现在要为安排一个和客户通话的时间表,要求满足

  • 每个人每天最多花费 N 个小时在参加会议和通话上
  • 第 i 个人每周最多花费 Li 个小时和客户通话
  • 每个人每天在 LTbegin 到 LTend 的时间内至少要有 1 小时没有被安排任务
  • 对于第 i 天第 j 个小时,恰好有 Ri, j 个人可以和客户通话

另外,开会的时间表已经安排完成,如果 Fk, i, j 是 0 则表示第 k 个人在第 i 天的第 j 个小时开会,如果是 1 则表示他可以被安排去和客户通话。

时限是 2s,最多 5 组数据,并且 1 \leq N \leq H \leq 70, 1 \leq D, P \leq 70, 1 \leq L_i \leq ND, 0 \leq R_{i, j} \leq 15, F_{k, i, j} \in \{0, 1\}, 1 \leq LT_{begin} \leq LT_{end} \leq N

(more…)

Read More

Tarjan算法寻找有向图的强连通分量

在图论中,一个有向图被成为是强连通的(strongly connected)当且仅当每一对不相同结点 u 和 v 间既存在从 u 到 v 的路径也存在从 v 到 u 的路径。有向图的极大强连通子图(这里指点数极大)被称为强连通分量(strongly connected component)。

tarjan-graph

比如说这个有向图中,点 1, 2, 4, 5, 6, 7, 8 和相应边组成的子图就是一个强连通分量,另外点 3, 9 单独构成强连通分量。

Tarjan 算法是由 Robert Tarjan 提出的用于寻找有向图的强连通分量的算法。它可以在\mathcal O(|V|+|E|)的时间内得出结果。

Tarjan 算法主要是利用 DFS 来寻找强连通分量的。在介绍该算法之前,我们先来介绍一下搜索树。先前那个有向图的搜索树是这样的:

search-tree

有向图的搜索树主要有 4 种边(这张图只有 3 种),其中用实线画出来的是树边(tree edge),每次搜索找到一个还没有访问过的结点的时候就形成了一条树边。用长虚线画出来的是反祖边(back edge),也被叫做回边。用短虚线画出来的是横叉边(cross edge),它主要是在搜索的时候遇到了一个已经访问过的结点,但是这个结点并不是当前节点的祖先时形成的。除此之外,像从结点1到结点6这样的边叫做前向边(forward edge),它是在搜索的时候遇到子树中的结点的时候形成的。

(more…)

Read More

CodeChef February 2016. Weird Sum 题解

这是 CodChef February Challenge 2016 上的一题,题目链接在这里:WRDSUM

给定一个整数 n(n \geq 2),考虑其质因子分解 n = p_1^{k_1}p_2^{k_2}\cdots p_r^{k_r}

g = \text{gcd}(k_1, k_2, \cdots, k_r)m_i = \frac{k_i}{g}

定义函数 F(n) = p_1^{m_1}p_2^{m_2}\cdots p_r^{m_r}

现在给定一个整数 N(N \leq 10^{2016}),你需要计算

 S = F(2) + F(3) + \cdots + F(N)

答案对 998244353 取模

(more…)

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