今天在知乎上看了一个问题，发现结果与自己起初的直觉大相径庭，同时回忆起了以前的一些统计知识。

问题

房间内有 100 人，每人有 100 块，每次随机给另一个人 1 块，最后这个房间内的财富分布怎样？

一开始直觉上告诉我，应该是平均100元左右，具体呈正态分布，但其实如果严谨一点的话，远远没有这么简单。

略微严谨一点的话，问题有两种可能（题目缺少条件）：

1. 允许财富为负值
2. 不允许财富为负值（即0块的人仅有可能收钱）

情况1进行10000次游戏

每人的最终财富状况如下

最终财富状况直方图如下

结果发现用户间的财富呈正态分布, 符合Centrel Limit Therom。
用户财富变化如下。（10人进行10000次游戏）

从该图可以直观的理解，如果所有人的投资理念一样（每次给钱和拿钱的概率相同），大概率穷人会更穷，而有钱人会更有钱。但是如果我们设定一个穷的阈值呢？比如余额不能是负数，这样余额为0元的人则无需给钱，但有可能拿钱。这样的情况下会呈现怎样的分布呢？这就是违反直觉的部分了

情况2进行10000次游戏

每人的最终财富状况如下

最终财富状况直方图如下

结果发现十分接近幂律分布（或Pareto Distribution），不由得联想起了与熵有关的知识（非专业，只是大概理解）：
幂律分布为平均值固定，熵最大的分布。
由于本游戏为一种零和博弈，故平均值固定，方差不断累积。
在知乎上也有看到其他通过随机游走的角度来解决问题，感兴趣的可以去找一下。

Reference
https://en.wikipedia.org/wiki/Pareto_distribution

library(tidyverse)
library(viridis)

give_money <- function(seq, n, allowNegative = FALSE){
  if(allowNegative == TRUE){
    if(seq[n] <= 0){
      return(seq)
      }
  }
  seq[n] <- seq[n] - 1
  index <- sample(1:length(seq), 1)
  seq[index] <- seq[index] + 1
  return(seq)
}

give_money_all_people <- function(seq, allowNegative = FALSE){
  for(i in 1:length(seq)){
    seq <- give_money(seq, i, allowNegative)
  }
  return(seq)
}

seq <- seq(100, by = 0, length.out = 100)

for(i in 1:10000){
  seq <- give_money_all_people(seq)
}

data <- cbind(data.frame(c(1:100)), data.frame(sort(seq)))
names(data) <- c('index', 'fortune')

ggplot(data = data, aes(x = index, y = fortune)) +
  geom_bar(stat='identity') +
  labs(x = 'Players #', y = 'Fortune') + 
  theme_bw()

ggplot(data = data, aes(x = fortune)) +
  geom_histogram() +
  labs(x = 'Fortune', y = 'Freq') + 
  theme_bw()


game_result <- function(seq, times, allowNegative = FALSE){
  data <- seq
  for(i in 1:times){
    seq <- give_money_all_people(seq, allowNegative)
    data <- c(data, seq)
  }
  mat <- matrix(data = data, nrow = times + 1, ncol = length(seq), byrow = TRUE)
  result <- data.frame(mat)
  times <- c(0:times)
  result <- cbind(data.frame(times), result)
  return(result)
}

seq <- seq(100, by = 0, length.out = 10)

data2 <- game_result(seq, 10000)

ggplot(data = data2) +
  geom_line(aes(x = times, y = X1), color = viridis(10)[1]) +
  geom_line(aes(x = times, y = X2), color = viridis(10)[2]) +
  geom_line(aes(x = times, y = X3), color = viridis(10)[3]) +
  geom_line(aes(x = times, y = X4), color = viridis(10)[4]) +
  geom_line(aes(x = times, y = X5), color = viridis(10)[5]) +
  geom_line(aes(x = times, y = X6), color = viridis(10)[6]) +
  geom_line(aes(x = times, y = X7), color = viridis(10)[7]) +
  geom_line(aes(x = times, y = X8), color = viridis(10)[8]) +
  geom_line(aes(x = times, y = X9), color = viridis(10)[9]) +
  geom_line(aes(x = times, y = X10), color = viridis(10)[10]) +
  labs(x = 'Times of Game', y = 'Fortune') +
  theme_bw()