Carmér‘s Theorem

我们知道,两个独立的正态随机变量加起来还是一个正态的随机变量。但是反过来呢?是不是两个独立的随机变量加起来是个正态,它们就都是正态呢?Carmér定理说的就是这样一件事情:假设X_1, \cdots, X_n是一组独立的随机变量,并且X_1+\cdots+X_n是服从正态分布的,那么X_1,\cdots,X_n都是服从正态分布的。

这是这学期概率论课上老师给的bonus问题,这个定理原始论文是“Über eine Eigenschaft der normalen Verteilungsfunktion”,一篇用德文写的。然后我看着它的公式大概猜了一下这篇文章在干啥…… 证明见下面正文,应该是对的吧。感觉这个证明把这学期实分析、复分析学的东西都用上了一些。

主要的思路是,先证明特征函数\mathbb E[e^{itX}]可以定义在整个复平面上,并且是一个解析函数。同时,证明这个函数的增长阶是2并且没有零点,之后用Hadamard因子分解定理说明它的形式。

还有一个叫做Raikov定理,和这个定理非常相似,把定理中的“正态分布”换成“Poisson分布”就是Raikov定理的内容。证明和这个定理也很像,不过需要先说明两个随机变量都是离散的。

(more…)

Read More

FPGA Console——硬件实现的VT220兼容终端

这是和@HarryChen在这学期数字逻辑课上做的项目。主要是利用FPGA模拟一个终端。项目放在https://github.com/miskcoo/fpga-virtual-console上。

现在大部分的操作系统都是内置有虚拟终端,就像Linux下的xterm、gnome-terminal。但是,在最初控制台和计算机是相互独立的。控制台一般是由键盘和屏幕组成,键盘将用户输入传入计算机,计算机将输出发送显示给用户。这整个通信的过程是有一套标准来规定的,例如VT100就是一个早期的标准,更进一步有VT220、xterm-256color等。由于计算机回显的数据并不是单纯地整个屏幕的像素数据,而是更高层次的类似控制命令的内容,例如表示光标移动到某个位置,在某个位置插入某个字符,删除某一行的字符等。这些标准主要是定义了计算机发送回的这些控制命令的行为,和键盘发送给计算机的控制命令的行为。

我们做的这个项目就是一个物理上的终端!通过FPGA接受键盘的输入,将输入转化为控制命令通过串口输出给计算机。同时也通过串口接受计算机传回的控制命令名且解析、执行,修改对应位置的字符,再将字符进行渲染通过VGA输出到屏幕。

主要支持的标准是VT220,以及xterm-256color下大部分影响最终显示效果的指令,对于大部分linux下的程序的输出都可以正确地解析。同时串口速度提升到了3M可以支持以大约170fps的帧率播放黑白字符视屏,和约30fps的帧率播放彩色字符视屏。简单来说,就是正常终端可以做的事情这个物理终端都可以做到。

比如下面就是一个调用笔记本上摄像头的例子(我真的不是直接把VGA连到电脑上了),这时候串口速度还没那么快所以看起来有点卡……

fpga-camera

具体的细节我就不打算说了,文档里都写得非常明白。

https://github.com/miskcoo/fpga-virtual-console/raw/v1.0/doc/report.pdf

Read More

Pólya’s Urn

Pólya’s Urn是这样一个模型:假设现在有一个盒子,里面有r个红球和g个绿球,现在不断地随机从盒子取出一个球,之后再放回去,外加放入额外的c个和取出的球颜色相同的球。现在问,盒子中绿色球的比例在最后是否会趋向与某个分布?

实际上最后这个比例会是一个参数是g/cr/c的Beta分布。

首先我们设S_n = r + g + nc是第n次操作后盒子中所有球的数目。再设G_n, R_n分别是第n次操作后盒子中绿色球、红色球的数目。以及令X_n = G_n / S_n, \mathcal F_n = \sigma(G_k, S_k, 0 \leq k \leq n)。计算条件期望

 \begin{aligned}
\mathbb E[G_{n+1} | \mathcal F_n] &= \frac{G_n}{S_n}(G_n + c) + \frac{R_n}{S_n}G_n \\
&= (c + S_n) X_n = S_{n+1} X_n
\end{aligned}

