# 坑 4 思维工具

## 4.3 统计思维

### 4.3.3 回归

# 代码源自animation包quincunx函数
op = par(mar = c(1, 0.1, 0.1, 0.1), mfcol = c(4, 1))
nmax = 10000+15-2
# 15层钉板
layers = 15
# 10000个小球
balls = 10000
# 模拟三层钉板
layerx = layery = newlayerx = newlayery = NULL
flagi = ifelse(layers%%2, 1, 0)
for (i in 1:layers) {
if (flagi) {
newi = ifelse(i%%2, 1, 2)
}
else {
newi = ifelse(i%%2, 2, 1)
}
layerx = c(layerx, seq(0.5 * (i + 1), layers - 0.5 *
(i - 1), 1))
layery = c(layery, rep(i, layers - i + 1))
if (i > 1) {
newlayerx = c(newlayerx, seq(0.5 * (newi + 1), layers -
0.5 * (newi - 1), 1))
newlayery = c(newlayery, rep(i, layers - newi + 1))
}
}
# 模拟小球过钉板
ballx = bally = matrix(nrow = balls, ncol = nmax)
finalx = newfinalx = numeric(balls)
# 第一层钉板
for (i in 1:balls) {
ballx[i, i] = (1 + layers)/2
if (layers > 2) {
tmp = rbinom(layers - 2, 1, 0.5) * 2 - 1
ballx[i, i + 1:(layers - 2)] = cumsum(tmp) * 0.5 +
(1 + layers)/2
}
bally[i, (i - 1) + 1:(layers - 1)] = (layers - 1):1
finalx[i] = ballx[i, i + layers - 2]
}
# 第二层钉板
rgx = c(1, layers)
rgy = c(0, max(table(finalx)))
newballx = ballx
diag(newballx) = finalx
for (i in 1:balls) {
tmp = rbinom(layers - 2, 1, 0.5) * 2 - 1
tmp = ifelse(newballx[i, i] + cumsum(tmp) < rgx[2], tmp,
-1)
tmp = ifelse(newballx[i, i] + cumsum(tmp) > rgx[1], tmp,
+1)
newballx[i, i + 1:(layers - 2)] = newballx[i, i] + cumsum(tmp) *
0.5
newfinalx[i] = newballx[i, i + layers - 2]
}
# 绘制钉板与分布图
plot(1:layers, type = "n", ann = FALSE, axes = FALSE)
points(layerx, layery, pch = 19)
hist(finalx,breaks = 1:layers,main = "", xlab = "", ylab = "", ann = FALSE, axes = FALSE, xlim = rgx, ylim = rgy)
plot(1:layers, type = "n", ann = FALSE, axes = FALSE)
points(newlayerx, newlayery, pch = 19)
hist(newfinalx,breaks = 1:layers,main = "", xlab = "", ylab = "", ann = FALSE, axes = FALSE, xlim = rgx, ylim = rgy)

### 4.3.4 因果

• A导致B，B导致C $$A\rightarrow B\rightarrow C$$
• A导致B跟C $$A\rightarrow B, A\rightarrow C$$
• A跟B共同导致C，控制C后A跟B相关 $$A\rightarrow C, B\rightarrow C, (C): A \rightarrow B$$
• 随机相关

$$A \rightarrow B \rightarrow C, D\rightarrow A, D\rightarrow C$$

$$A \rightarrow B \rightarrow C \leftarrow D \rightarrow E$$

$$遗传 \rightarrow 身高$$

$$吸烟 \rightarrow 出生体重 \rightarrow 婴儿死亡率$$

$$吸烟 \rightarrow (出生体重) \leftarrow 先天缺陷 \rightarrow 婴儿死亡率$$

$$A \rightarrow X \rightarrow B$$

$$A \rightarrow X\rightarrow B, C \rightarrow X, C \rightarrow B$$

$$A \rightarrow (X) \leftarrow C \rightarrow B$$

$ACE = E[Y{i,1}] −E[Y{i,0}]$

$$X \rightarrow A \rightarrow B, C \rightarrow A, C \rightarrow B$$

X 与 B 之间没有直接因果，在效应计算中，X 与 B 的相关只会来自于 A 的中介，也就是 A 与 X 的相关性与 A 与 B 的相关性乘积，如果我们打算估计 A 与 B 因果效应，就是用 B 与 X 的回归系数去除 A 与 X 的回归系数。这个过程中因为 C 跟 X 间没有因果通路，所以相当于跳过了混杂因素的讨论。这个路径影响等于回归系数乘积的代数技巧使得计算上非常清晰，但还是牢记工具变量的选择不是探索的而是因果的，需要背景知识。有了工具变量，研究人员就从无穷的相关与混杂中解放出来而可以通过因果图构建来求解目标因果关系。当因果图很复杂时，你需要去控制一些变量或引入工具变量实现你关心变量的 do 分离，进而估计净效应。

## 4.4 模型思维

### 4.4.1 编程

• 条件分支：函数中出现需要对数据子分类进行不同运算时的设计，不同子分类用不同条件语句进行逻辑判断，例如数值求绝对值要先判断正负。

• 循环：同样的操作要对不同的可索引或满足特定条件的数据进行运算，这种情况要设计循环结构，例如按数据行/列求值。有些循环循环数是知道的，有些则要对数据运行结果进行判断，满足特定条件时可跳出或继续循环。

• 递归：比较特殊的条件与循环结构，当数据不满足某条件时就执行函数本身直到满足条件，例如求解斐波那契数列之和就可以设计递归结构循环执行本身直到数据可计算的起点。递归的效率一般不高，但递归结构有助于简化思考问题的步骤。

• 正则表达式：正则表达式是字符串处理时常用的模式识别工具，灵活使用正则表达式与条件分支可以有效处理真实数据中的混杂，强烈推荐学习掌握。

