
> 本文也是马克思·普朗克复杂物理研究所《非平衡态集体过程》第9讲的笔记。
> 

最近在蹭一门数学系的读书课，题目叫“非线性时间序列分析”。本来以为会是数学系的人抢物理系动力学方向的饭碗用的（可能实际上也确实是，或者更有可能是反过来），但是老师在前几讲一直把重点放在何种定理如何证明上面。数学系的嘛，这种对公理化系统的喜爱可以理解。但是发现我们几个上课的对这些证明都不太感兴趣之后，~~直接一个发卡弯漂移，教我们用 R 语言分析现实数据，现在新课还没上，不过可能是以股票数据做例子，我的老天鹅，这也太“经世济民”了吧……~~ 最后期末作业发现大家居然又都选择做 PPT 讲中心极限定理的证明，经世济民计划无疾而终。

以下是为了新课程，我自己提前做的准备。

## 下载、安装、环境配置（精神病版）

- [https://cloud.r-project.org/](https://cloud.r-project.org/)
- 看这一个链接就够了：[https://marketplace.visualstudio.com/items?itemName=Ikuyadeu.r](https://marketplace.visualstudio.com/items?itemName=Ikuyadeu.r)
- [https://github.com/randy3k/radian](https://github.com/randy3k/radian)
- [https://github.com/nx10/httpgd](https://github.com/nx10/httpgd)

`Fedora Linux` + `R` + `radian` + `vscode` + `R extension for vscode`

1. 在官网（[https://cloud.r-project.org/](https://cloud.r-project.org/)）下载并安装 R [解释器](https://python-interpreter)。
2. 在命令行输入 `R` 打开解释器，在 R 中输入 `install.packages("languageserver")` 
3. 在 vscode 的市场页面搜索 “R”，安装 R 语言插件。
4. 安装 radian：
    1. 在合适的 python [虚拟环境](https://pyhton-virtualenv) 里输入 `pip intall radian`
    2. 用 radian 取代 R：编辑 `~/.bashrc` , 在其中加入一句`alias r="radian"` ，重启命令行
    3. 在 vscode 的设置中找到 `R>Rterm: Linux`, 输入 radian 的安装路径（在命令行中输入 `which R` 可以找到）
5. 在 vscode 的市场页面搜索 “R debugger”，安装 `R debugger for vscode`
6. 在命令行输入 `R` 打开解释器，在 R 中输入 `install.packages("httpgd")`，安装可视化工具 `httpgd`

## 简单语法

来自马克思·普朗克复杂物理研究所《非平衡态集体过程》第9讲，这老师的讲课顺序简直了……

R 的语法也简直了……

### 简单数据结构 Some data types

```r
# vector
a <- c(1,2,3,NA)
b <- c("m","f","f","m")
# using a value
b[2]

# list
list.ab <- list(number=a, gender=b)
list.ab$number
list.ab[[1]]

# factors (categorial variable)
f <- factor(b)
f
# [1] m f f m f m m
# levels: f m
levels(f) <- c("female", "male")

# dataframe (lists of vectors of the same length)
x <- c(1,2,3)
y <- c("aa","bb","cc")
z <- c(TRUE, FALSE, TRUE)
df <- data.frame(first=x, second=y, third=z)
View(df)
# |      | first | second | third |  
# |    1 |     2 |     aa |  TRUE |
# |    2 |     3 |     bb | FALSE |
# |    3 |     5 |     cc |  TRUE |
```

### 控制流 Loops, conditional statements

```r
# for loop
for(i in 1:5) {
	print(i)
}
# while loop
while(!finished) {
	print("Hello")
}
# if statement
if(i<5) {
	print("Hello")
} else {
	print("Not Hello")
}
```

### 函数 Functions

```r
# calling a function
rnorm(5,mean=1,sd=1)

# defining a function
mysum <- function(a,b,c=1) { a + b + c }
mysum(1,1)
```

### 复杂数据类型 `data.table`

Some comparison:

- `data.table`: R package, fast and memory efficient
- `python.pandas`: python implementation of data frames
- `dplyr`: highly popular, easy to learn

```r
library(data.table)
dt <- as.data.table(df)
```

```r
# read and write
flights <- fread("path/to/your/flights.txt")
weather <- fread("path/tp/your/weather.txt")
```

Making data tidy can simplify the following analysis.

- Every column is an observable.
- Every row is an observation.

```r
# making data tidy
melt( dt,
      id.vars = "ID",
      value.name = "expression",
      variable.name = "cell"
)
# making data messy
dcast(dt, ID ~ cell+expression) # not understand
```

Get item

```r
d[i, j, by] # take `d`, subset rows using `i`, then **calculate** `j` grouped by `by`
# examples
planes[engines == 4]

# slice the 2 columns
planes[, .(tailnum, year)] 
# groupby and calculate
flights[, .(mean_delay = mean(dep_delay, na.rm=T)), by=carrier]
flights[time_hour>20, .(mean_delay=mean(dep_delay,na.rm=T)), by=.(month, origin)]
# math calculations
flights[, speed_kmh := 60*1.61*distance/air_time]
flights[, resc_distance := distance/mean(distance), by=carrier]
```

Merge 2 tables

```r
# left join
merge(a,b,all.x=T)
b[a] # same as above

# right join
merge(a,b,all.y=T)
a[b] # same as above

# inner join
merge(a,b, all=F)
a[b, nomatch=0]

# full join
merge(a,b,all=T)
```

Chaining operations

```r
# data.table way
weather[, ws_kmh:=1.61*wind_speed][, .(mean_ws=mean(ws_kmh)), by=month]
# use operator %>% to take the result on the left as the 1st argument on the right
library(magrittr)
weather[, ws_kmh:=1.61*wind_speed] %>%
.[, .(mean_ws = mean(ws_kmh)), by=month] %>%
head()
```


.r | R语言入门笔记


> 第十讲是由校外邀请来的嘉宾来讲数据降维，所以课程网站上没有课件。时间轴和语音转录相对来说不那么有意义了。笔记为英文，请点击页面右下角的语言选择英语。



.mpipks-note | 10. Dimension Reduction


上周三老板和艾萨鸽了，但是我们还是照常进行了讨论，然后讨论的主要内容很让我意外地落回了各种数学细节里，可能是本科生不在，我们都不好意思形而上的缘故吧。

## 疾病传播的朗之万方程——随机指数传播

两个要点：
1. 在有随机项的情况下，$$ d[log(x(t))] \neq (1/x)(dx/dt) $$
2. 前面 Ito formula 里面的 $$a$$, $$b$$ 和后面 $$a=\mu+\sigma\xi(t)$$ 里的 $$a$$ 容易弄混。

## 维纳过程 (Wiener process)

讨论过程中他们问我有什么问题，我当时正在走神（呃），于是随口问了一句：讲义第4页左侧的式子为什么等于零……老实把表达式展开就知道，是后面一项 $$\left<W_{t_{j+1}}-W_{t_{j}}\right>=0$$

讲义中维纳过程是作为第一讲中的高斯噪音的积分出现的，也就是：$$dW(t)=\xi(t)dt$$，其中高斯噪声满足以下两个条件：
- (1).  $$\left<\xi(t)\right>=0$$
- (2).  $$\left<\xi(t_1)\xi(t_2)\right>=\delta(t_2-t_1)$$

但是[维基百科中的介绍](https://zh.wikipedia.org/wiki/%E7%BB%B4%E7%BA%B3%E8%BF%87%E7%A8%8B)是直接定义了维纳过程本身，它应该是满足如下三个条件的随时间改变的随机变量 $$W(t)$$：
- (1). $$W(0)=0$$
- (2). 映射 $$t\mapsto W(t)$$ 在正实数轴上几乎处处连续
- (3). 对于两个不重叠的时间段 $$0 \leq s_1 < t_1 \leq s_2 < t_2$$, $$W(t_1-s_1)$$ 和 $$W(t_2-s_2)$$ 独立，$$W(t-s)\sim N(0,t-s)$$

这两种定义等价吗？好像是，没仔细想。不过维纳过程的存在是为了对布朗运动进行建模，布朗运动可以看作随机行走在步长趋近于0时的极限；在一年级学过的数学物理方法里，扩散方程的解也可以看作随机行走者经过时间t之后在不同位置出现的概率，也是取步长为0的极限。所以 $$\xi(t)$$ __似乎__ 可以导出 $$W(t)$$，甚至可以不要求 $$\xi(t)$$ 具有高斯分布。

以上似乎有类比谬误的嫌疑，而且假如思路正确的话，会引出一个新的问题——扩散方程的推导来自 Fick 定律，这个过程在何处引入了对随机变量在时间上的积分？

好吧，扯完了之后发现了这个链接：[http://physicallensonthecell.org/advanced-diffusion-and-fokker-planck-picture](http://physicallensonthecell.org/advanced-diffusion-and-fokker-planck-picture)

##  随机变量 vs. 时间变量

对于 $$\xi(t)$$ 函数的第二条性质，我在草纸上顺手写了如下等式：

$$ \left<\xi(t)\xi(t')\right>=\delta(t-t')=\frac{\int_{start}^{end}\int_{start}^{end}\xi(t)\xi(t')dtdt'}{\int_{start}^{end}\int_{start}^{end}dt\ dt'} $$

但是实际上，此处的 $$\left<\right>$$ 应该是随机变量的期望值，也就是

$$ \left<\xi(t)\xi(t')\right>=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}x_1 x_2 p(\xi(t)=x_1,\xi(t')=x_2)dx_1dx_2 $$

跟 t 没什么关系。在平衡态统计中，系综满足 ergodicity，两者应该是相等的。但在非平衡态下没有这种保证。

随机变量说是变量，实际上需要一个函数来描述，而这个函数的“自变量”是这个随机变量的所有可能取值，这个函数的输入和输出变量之间不是因果关系，这究竟是哥本哈根学派战胜隐变量理论的过人之处，还是新的隐变量理论可能的突破口呢？再往下胡扯存在被当成民科的可能，暂且打住。

## 主方程解得指数分布：随机变量的特征函数 (characteristic function in probability theory)

讲义第 9 页掉落了一件重要装备——解为指数分布的主方程：

$$ \frac{d}{dt}p(n,t) = \lambda \left[p(n-1,t)-p(n,t)\right] $$

上面说随机变量需要用一个函数来描述，但是究竟用哪个函数，却可以有多种选择。最直观的就是中学就学的，离散变量的概率分布函数和连续变量的概率分布密度函数。好处是意思直观，缺点是对于离散变量和连续变量需要分别讨论，两种函数的量纲不同，概率密度的量纲等于随机变量量纲的倒数。

为了统一量纲，所以有了累积分布函数，对离散和连续变量 $$X$$ 的定义都是一样的：$$F(x)=prob(X<x)$$。离散变量的概率分布函数是累积分布函数的差分 (difference)，连续随机变量的概率密度则是累积分布函数的微分 (differential)。

在本节第 9 页的例子里，给出了第三种函数——[特征函数](https://en.wikipedia.org/wiki/Characteristic_function_(probability_theory))，对于一个随机变量 $$X$$ 与其所有的可能取值 $$x$$，特征函数引入新的自变量 $$s$$，由 $$e^{isX}$$ 的期望值给出，也就是说，特征函数是对随机变量的一个 __变换__。

$$ G_X(s) = \left<e^{isX}\right> = \sum_x prob(X=x)e^{ixs}$$

因为指数函数的导数还是指数函数，所以在微分方程里很有用，下面也确实是这么用的。用变换解微分方程的套路，上一节已经遇到了，只不过研究的是动力学，总是拖着一个并未参战的 $$t$$，而且笔记里的形式参数和实际参数并不明确，要是在 C 语言里这么写，是会被测试打死的:
$$G(s,t)_{[笔记里的写法]}:=G_{N或者n[被变换的随机变量]}(s_{[特征函数真正的自变量]},t_{[时间，打酱油的]})$$

## 福克-普朗克方程 (Fokker-Planck equation)

第 10 页，Kramers-Moyal 展开的写法是：

$$ \partial_tp(x,t) = \sum_{n=1}^{\infty}\frac{(-1)^n}{n!}\partial_x^n\alpha_n(x)p(x,t) $$

维基百科上的写法是

$$ \partial_tp(x,t) = \sum_{n=1}^{\infty}\frac{(-1)^n}{n!}\partial_x^n\left[\alpha_n(x)p(x,t)\right] $$

从后一页的算式来看应该是后者，前者写法不规范。





.mpipks-note | 02. 朗之万方程和主方程


来来回回地协调了好几次时间之后，我们最终确定每周三下午进行一次讨论。刚刚过去的这一周讨论了第一讲，时长超过了一个小时。

在讨论会之前，我的主要精力放在了仅有的几张有公式的幻灯片上，结果这一部分完全没有出现在讨论当中。讨论的主要内容是组里本科生艾萨提出的问题：“细致平衡条件 (detailed balance) 和平衡 (equilibrium) 的关系是什么”。这个问题其实挺有质量的，而且我们确实也没有给出一个令人满意的回答。我觉得要想说清楚这个问题，其实需要的是逻辑学和科学哲学，我现在在这方面的水平还不够，这篇文章姑且尝试解释一下。不过不谈数学的原因恐怕不是在这方面没有疑问，毕竟一不用做作业，二没有考试，完全出于兴趣的情况下，直面数学是需要热情和勇气的。

## 朗之万方程和响应函数

### 朗之万方程

讲义第12页到15页，讲解了热力学平衡条件下的朗之万方程。所谓热力学平衡条件，真正用到的性质就只有能量均分定理。朗之万方程听起来高大上，实际上只不过是对牛顿第二定律进行了一下改写，明确指明了这个系统受到两种力的作用，下面的公式里，等号右边第一项是正比于速度的衰减项，防止方程的解在无限远处发散（"所以你承认你们科学家是在瞎编咯？"），另外一项是一个关于时间的随机函数，用来刻画不确定性的涨落。

$$ m \frac{d^2x(t)}{dt^2} = -\lambda \frac{dx(t)}{dt} + \xi(t) $$

朗之万是和爱因斯坦以及居里夫人（咳咳）同时代的人，依然可以用牛顿时代的东西青史标名，可见物理学果然还是内卷不足，广阔天地，大有可为。（逃）

### 粒子的位置是粒子之前受力作用的积累：响应函数是卷积的核

所谓牛顿第二定律，对于高中生来说是公式，对于大学生来说是方程，（“公式和方程有啥区别？”），对于刚上大学的高中毕业生来说是噩梦—— $$ F=ma $$  是方程？这玩意需要解？这玩意也能解？这玩意怎么解？三个字母倒真不算方程，但是物理受到的力往往取决于其所在的位置 $$ x(t) $$ ，所以F其实是 $$ F(x(t)) $$ ；$$ a $$  是加速度，位置的二阶导数，也就是 $$ d^2x(t)/dt^2 $$ ，这里的未知量是 $$ x(t) $$，这是一个微分方程，常微分方程往往是本科阶段要上一个学期的课程。

响应函数用的是另一套思路，先不从微分方程的角度入手，直接来猜猜解应该长什么样子。此刻的位置与一段时间以前的受力以及当时的那个力对此刻残留的影响有关，这个“一段时间”究竟是多长时间呢？不知道，索性考虑所有时间长度，把这些影响全都加起来。学过微积分应该知道，对无限个无穷小量求和，再取个极限，这就是积分（“快把那个要先证明积分存在并收敛的数学家嘴堵上”），写作：

$$ x(t) = \int_{-\infty}^{t}h(t-\tau)\xi(\tau)d\tau $$

熟悉数学的朋友知道，等号右边的积分叫做卷积。一般卷积的定义里形式稍有不同，积分上下限分别是正负无穷，这对物理学家不是问题，只要令  $$ \tau > t $$  的时候  $$ h(t-\tau) = 0 $$   就行了。

一方面，卷积既是一种定义在两个函数（$$ \xi $$ 和 $$ h $$）之间的运算，可以简写作

 $$ x = h * \xi $$  

另一方面，也可以看作是对 $$ \xi $$ 的变换 (transform)，  $$ h(t-\tau) $$ 是这个变换的核 (kernal)，也就是讲义中所讲的 response function。当一个具体问题确定的时候，方程的解就和响应函数一一对应。

在方程实际被求解之前，响应函数和待求解的函数一样，是一个未知量。因为牛顿第二定律是一个二阶微分方程，要想求解，所以很自然的思路就是将方程两边同时做傅里叶变换，将导数转化为含有 iw 的代数方程，然后再将解做傅里叶逆变换。因为两个函数的卷积的傅里叶变换，等于两个函数各自的傅里叶变换的乘积，所以在傅里叶倒空间里，方程的解、响应函数、随机涨落三项是很简单的乘积形式。

### 关于朗之万方程的解，讲义中不清楚的地方

第14页右侧蓝字似乎是 "spectral density".

第15页 $$ \tilde{C}(\omega) $$ 似乎是 $$ \left<\tilde v(\omega)^2\right> $$  

## 平衡 (equilibrium) 和细致平衡 (detailed balance)

先贴几个相关的链接：
- [维基百科：thermal equilibrium](https://en.wikipedia.org/wiki/Thermal_equilibrium)
- [维基百科：thermaldynamic equilibrium](https://en.wikipedia.org/wiki/Thermodynamic_equilibrium)
- [维基百科：detailed balance](https://en.wikipedia.org/wiki/Detailed_balance)
- [Physical Lens on the Cell: Equilibrium Means Detailed Balance](http://physicallensonthecell.org/chemical-physics/equilibrium-means-detailed-balance)

目前直接在谷歌搜索 "equilibrium" 并进入维基百科的话会发现找不到和物理有关的词条，因为 “mechanical equilibrium”, "thermal equilibrium", "thermaldynamic equilibrium" 是几个独立的词条，而且都没有链接到 "equilibrium" 的消歧义页面上。

所谓 “mechanical equilibrium”，就是受力平衡，分为三种：稳定平衡、不稳定平衡、随遇平衡。实际上就是牛顿第二定律作为一个动力学方程，其中的稳定、不稳定和中性三种驻点 (fixed point)，和这一讲的关系不大。

热平衡简单。**两个** 系统处于热平衡，那么两个系统的温度相同。**一个** 系统处于热平衡，指的是这个系统温度时间上恒定，空间上均一。虽然这样说有点鸡贼了，盗用了读者在日常生活中对“温度”的理解。

热力学平衡的条件严格一些，不止热平衡，还要达到机械平衡，化学平衡，所有宏观上的热力学物理量（温度、质心坐标、粒子数……）都没有变化。

据此，大概可以归纳物理中一般性的“平衡”概念：某种物理性质在一个系统中的平衡，指的是这种物理量在系统内时间上的恒定，空间上的均匀；在多个系统之间的平衡，指的是这一物理量在这些系统之间没有流动。

其实上面这一段的草稿里，“指的是这个系统温度时间上恒定，空间上均一”后面想当然地来了一句“也就是任意两个子系统之间相互平衡”。当你听到了两个系统之间的平衡之后，再想到一个系统可以划分为多个子系统，思维的惯性就会用类比的方法来如此定义一个系统的平衡。

但是其实仔细一想，“时间上恒定，空间上均一”并不一定要求任意两个子系统之间相互平衡，但是任意两个子系统互相平衡是一个很有趣的状态，我们给他起名叫“细致平衡”。当任意两个子系统互相平衡的时候，可以确定地达到时间上恒定空间上均一的结果。

**所以这是我对艾萨问题的回答：**细致平衡条件是平衡的充分不必要条件。之所以常有这样的疑惑，是因为平衡的定义不能递归，一个系统的平衡和多个系统之间的平衡定义不同。

上面最后一个链接可以说是直接回答了艾萨的问题，但是这个页面的主标题是 "Equilibrium Means Detailed Balance", 字面意思是“平衡意味着细致平衡”，也就是说平衡是细致平衡的充分条件，和我前面说的相反。但是文章的内容实际并没有讨论究竟是谁蕴含谁，所以姑且可以无视这个标题。

细致平衡条件在维基百科词条的正文中称作 "principle of detailed balance"，也就是说这是一个“原理”，和“定律”(law) 一样，这是科学 **理论** 这一逻辑体系的公理，是没办法在理论内部通过逻辑推导来证明或者证伪的。我们之所以采用这些原理（也就是认为他们是真命题），是因为现有的 **实验** 全都没有证明这些原理不正确。

多说两句，科学虽然存在理论，但归根结底是一种经验科学，必须落地于实验和观测。（这里的“经验”对应于“理性”，是有严格定义的哲学概念，不是我们日常用语中所说的那个经验。）而观测和实验总有认知的的边界，边界以外的事情科学家做不了保证，所以实验总没有办法一劳永逸地说某某是对的，只能说某某还没错。只要证而不伪，科学就说这是一个真命题。科学与非科学的界限，根据波普尔的科学哲学理论，就在于[可证伪性](https://program-think.blogspot.com/2015/10/What-is-Science.html)。

另外讨论课里完全没有涉及的一点是，这里的“子系统”必须依然是“宏观的”(macroscopic)，各个子系统只能是小一号的总系统，（用软件工程的话说，定义的接口是一样的），不涉及组成这个系统的微观机制（不能是另一类对象）。比如说一盒理想气体，子系统只能是给盒子里加上隔板后的几小盒理想气体，不能精确到具体的气体原子分子，否则此时连机械平衡都满足不了。

这和经典力学里的“宏观无穷小”这一扭曲的概念有异曲同工之妙。归根结底，是因为宏观的科学理论建立得不够早，已经用上了微积分；但又太早，没考虑也不知道更低层次的物理规律。经典力学的更低层次已经形成了量子力学，而热力学的底层统计物理就只能和热力学挤在一起，和理论力学、电动力学、量子力学并列为四大力学，逃不出“N 大某某有 (N+1) 个”的魔咒。


2021 年 5 月 15 日.mpipks-note | 10. Dimension Reductionmpipks-note

.mpipks-note | 10. Dimension Reduction

.mpipks-note马克思·普朗克研究所《非平衡态统计》笔记