可以看出\mathbb E[X_{n+1}|\mathcal F_n] = X_n\{X_n\}是一个鞅。再由于0\leq X_n \leq 1,我们知道X_n几乎处处收敛到一个随机变量X。现在考虑其分布。

首先可以计算出在第n步及以前,有m \leq n次取出绿色球的概率,并且写为紧凑的形式如下

 \begin{aligned}
&\mathbb P\left[ X_n = \frac{mc + g}{nc + g + r} \right] \\&= \frac{n!}{m!(n-m)!}\cdot \frac{\displaystyle\prod_{k=0}^{m-1}(k + g/c) \prod_{k=0}^{n - m - 1}(k + r/c)}{\displaystyle\prod_{k=0}^{n - 1}(\frac{g+r}{c} + k)}\\&= \frac{\Gamma(n + 1)}{\Gamma(m + 1)\Gamma(n - m + 1)}\cdot \frac{\Gamma(m + g/c)\Gamma(n - m + r/c)}{\Gamma(n + (g + r)/c)\text{B}(r/c, g/c)}
\end{aligned}

其中\Gamma\text{B}分别是Gamma函数和Beta函数。

接下来利用Stirling公式\Gamma(x) = \sqrt{2\pi/x} (x/e)^x (1 + O(1/x)),可以得到

 \frac{\Gamma(x + a)}{\Gamma(x + b)} = x^{a - b} \left( 1 + O\left( \frac{1}{x} \right) \right)\text{ as } x \to \infty

那么对于固定的\alpha, \beta, m满足0 < \alpha < \beta < 1以及\alpha n < m < \beta n,可以有

 \begin{aligned}
&\mathbb P\left[ X_n = \frac{mc + g}{nc + g + r} \right] \\&= \frac{1}{\text{B}(r/c, g/c)} n^{1 - (g + r)/c} (n - m)^{r/c - 1} m^{g/c - 1} \left( 1 + O\left( \frac{1}{n} \right) \right)\\
&= \frac{1}{n \text{B}(r/c, g/c)} \left (1 - \frac{m}{n} \right )^{r/c - 1} \left( \frac{m}{n} \right)^{g/c - 1} \left( 1 + O\left( \frac{1}{n} \right) \right)
\end{aligned}

对于所有满足\alpha n < m < \beta nm求和

 \lim_{n\to \infty} \mathbb P[\alpha \leq X_n \leq \beta] = \frac{1}{\text{B}(r/c, g/c)}\int_{\alpha}^\beta (1 - x)^{r/c - 1}x^{g/c - 1} dx

最后可以得到X的分布是\text{Beta}(g/c, r/c)

Read More

一道有意思的题目

前一段时间在做题的时候看到一道比较有意思的题目,大意是这样的

在单位圆周C上任意给出n个点w_1, \cdots, w_n,要求证明在这个圆周C上存在一个点w使得它到前面n个点距离的乘积为1。也就是

 \prod_{i=1}^n | w - w_i | = 1

这道题目是从 Stein 的 Complex Analysis 上来的,一眼看过去似乎和复分析没有任何关系… 但是却可以利用极大模原理来证明w的存在性。

大概是首先将w_1, \cdots, w_n放在复平面上,看成一个复数。然后定义函数

 f(z) = \prod_{i=1}^n (z - w_i)

这个函数可以定义在整个\mathbb C上并且是一个整函数。考虑其在z = 0处的模|f(0)| = \prod_{i=1}^n |w_i| = 1,那么按照极大模原理,f(z)在开单位圆\mathbb D内部不能达到极大模(除非它是常数),因此边界C = \partial \mathbb D上一定存在一个点z_0使得|f(z_0)| \geq 1。之后再利用|f(z)|的连续性,以及|f(w_1)|=1和介值定理就可以说明w的存在性。

证明还是挺有意思的,我想了一会儿似乎没想到直接通过几何的办法该怎么证…

Read More

平面等距变换、三维空间的旋转和四元数

之前寒假从Artin的Algebra上看到了一些有意思的东西,主要是关于对称和旋转这样的变换的。本来很早就想写了但是拖到现在,还有一部分内容是关于平面无限密铺的,估计这学期要被两门分析支配了就以后再说了。有兴趣可以自行阅读Artin的书。

首先可以考虑一个这样的问题,将平面上一点绕着某点P旋转\theta角后再绕Q旋转\varphi角,最终得到的变换是什么?能否表示成为绕另外某个点R旋转某个特定角度的形式?