• 数据结构：通常不同数据按照实际需求会有不同的格式，不同格式的数据处理方式会不一样，一般函数都会先验证数据结构，如果不能处理则返回错误。

• 数据表：常见的数据处理格式，一般不同行表示不同样品，不同列表示不同样品属性且数据类型一致，数据值可以是数值、字符或逻辑值但不能是数据表。由于数据处理算法大都基于数据表开发，这类格式数据比较容易找到现成的算法函数/库/扩展包来进行处理。

• 字典：很多程序语言支持字典，字典是一种对应关系，字典中的元素是键值-数值对，通过键值索引数值，也可以反查。数据表中搜索元素是按照数值索引顺序索引的，字典则可以用哈希表快速索引。字典可以在编程中用来构建基于输入的数据库，方便进一步查询。

• 列表：列表属于数据表与字典的泛化，列表元素可以是数据表或列表，因此列表的数据结构不是平行的而是具备层级，有的元素可以进一步展开。列表常用来表示一组关联概念且可以数值索引，例如在回归分析的返回值中，就会包括拟合值、回归系数、残差等数据表或数值。

• 类型：通常列表可被定义成一种新通用类型，算法可基于这个类型进行开发或泛化，例如当你调用画图程序时，其程序会首先判断你输入数据的类型，如果有对应方法则直接调用，没有则用通用方法或返回错误。有些语言中列表是不能直接操作的，这样设计就是为了防止类型不兼容而强制定义格式。

C的核心地位还体现在很多高级语言的编译器是构建在C之上的，或者是C++之上，多数高级语言都通过限制自由度（例如不让操作内存、功能模块化等）来实现上手容易与较高的开发效率，但只要关注程序性能，肯定回去学C或C++的。GNU的GCC编译器是一个相对通用的编译器合集，可以用来编译包括C在内的多种语言。然而，GCC也是有历史包袱的，所有有人就另起炉灶单独针对C或C++重写了效率更高的编译器及其后台，这就是苹果的LLVM项目与clang编译器。但要注意的是LLVM支持的语言不如GCC多，所以如果你还要用到fortain或java编译器，那就还是老老实实用GCC吧，或者cmake的时候分别指定编译器，只要你不嫌麻烦。一般而言，效率与性能往往不能兼得。

### 4.4.3 个体-整体模型

# 谢林隔离模型，代码源于 https://simulatingcomplexity.wordpress.com/2016/01/06/building-a-schelling-segregation-model-in-r/
# 生成单元格矩阵，两种颜色随机分布
number<-2000
group<-c(rep(0,(51*51)-number),rep(1,number/2),rep(2,number/2))
grid<-matrix(sample(group,2601,replace=F), ncol=51)
# 找邻居函数
get_neighbors<-function(coords) {
n<-c()
for (i in c(1:8)) {

if (i == 1) {
x<-coords[1] + 1
y<-coords[2]
}

if (i == 2) {
x<-coords[1] + 1
y<-coords[2] + 1
}

if (i == 3) {
x<-coords[1]
y<-coords[2] + 1
}

if (i == 4) {
x<-coords[1] - 1
y<-coords[2] + 1
}

if (i == 5) {
x<-coords[1] - 1
y<-coords[2]
}

if (i == 6) {
x<-coords[1] - 1
y<-coords[2] - 1
}

if (i == 7) {
x<-coords[1]
y<-coords[2] - 1
}

if (i == 8) {
x<-coords[1] + 1
y<-coords[2] - 1
}

if (x < 1) {
x<-51
}
if (x > 51) {
x<-1
}
if (y < 1) {
y<-51
}
if (y > 51) {
y<-1
}
n<-rbind(n,c(x,y))
}
n
}
par(mfrow=c(1,2))
image(grid,col=c("black","blue","yellow"),axes=F)
# 设定喜好度，此处设定为50%，也就是单元格里人周围少于四个单元格的人与自己颜色相同时，这个人就会随机搬家到另一个单元格
alike_preference<-0.5
# 模拟搬家十次
for (t in c(1:10)){
happy_cells<-c()
unhappy_cells<-c()
for (j in c(1:51)) {
for (k in c(1:51)) {
current<-c(j,k)
value<-grid[j,k]
if (value > 0) {
like_neighbors<-0
all_neighbors<-0
neighbors<-get_neighbors(current)
for (i in c(1:nrow(neighbors))){
x<-neighbors[i,1]
y<-neighbors[i,2]
if (grid[x,y] > 0) {
all_neighbors<-all_neighbors + 1
}
if (grid[x,y] == value) {
like_neighbors<-like_neighbors + 1
}
}
if (is.nan(like_neighbors / all_neighbors)==FALSE) {
if ((like_neighbors / all_neighbors) < alike_preference) {
unhappy_cells<-rbind(unhappy_cells,c(current[1],current[2]))
}
else {
happy_cells<-rbind(happy_cells,c(current[1],current[2]))
}
}

else {
happy_cells<-rbind(happy_cells,c(current[1],current[2]))
}
}
}
}
rand<-sample(nrow(unhappy_cells))
for (i in rand) {
mover<-unhappy_cells[i,]
mover_val<-grid[mover[1],mover[2]]
move_to<-c(sample(1:51,1),sample(1:51,1))
move_to_val<-grid[move_to[1],move_to[2]]
while (move_to_val > 0 ){
move_to<-c(sample(1:51,1),sample(1:51,1))
move_to_val<-grid[move_to[1],move_to[2]]
}
grid[mover[1],mover[2]]<-0
grid[move_to[1],move_to[2]]<-mover_val
}
}
# 观察此时的单元格颜色分布状况
image(grid,col=c("black","blue","yellow"),axes=F)