直接从几何上考虑这个问题似乎比较困难,通过分析平面上的等距变换(即由旋转、平移和对称组成的变换)的结构,可以让我们从代数角度对这样的变换进行化简,将几何问题转化为解线性方程组的问题。

另外,三维空间中的旋转可以通过旋转矩阵来进行描述,但是旋转矩阵的结构较为复杂。我们通过分析旋转的特性,以及二维的酉矩阵P \in SU_2的性质,推导了四元数表示旋转的方法。如果对四元数感兴趣可以直接跳到对应部分,这篇对四元数的介绍的风格可能比较奇特。

对于文中的部分内容,需要基本的线性代数的知识,特别是需要知道特征值等方面的知识。对于证明直接略过是并不影响阅读的。同时,平面等距变换和三维空间中的旋转是两个相对独立的部分,可以分别阅读。

(more…)

Read More

实序列离散傅里叶变换的快速算法

\newcommand{\dft}[1]{\text{DFT}\left[#1\right]}从多项式乘法到快速傅里叶变换中我们介绍了快速傅立叶变换的算法。由于在计算多项式乘法或者大整数乘法时,我们所使用的都是实数,而快速傅立叶变换并不仅适用于实序列,对于复序列也是可以计算的,这让我们考虑是否可以利用虚部来保存额外的信息而减少计算的次数。

事实上,对于实序列x(n),其DFT是共轭对称的,也就是X^*(n) = X(N - n),如果将两个实数序列x(n)y(n)合并称为新的z(n) = x(n) + i\cdot y(n),对其进行FFT,我们可以从中复原出分别对两个序列进行FFT的结果。也就是,可以通过两次FFT就进行多项式乘法的计算(原先是需要三次的)。

沿用其中的记号,先回忆一下离散傅里叶变换的公式。

考虑长度为N的序列x(n)及其N点离散傅里叶变换X(n) := \dft{x(n)}

\begin{equation} X(n) = \sum_{k=0}^{N - 1} x(k)\cdot\omega_N^{kn},~ n = 0, 1, \cdots, N - 1 \end{equation}

其中\omega_N = e^{-2\pi i/N},即N次单位根。

现在来考虑X(n)X(N - n)之间的关系,我们用X^*(n)\overline{X(n)}来表示X(n)的共轭

\begin{equation} X^*(n) = \sum_{k=0}^{N - 1} x^*(k)\cdot\omega_N^{-kn} = \sum_{k=0}^{N - 1} x^*(k)\cdot\omega_N^{k(N-n)} \end{equation}

这是因为 \overline{\omega_N^{kn}} = e^{2\pi ikn/N} = e^{-2\pi ik(N - n)/N} = \omega_N^{k(N - n)}

那么,对于实序列,x^*(n) = x(n),我们可以得到 X^*(n) = X(N - n)。也就是说,实序列的DFT是共轭对称的。

现在来考虑两个实序列x(n)y(n),其对应的DFT为X(n)Y(n),定义新的序列z(n) = x(n) + i\cdot y(n),其对应DFT为Z(n),那么

\begin{equation} Z(n) = \dft{x(n) + i\cdot y(n)} = X(n) + i\cdot Y(n) \end{equation}

注意这里X(n)Y(n)本身并不是一个实数,因此它们也不是Z(n)的实部和虚部,我们设Z(n)的实部和虚部分别为A(n)B(n)

\begin{equation} \begin{aligned} \overline{X(N - n) + i\cdot Y(N - n)} &= \overline{A(N - n) + i\cdot B(N - n)} \\ \Rightarrow X(n) - i\cdot Y(n) &= A(N - n) - i\cdot B(N - n) \end{aligned} \end{equation}

这样就可以得到

\begin{equation} \begin{aligned} 2X(n) &= [A(n) + A(N - n)] + i[B(n) - B(N - n)] \\ 2Y(n) &= [B(n) + B(N - n)] - i[A(n) - A(N - n)] \end{aligned} \end{equation}

这样,对于两个实序列的DFT可以只利用一次FFT以及一点额外的运算进行计算而不必分别计算两次FFT。

下面来考虑一个长度为 2N 的实序列 x(n)。我们将其按照下标的奇偶分解为两个长度为N的序列

\begin{equation} \begin{aligned} x_1(n) &= x(2n)\\ x_2(n) &= x(2n + 1) \end{aligned} \end{equation}

现在可以利用刚刚的方法通过一次N点的FFT来计算它们分别的DFT序列X_1(n)X_2(n),通过类似于FFT的办法将其合并

\begin{equation} \begin{aligned} X_1(n) + \omega_{2N}^n X_2(n) &= \sum_{k=0}^{N-1}x_1(k)\omega_N^{kn}+\omega_{2N}^n\sum_{k=0}^{N-1}x_2(k)\omega_N^{kn} \\ &= \sum_{k=0}^{N-1}x(2k)\omega_{2N}^{2kn}+\sum_{k=0}^{N-1}x(2k + 1)\omega_{2N}^{(2k+1)n} \\ &= X(n) \end{aligned} \end{equation}

这里只有前N项,对于其余项可以利用共轭对称这个性质直接计算。 对于大整数乘法,我们也可以利用两个实序列DFT的快速算法进行计算,这样只需要两次FFT即可,下面我提供一个多项式乘法的代码(UOJ#34),FFT的模板是很久以前写的就别看了,可以看看下面处理的部分。

(more…)

Read More

程序设计训练 - 数独、国际跳棋、人物信息检索系统

这次程序设计训练小学期一共有三个大作业,分别是数独、在线国际跳棋和人物信息检索系统。这里是简单的一些介绍,代码放在 https://github.com/miskcoo/programming-training

Week 1: Sudoku Game

简介

Sudoku 是一款利用 Qt 实现的数独游戏,提供了多达 10 个难度的关卡选择,同时还有丰富的功能来帮助玩家更加高效地求解数独问题,例如候选数、高亮相同数字、高亮选中的行列、撤销当前操作以及提示等功能。玩家还可以手动输入数独题目利用 Sudoku 帮助求解。

除了传统 9x9 的数独游戏以外,还提供了更高难度的 16x16 的数独游戏。

基本功能

Sudoku提供了多个方便的按钮:

  • 新游戏:玩家可以开始一局新的游戏。
  • 重玩:玩家可以重新开始本局游戏。
  • 暂停:玩家可以暂停该局游戏(即暂停计时)。
  • 提示:如果当前已经确定的数都是正确的,玩家将会得到一个未填空格的正确数字;如果当前已经确定的数和答案矛盾,导致整个数独无解,那么所有与答案矛盾的数字将会被粗体标出。
  • 清除:清除当前选中格子的所有数字。
  • 撤销:撤销前一步的操作,以及取消撤销(最多可支持 50 步撤销)。

同时可以通过菜单来实现多达 10 种难度的游戏选择,可以求解任意用户输入的数独问题。

玩家在空格中填入的数字分为确定数(采用正常大小的字体表示)以及候选数(采用小号字体表示)。确定数只能存在一个,候选数可以存在多个,并且确定数和候选数不能同时存在。

玩家可以通过多种方式进行格子的选择:

  • 鼠标直接点击格子进行选中。
  • 键盘方向键来进行当前选中格子的切换。
  • Tab键快速切换到下一个格子。

玩家可以用鼠标右键点击格子来对其进行标记。

玩家在选中一个格子之后,如果该格子的数字已经确定,那么所有与该数相同的数将会被加粗显示,以方便确认是否满足数独的条件。同时,所有与该格在相同行或列的格子都会被高亮显示。此外,无论该格子数字确定与否,右侧数字列表中该格子的所有数都会被加粗。

玩家可以通过多种方式在空格中填入数字:

  • 键盘直接输入数字将填入对应的确定数,如果 Ctrl 键被按下,那么填入的将会是候选数。
  • 鼠标点击右侧数字列表填入。使用鼠标左键点击将会填入确定数,右键点击则会填入候选数。

玩家也可以通过多种方式在空格中填入数字:

  • 键盘删除键删除格子中的一个数字。
  • 鼠标点击右侧数字列表已经选中的数进行删除。
  • “清除”功能键来清除该格子所有的数字。

(more…)

Read More

多维正态分布中的边际分布、条件分布及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就会自动变成如下

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

我比较喜欢的是类似 \alpha 的命令可以用 `a 来快速输入,还有就是 \frac{}{} 这样的命令可以用 `/ 或者输入 frac + <F7> 来补全在输入数学公式的时候特别方便。

具体请戳: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