shadowsocks服务端以及客户端配置

发表于 2019-03-04 | 更新于 2019-12-17 | 分类于 linux

服务器端配置

首先需要有一个VPS账号，vultr,digitalocean,搬瓦工等等都行。
首先到下面两个网站检测22端口是否开启，如果关闭的话，vps换个ip把。。
http://tool.chinaz.com/port
https://www.yougetsignal.com/tools/open-ports/

启用BBR加速

~#:apt update
~#:apt upgrade
~#:echo “net.core.default_qdisc=fq” >> /etc/sysctl.conf
~#:echo “net.ipv4.tcp_congestion_control=bbr” >> /etc/sysctl.conf
~#:sysctl -p
上述命令就完成了BBR加速，执行以下命令验证：
~#:lsmod |grep bbr
看到输出包含tcp_bbr就说明已经成功了。

搭建shadowsocks server

安装shadowsocks server

~#:apt install python-pip
~#:pip install shadowsocks
需要说一下的是，shadowsocks目前还不支持python3.5及以上版本，上次我把/usr/bin/python指向了python3.6，就是系统默认的python指向了python3.6，然后就gg了。一定要使用Python 2.6,2.7,3.3,3.4中的一个版本才能使用。。

创建shadowsocks配置文件

如果你的VPS支持ipv6的话，那么可以开多进程分别运行ipv4和ipv6的shadowsocks server。本地只有ipv4的话，可以用本地ipv4访问ipv6，从而访问byr等网站，但是六维空间对此做了屏蔽。如果本地有ipv6的话，还可以用本地的ipv6访问ipv6实现校园网不走ipv4流量。

ipv4配置

~#:vim /etc/shadowsocks_v4.json
配置文件如下

{
"server":"0.0.0.0",
"server_port":"你的端口号",
"local_address":"127.0.0.1",
"local_port":1080,
"password":"你的密码",
"timeout":600,
"method":"aes-256-cfb"
}

ipv6配置

~#:vim /etc/shadowsocks_v6.json
配置文件如下

{
"server":"::",
"server_port":"你的端口号",
"local_address":"127.0.0.1",
"local_port":1080,
"password":"你的密码",
"timeout":600,
"method":"aes-256-cfb"
}
``` 
注意这两个文件的server_port一定要不同，以及双引号必须是英文引号。
##### 1.2.2.3.手动运行shadowsocks server
~#:ssserver -c /etc/shadowsock_v4.json -d start --pid-file ss1.pid
~#:ssserver -c /etc/shadowsock_v6.json -d start --pid-file ss2.pid
注意这里要给两条命令分配不同的进程号。

### 设置shadowsocks server开机自启
如果重启服务器的话，就需要重新手动执行上述命令，这里我们可以把它写成开机自启脚本。
~#:vim /etc/init.d/shadowsocks_v4
内容如下：
``` shell
#!/bin/sh
### BEGIN INIT INFO
# Provides:          apache2
# Required-Start:    $local_fs $remote_fs $network $syslog
# Required-Stop:     $local_fs $remote_fs $network $syslog
# Default-Start:     2 3 4 5
# Default-Stop:      0 1 6
# Short-Description: apache2 service
# Description:       apache2 service daemon
### END INIT INFO
start(){
  ssserver -c /etc/shadowsocks_v4.json -d start --pid-file ss2.pid
}
stop(){
  ssserver -c /etc/shadowsocks_v4.json -d stop --pid-file ss2.pid
}
case "$1" in
start)
  start
  ;;
stop)
  stop
  ;;
restart)
  stop
  start
  ;;
*)
  echo "Uasage: $0 {start|reload|stop}$"
  exit 1
  ;;
esac

~#:vim /etc/init.d/shadowsocks_v6
内容如下：

#!/bin/sh
### BEGIN INIT INFO
# Provides:          apache2
# Required-Start:    $local_fs $remote_fs $network $syslog
# Required-Stop:     $local_fs $remote_fs $network $syslog
# Default-Start:     2 3 4 5
# Default-Stop:      0 1 6
# Short-Description: apache2 service
# Description:       apache2 service daemon
### END INIT INFO
start(){
  ssserver -c /etc/shadowsocks_v6.json -d start --pid-file ss1.pid
}
stop(){
  ssserver -c /etc/shadowsocks_v6.json -d stop --pid-file ss1.pid

}
case "$1" in
start)
  start
  ;;
stop)
  stop
  ;;
restart)
  stop
  start
  ;;
*)
  echo "Uasage: $0 {start|reload|stop}$"
  exit 1
  ;;
esac

然后执行下列命令即可：
~#:chmod a+x /etc/init.d/shadowsocks_v4
~#:chmod a+x /etc/init.d/shadowsocks_v6
~#:update-rc.d shadowsocks_v4 defaults
~#:update-rc.d shadowsocks_v6 defaults

至此，服务器端配置完成。

服务端自动配置脚本

地址
首先将该文件中所有文件复制到vps上，然后执行
~#:sh install_ss_server.sh
即可

补充说明

该文件夹共包含五个文件
shadowsocks_v4.json为ipv4 ss配置文件，可根据自己的需要修改端口号和密码
shadowsocks_v6.json为ipv6 ss配置文件，可根据自己的需要修改端口号和密码
shadowsocks_v4为ipv4 ss自启动文件，无需修改
shadowsocks_v6为ipv6 ss自启动文件，无需修改
install_ss_server.sh为安装脚本，该脚本同时配置ipv4和ipv6 ss server。可根据自己需要自行选择。

客户端配置

Windows客户端配置

安装shadowsock客户端

到该网址 https://github.com/shadowsocks/shadowsocks-windows/releases 下载相应的windows客户端程序。
然后配置服务器即可～

Linux客户端配置

安装shadowsocks程序

~$:sudo pip install shadowsocks

运行shadowsocks客户端程序

~$:sudo vim /etc/shadowsocks.json
填入以下配置文件
{
“server”:“填上自己的shadowsocks server ip地址”,
“server_port”:“8888”,//填上自己的shadowsocks server 端口"
“local_port”:1080,
“password”:“mxxhcm150929”,
“timeout”:600,
“method”:“aes-256-cfb”
}

接下来可以执行以下命令运行shadowsocks客户端：
~$:sudo sslocal -c /etc/shadowsocks.json
然后报错：

INFO: loading config from /etc/shadowsocks.json
2019-03-04 14:37:49 INFO loading libcrypto from libcrypto.so.1.1
Traceback (most recent call last):
File “/usr/local/bin/sslocal”, line 11, in
load_entry_point(‘shadowsocks==2.8.2’, ‘console_scripts’, ‘sslocal’)()
File “/usr/local/lib/python2.7/dist-packages/shadowsocks/local.py”, line 39, in main
config = shell.get_config(True)
File “/usr/local/lib/python2.7/dist-packages/shadowsocks/shell.py”, line 262, in get_config
check_config(config, is_local)
File “/usr/local/lib/python2.7/dist-packages/shadowsocks/shell.py”, line 124, in check_config
encrypt.try_cipher(config[‘password’], config[‘method’])
File “/usr/local/lib/python2.7/dist-packages/shadowsocks/encrypt.py”, line 44, in try_cipher
Encryptor(key, method)
File “/usr/local/lib/python2.7/dist-packages/shadowsocks/encrypt.py”, line 83, in init
random_string(self._method_info[1]))
File “/usr/local/lib/python2.7/dist-packages/shadowsocks/encrypt.py”, line 109, in get_cipher
return m[2](method, key, iv, op)
File “/usr/local/lib/python2.7/dist-packages/shadowsocks/crypto/openssl.py”, line 76, in init
load_openssl()
File “/usr/local/lib/python2.7/dist-packages/shadowsocks/crypto/openssl.py”, line 52, in load_openssl
libcrypto.EVP_CIPHER_CTX_cleanup.argtypes = (c_void_p,)
File “/usr/lib/python2.7/ctypes/init.py”, line 379, in getattr
func = self.getitem(name)
File “/usr/lib/python2.7/ctypes/init.py”, line 384, in getitem
func = self._FuncPtr((name_or_ordinal, self))
AttributeError: /usr/lib/x86_64-linux-gnu/libcrypto.so.1.1: undefined symbol: EVP_CIPHER_CTX_cleanup

按照参考文献4的做法，是在openssl 1.1.0版本中放弃了EVP_CIPHER_CTX_cleanup函数

EVP_CIPHER_CTX was made opaque in OpenSSL 1.1.0. As a result, EVP_CIPHER_CTX_reset() appeared and EVP_CIPHER_CTX_cleanup() disappeared.
EVP_CIPHER_CTX_init() remains as an alias for EVP_CIPHER_CTX_reset().

将openssl库中的EVP_CIPHER_CTX_cleanup改为EVP_CIPHER_CTX_reset即可。
再次执行以下命令，查看shadowsocks安装位置
~#:pip install shadowsocks
Requirement already satisfied: shadowsocks in /usr/local/lib/python2.7/dist-packages
~#:cd /usr/local/lib/python2.7/dist-packages/shadowsocks
~#:vim crypto/openssl.py
搜索cleanup，将其替换为reset
具体位置在第52行libcrypto.EVP_CIPHER_CTX_cleanup.argtypes = (c_void_p,)和第111行libcrypto.EVP_CIPHER_CTX_cleanup(self._ctx)

手动运行后台挂起

将所有的log重定向到~/.log/sslocal.log文件中
~$:mkdir ~/.log
~$:touch ~/.log/ss-local.log
~$:nohup sslocal -c /etc/shadowsocks_v6.json </dev/null &>>~/.log/ss-local.log &

开机自启shadowsocks client

但是这样子的话，每次开机都要重新运行上述命令，太麻烦了。可以写个开机自启脚本。执行以下命令：
~$:sudo vim /etc/init.d/shadowsocks
内容为以下shell脚本

#!/bin/sh

### BEGIN INIT INFO
# Provides:          shadowsocks local
# Required-Start:    $local_fs $remote_fs $network $syslog
# Required-Stop:     $local_fs $remote_fs $network $syslog
# Default-Start:     2 3 4 5
# Default-Stop:      0 1 6
# Short-Description: shadowsocks service
# Description:       shadowsocks service daemon
### END INIT INFO
start(){
　　  sslocal -c /etc/shadowsocks.json -d start
}
stop(){
　　  sslocal -c /etc/shadowsocks.json -d stop
}
case “$1” in
start)
　　　start
　　　;;
stop)
　　　stop
　　　;;
reload)
　　　stop
　　　start
　　　;;
\*)
　　　echo “Usage: $0 {start|reload|stop}”
　　　exit 1
　　　;;
esac

然后执行以下命令即可：
~$:sudo chomod a+x /etc/init.d/shadowsocks
~$:sudo update_rc.d shadowsocks defaults
上述命令执行完成以后，进行测试
~$:sudo service shadosowcks start

配置代理

上一步的目的是建立了shadowsocks服务的本地客户端，socks5流量会走该通道，但是浏览器的网页的流量是https的，我们需要配置相应的代理，将https流量转换为socks5流量，走ss客户端到达ss服务端。当然，也可以把其他各种流量，如tcp,udp等各种流量都转换为socks5流量，这个可以通过全局代理实现，也可以通过添加特定的代理规则实现。

配置全局代理

如下图所示，添加ubuntu socks5系统代理：

然后就可以成功上网了。

使用SwitchyOmega配置chrome代理

首先到 https://github.com/FelisCatus/SwitchyOmega/releases 下载SyitchyOmega.crx。然后在chrome的地址栏输入chrome://extensions，将刚才下载的插件拖进去。
然后在浏览器右上角就有了这个插件，接下来配置插件。如下图：
mxx
直接配置proxy，添加如图所示的规则，这样chrome打开的所有网站都是走代理的。

使用privoxy让terminal走socks5

~$:sudo apt install privoxy
~$:sudo vim /etc/privoxy/config
取消下列行的注释，或者添加相应条目
forward-socks5 / 127.0.0.1:1080 . # SOCKS5代理地址
listen-address 127.0.0.1:8118 # HTTP代理地址
forward 10.*.*.*/ . # 内网地址不走代理
forward .abc.com/ . # 指定域名不走代理
重启privoxy服务
~$:sudo service privoxy restart
在bashrc中添加如下环境变量
export http_proxy="http://127.0.0.1:8118"
export https_proxy=“http://127.0.0.1:8118”

~$:source ~/.bashrc
~$:curl.gs

参考文献

DQN

发表于 2019-03-02 | 更新于 2019-12-17 | 分类于强化学习

背景

Atari 2600是一个RL benchmark，有2600个游戏，每个agent会得到一个图像输入(60Hz的210 x 160 RGB视频)。本文的目标是设计一个NN架构尽可能学会更多游戏，网络的输入只有视频信息，reward和terminal信号以及可能采取的action，和人类玩游戏时得到的信息是一样的。
Agent与Atari模拟器不断交互，agent不能观测到模拟器的内部状态，只能得到当前屏幕信息的一个图片。这个task可以认为是部分可观测的，因为仅仅从当前的屏幕图像$x_t$上是不能完全理解整个游戏状况的。所有的序列都认为在有限步骤内是会结束的。
注意agent当前的得分取决于整个sequence的action和observation。一个action的feedback可能等到好几千个timesteps之后才能得到。
agent的目标最大化累计reward。定义$t$时刻的回报return为$R_t = \sum^T_{t’=t} \gamma^{t’-t}r_{t’}$，其中$\gamma$是折扣因子，$T$是游戏终止的时间步。
定义最优的动作值函数$Q^{*}(s,a)$是遵循最优策略在状态$s$处采取动作$a$能获得的最大的期望回报，$Q^{*}(s,a) = \max_{\pi}E[R_t|s_t=s,a_t=a,\pi]$。
最优的动作值函数遵循Bellman optimal equation。如果在下个时间步的状态$s’$处，对于所有可能的$a’$，$Q^{*}(s’,a’)$的最优值是已知的（这里就是对于每一个$a’$，都会有一个最优的$Q(s’,a’)$，最优的策略就是选择最大化$r+Q^{*}(s’,a’)$的动作$a’$：
$$Q^{*}(s,a) = E_{s\sim E}[r+ \gamma \max_{a’} Q^{*}(s’,a’)|s,a], \tag{1}$$
强化学习的一个思路就是使用Bellman optimal equation更新动作值函数，$Q_{i+1}(s,a) = E[r + \gamma Q_i(s’,a’)|s,a]$，当$i\rightarrow \infty$时，$Q_i \rightarrow Q^{*}$。
上述例子是state-action pair很少的情况，当有无穷多个的时候，是无法精确计算的。这时候可以采用函数来估计动作值函数，$Q(s,a;\theta) \approx Q^{*}(s,a)$。一般来说，通常采用线性函数进行估计，当然可以采用非线性的函数，如神经网络等等。这里采用的是神经网络，用$\theta$表示网络的参数，这个网络叫做Q网络，Q网络通过最小化下列loss进行训练：
$$L_i(\theta_i) = E_{s,a\sim \rho(\cdot)}\left[(y_i - Q(s,a;\theta_i))^2\right]\tag{2}$$
其中$y_i = E_{s’\sim E}[r+\gamma \max_{a’}Q(s’,a’;\theta_{i-1})]$是第$i$次迭代的target值，其中$\rho(s,a)$是$(s,a)$服从的概率分布。
注意在优化$L_i(\theta_i)$时，上一次迭代的$\theta_{i-1}$是不变的，target取决于网络参数，和监督学习作对比，监督学习的target和网络参数无关。
对Loss函数进行求导，得到下列的gradient信息：
$$\nabla_{\theta_i}L_i(\theta_i) = E_{s,a\sim \rho(\cdot),s’\sim E}\left[(r+\gamma \max_{a’}Q(s’,a’;\theta_{i-1})-Q(s,a;\theta_i))\nabla_{\theta_i}Q(s,a;\theta_i)\right]\tag{3}$$
通过SGD优化loss函数。如果权重是每隔几个timestep进行更新，并且用从分布$\rho$和环境$E$中采样得到的样本取代期望，就可以得到熟悉的Q-learning算法[2]。(这个具体为什么是这样，我也不清楚，可以看参考文献2)
什么是on-polciy算法：

On-policy methods attempt to evaluate or improve the policy that is used to make decisions, whereas off-policy methods evaluate or improve a policy different from that used to generate the data.

Sarsa和Q-learning的区别在于更新Q值时的target policy和behaviour policy是否相同。我觉的是是policy evaluation和value iteration的区别，policy evaluation使用动态规划算法更新$V(s)$，但是并没有改变行为策略，更新迭代用的数据都是利用之前的行为策略生成的。而值迭代是policy evaluation+policy improvement，每一步都用贪心策略选择出最大的$a$更新$V(s)$，target policy（greedy）和behaviour policy（$\varepsilon$-greedy）是不同的。

强化学习需要解决的问题

大量有标记的训练数据。
delayed-reward。这个delay存在于action和reward之间，可以达到几千个timesteps那么远，和supervised learnign中输入和输入之间直接的关系相比要复杂的多。
大多数深度学习算法假设样本之间都是独立的，然而强化学习的一个sequence(序列)通常是高度相关的。
强化学习算法学习到的policy变化时，数据服从的分布通常会改变，然而深度学习通常假设数据服从一个固定的分布。

DQN

论文名称Playing Atari with Deep Reinforcement Learning

概述

DQN算法使用卷积神经网络代替Q-learning中tabular的值函数，并提出了几个trick促进收敛。DQN agnet的输入是原始的图片，输出是图片表示的state可能采取的action的$Q$值。

dqn是Model-Free的，它直接从环境$E$中采样，并没有显式的对环境进行建模。
dqn是一个online的方法，即训练数据不断增加；offline是训练数据固定。
dqn是一个off-policy算法，target policy 是greedy policy，behaviour policy是$\varepsilon$-greedy policy，target policy和greedy policy策略不同。
DQN是不收敛的。

解决方案

Experience replay

DQN使用了experience replay，将多个episodes中的经验存储到一个大小为$N$的replay buffer中。在更新$Q$值的时候，从replay buffer中进行采样更新。behaviour policy是$\varepsilon$-greedy策略，保持探索。target policy是$\varepsilon$ greedy 算法，因为replay buffer中存放的都是behaviour policy生成的experience，所以是off-policy算法。
采用experience replay的DQN和Q-learning算法相比有三个好处，第一个是每一个experience可以多次用来更新参数，提高了数据训练效率；第二个是直接从连续的样本中进行学习是低效的，因为样本之间存在强关联性。第三个是on-policy的学习中，当前的参数决定下一次采样的样本，就可能使学习出来的结果发生偏移。
replay buffer中只存储最近N个experience。

Data preprocess

原始图像是$210\times 160$的RGB图像，预处理首先将它变为灰度图，并进行下采样得到一个$110\times 84$的图像，然后从这个图像中截取一个$84\times 84$的图像。
作者使用预处理函数$\phi$处理连续四张的图像而不是一张，然后将这个预处理后的结果输入$Q$函数。
预处理函数$\phi$是一个卷积神经网络，输入是$84\times 84\times 4$的图像矩阵，经过$16$个stride为$4$的$8\times 8$filter，经过relu激活函数，再经过$32$个stride为$2$的$4\times 4$filter，经过relu激活函数，最后接一个256个单元的全连接层。输出层的大小根据不同游戏的动作个数决定。
$Q$网络的输入是预处理后的图像state，输出是所有当前state可能采取的action的$Q$值。

网络结构

输入：[batch_size, 84, 84, 4]
第一个隐藏层：16个步长为$4$的$8\times 8$的filters
第二个隐藏层：32个步长为$2$的$4\times 4$的filters
全连接层：256个units
输出层：softmax

算法

算法 1 Deep Q-learning with Experience Replay
Initialize replay memory D to capacity N
Initialize action-value function Q with random weights
for episode = $1, M$ do
$\ \ \ \ \ \ \ \ $Initialize sequence $s_1 = {x_1}$ and preprocessed sequenced $\phi_1 = \phi(s_1)$
$\ \ \ \ \ \ \ \ $for $t = 1,T$ do
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $With probability $\varepsilon$ select a random action $a_t$
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $otherwise select $a_t = \max_a Q^{∗}(\phi(s_t), a; θ)$
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $Execute action $a_t$ in emulator and observe reward $r_t$ and image $x_{t+1}$
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $Set $s_{t+1} = s_t, a_t, x_{t+1}$ and preprocess $\phi_{t+1} = \phi(s_{t+1})$
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $Store transition $(\phi_t, a_t, r_t, \phi_{t+1})$ in D
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $Sample random minibatch of transitions $(\phi_j, a_j, r_j, \phi_{j+1})$ from D
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $Set $y_j = \begin{cases}r_j&\ \ \ \ for\ terminal\ \phi_{j+1}\\r_j+\gamma \max_{a’}Q(\phi_{j+1},a’|\theta)&\ \ \ \ for\ non-terminal\ \phi_{j+1}\end{cases}$
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $Perform a gradient descent step on $(y_j − Q(\phi_j, a_j|θ))^2$
$\ \ \ \ \ \ \ \ $end for
end for

Experiments

Datasets

七个Atari 2600 games: B.Rider, Breakout, Enduro, Pong, Q bert, Seaquest, S.Invaders。
在六个游戏上DQN是SOTA，在三个游戏上DQN的表现超过了人类。

Settings

不同游戏的reward变化很大，这里把正的reward全部设置为$1$，把负的reward全部设置为$-1$，reward为$0$的保持不变。这样子在不同游戏中也可以统一学习率。
采用RMSProp优化算法，batch size为$32$，behaviour policy采用的是$\varepsilon$-greedy，在前$100$万步内，$\varepsilon$从$1$变到$0.1$，接下来保持不变。
使用了fram-skip技术，每隔$k$步，agent才选择一个action，在中间的$k-1$步中，保持原来的action不变。这里选择了$k=4$，有的游戏设置的为$k=3$。
超参数设置没有说

Metrics

每个agent训练$10$ millions帧，replay buffer size是$1$ million。每个epoch进行$50000$个minibatch weight updates或者大约$30$分钟的训练（这里有些不理解）。然后使用$\epsilon$-greedy($\epsilon=0.05$) evaluation $10000$个steps。

average total reward

第一个metric是在一个episode或者一次游戏内total reward的平均值。这个metric带有很大噪音，因为policy权值一个很小的改变可能就会对policy访问states的分布造成很大的影响。

action value function

第二个metric是估计的action-value function，这里作者的做法是在训练开始前使用random policy收集一个固定的states set，然后track这个set中states最大预测$Q$值的平均。尽管缺乏理论收敛保证，DQN看起来还不错。

Baselines

Sarsa
Contingency
DQN
Human

代码

https://github.com/devsisters/DQN-tensorflow

Nature DQN

非线性拟合函数不收敛的原因

序列中状态的高度相关性。
$Q$值的一点更新就会对policy改变造成很大的影响，从而改变数据的分布。
待优化的$Q$值和target value(目标Q值)之间的关系，每次优化时的目标Q值都是固定上次的参数得来的，优化目标随着优化过程一直在变。
前两个问题是通过DQN中提出的replay buffer解决的，第三个问题是Natura DQN中解决的，在一定时间步内，固定target network参数，更新待network的参数，然后每隔固定步数将network的参数拷贝给target network。

This instability has several causes: the correlations present in the sequence of observations, the fact that small updates to Q may significantly change the policy and therefore change the data distribution, and the correlations between the action-values (Q) and the target values $r+\gamma \max_{a’}Q(s’,a’)$.
We address these instabilities with a novel variant of Q-learning, which uses two key ideas. First, we used a biologically inspired mechanism termed experience replay that randomizes over the data, thereby removing correlations in the observation sequence and smoothing over changes in the data distribution. Second, we used an iterative update that adjusts the action-values (Q) towards target values that are only periodically updated, thereby reducing correlations with the target.

解决方案

预处理的结构变了,CNN的层数增加了一层，
加了target network，
将error限制在$[-1,1]$之间。

clip the error term from the update $r + \gamma \max_{a’} Q(s’,a’;\theta_i^{-} - Q(s,a;\theta_i)$ to be between $-1$ and $1$. Because the absolute value loss function $|x|$ has a derivative of $-1$ for all negative values of $x$ and a derivative of $1$ for all positive values of $x$, clipping the squared error to be between $-1$ and $1$ corresponds to using an absolute value loss function for errors outside of the $(-1,1)$ interval.

框架和网络结构

框架

Nature-DNQ的框架如下所示
ndqn

网络结构

输入:[batch_size, 84, 84, 4]
三个卷积层，两个全连接层（包含输出层）
第一个隐藏层：$32$个步长为$4$的$8\times 8$filters，以及一个relu
第二个隐藏层：$64$个步长为$2$的$4\times 4$filters，以及一个relu
第三个隐藏层：$64$个步长为$1$的$3\times 3$filters，以及一个relu
第四个隐藏层：$512$个units
输出层：softmax，输出每个action对应的$Q$值

算法

算法 2 deep Q-learning with experience replay, target network
Initialize replay memory D to capacity N
Initialize action-value function Q with random weights $\theta$
Initialize target action-value function $\hat{Q}$ with weights $\theta^{-}=\theta$
for episode = $1, M$ do
$\ \ \ \ \ \ \ \ $Initialize sequence $s_1 = {x_1}$ and preprocessed sequenced $\phi_1 = \phi(s_1)$
$\ \ \ \ \ \ \ \ $for $t = 1,T$ do
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $With probability $\varepsilon$ select a random action $a_t$
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $otherwise select $a_t = \max_a Q^{∗}(\phi(s_t), a; θ)$
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $Execute action $a_t$ in emulator and observe reward $r_t$ and image $x_{t+1}$
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $Set $s_{t+1} = s_t, a_t, x_{t+1}$ and preprocess $\phi_{t+1} = \phi(s_{t+1})$
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $Store transition $(\phi_t, a_t, r_t, \phi_{t+1})$ in D
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $Sample random minibatch of transitions $(\phi_j, a_j, r_j, \phi_{j+1})$ from D
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $Set $y_j = \begin{cases}r_j&\ \ \ \ for\ terminal\ \phi_{j+1}\\r_j+\gamma \max_{a’}Q(\phi_{j+1},a’|\theta^{-})&\ \ \ \ for\ non-terminal\ \phi_{j+1}\end{cases}$
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $Perform a gradient descent step on $(y_j − Q(\phi_j, a_j|θ))^2$ with respect to the network parameters $\theta$
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $Every $C$ steps reset $\hat{Q} = Q$
$\ \ \ \ \ \ \ \ $end for
end for

Experiments

Settings

batch-size: 32
replacy memory size: 1000000 frames
target network update frequency: 10000
discount factor: 0.99
action repeat: 4 # 就是frame skip
history length: 4 # 使用最近的几帧重叠，实际上是16帧
paramteter update frequency: 4 # 执行sgd train的frequency
learning rate: 0.00025
gradient momentum: 0.95
squared gradient momentum: 0.95
min squared gradient: 0.01
initial exploration: 1
final exploration 0.1
final exploration frame: 1000000
replay start size: 50000
no-op max: 30
reward: clipped to [-1, 1]
total train frames: 50 millons frame（实际上是200 millions emulated frames，因为有设置为$4$ frame skip）。
选择random作为一个baseline，因为人类的极限是$10$hz，为了公平起见，random baseline以$10$hz的频率随机选择一个action，atari视频的频率是$60$hz，所以每隔$6$帧，随机选择一个action，在选择action中间的帧中保持这一个action。

Experiments

Average score和average action value
每个training epoch之后进行一次evaluation，记录evaluation过程中average episode reward。总共train $50$million frames，大概有$200$个epoch，也就是一个epoch是$25$万frames，在每个epoch后使用$\epsilon$-greedy($\epsilon =0.05$)策略evaluate $520k$个frames。
Main-Evaluation: Compartion between DQN and other baselines
Baselines有Random play, best linear learner，SARSA，Huamn等。
DQN总共训练了$50$ million frames，replay buffer存放最近的$1$ million framems。在完成训练后，至多执行$30$次no-op，产生随机初始状态，使用$\epsilon$-greedy($\epsilon=0.05$)玩$5$分钟，对多次结果取平均。
Human的数据是在玩家首先进行了$2$个小时训练后，然后玩大约20 episodes，每个episode最长5 min的 average reward。
表格中最后一列还给出了一个百分比，$100\times \frac{\text{DQN score} - \text{random play score}}{\text{human score} - \text{random play score}}$。
Replay buffer和target network的abalation实验
在三个不同的learning rate下使用standard hyperparameters训练DQN $10$ million frames。每隔$250,000$ training frames对每个agent进行$135,000$ frames的validation，记录最高的average episode score。这些valuation episodes并没有在$5$ min的时候截断，这个实验中Enduro得到的score要比main evaluation中高。这个实验中training的帧数($10$ million frames)要比baseline中training的frames($50$ million frames)少。
DQN和linear function approximator比较
除了function approximator由CNN变成linear的，其他都没有变。
在$5$个validation games上，每个agent在三个不同的learning rates使用标准的参数训练了$10$ million frames。每隔$250000$个training frames对agent进行$135,000$ frames的validation，reported 最高的average episode score。这些valuation episodes并没有在$5$ min的时候截断，这个实验中Enduro得到的score要比main evaluation中高。这个实验中training的帧数($10$ million frames)要比baseline中training的frames($50$ million frames)少。

Gorila DQN

论文名称：
Massively Parallel Methods for Deep Reinforcement Learning
下载地址：
https://arxiv.org/pdf/1507.04296.pdf

Experiments

Settings

网络结构和nature DQN一样。
使用了frame-skip，设置为$4$。
Replay memory是$1$M>

Evaluation

Evaluation有两种：

null op starts

每个agent在它训练的游戏上evaluated $30$个episodes，每个episode随机的至多执行$30$次no-op之后，评估$5$min的emulator时间($18000$ frames)。然后取这$100$次的平均值。

human starts

Human starts是用来衡量它对于agent可能没有遇到过的state的泛化能力。对于每一个游戏，从一个人类玩家的gameplay中随机取$100$个开始点，使用$\epsilon$-greedy policy玩三十分钟emulator时间（即$108000$frames)。

Double DQN

DQN中的overestimate问题

解决overestimate问题，Q-learning中在estimated values上进行了max操作，可能会导致某些更偏爱overestimated value而不是underestimated values。
本文将Double Q-learning的想法推广到了dqn上形成了double-dqn。实验结果表明了overestimated value对于policy有影响，double 会产生好的action value，同时在一些游戏上会得到更高的scores。

Contributions

解释了在large scale 问题上，Q-learning被overoptimistic的原因是学习固有的estimation errors。
overestimation在实践中是很常见，也很严重的。
Double Q-learning可以减少overoptimism
提出了double-dqn。
double-dqn在某些游戏上可以找到更好的policy。

Double Q-learning

Q-learning算法计算target value $y$的公式如下：
$$y = r + \gamma \max_a’ Q(s’, a’|\theta_t)\tag{4}$$
在计算target value的时候，使用同一个网络选择和评估action $a’$，这可能会让网络选择一个overestimated value，最后得到一个overoptimistic value estimate。所有就有了double Q-learning，计算公式如下：
$$y = r + \gamma Q(s’, \arg\max_a’ Q(s’,a;\theta_t);\theta’_t)\tag{5}$$
target policy还是greedy policy，通过使用$\theta$对应的网络选择action，然后在计算target value的时候使用$\theta’$对应的网络。
原有的公式可以写成下式，
$$y = r + \gamma Q(s’, \arg\max_a’ Q(s’,a;\theta_t);\theta_t)\tag{6}$$
即选择action和计算target value都是使用的同一个网络。

Double DQN

double-dqn
Double Q-learnign的做法是分解target action中的max opearation为选择和evaluation。而在Nature-dqn中，提出了target network，所以分别使用network和target network去选择和evaluation action是一个很好的做法，这样子公式就变成了
$$y = r + \gamma Q(s’, \arg\max_a’ Q(s’,a;\theta_t);\theta^{-}_t)\tag{7}$$
和Q-learnign相比，将$\theta’$换成了$\theta^{-}$ evaluate action，target network的更新和nature-dqn一样，过一段时间复制network的参数。

Double Q learning vs Q-learning

可以在数学上证明，Q-learning是overestimation的，但是double q leraing是无偏的。。。证明留待以后再说。

网络结构

网络结构和nature DQN一样。

算法

算法 3: Double DQN Algorithm.
输入: replay buffer $D$, 初始network参数$\theta$,target network参数$\theta^{-}$
输入 : replay buffer的大小$N_r$, batch size $N_b$, target network更新频率$N^{-}$
for episode $e \in {1, 2,\cdots, M}$ do
$\qquad$初始化frame sequence $\mathbf{x} \leftarrow ()$
$\qquad$for $t \in {0, 1, \cdots}$ do
$\qquad\qquad$设置state $s \leftarrow \mathbf{x}$, 采样 action $a \sim\pi_B$
$\qquad\qquad$给定$(s, a)$，从环境$E$中采样接下来的frame $x_t$,接收reward $r$,在序列$\mathbf{x}$上拼接$x$
$\qquad\qquad$if $|\mathbf{x}| \gt N_f$
$\qquad\qquad$then
$\qquad\qquad\qquad$从$\mathbf{x}$中删除最老的frame $x_{t_min}$
$\qquad\qquad$设置$s’ \leftarrow \mathbf{x}$,添加transition tuple (s, a, r, s 0 ) 到buffer D中，如果$|D| \ge N_r$替换最老的tuple
$\qquad\qquad$采样$N_b$个tuples $(s, a, r, s’) \sim Unif(D)$
$\qquad\qquad$计算target values, one for each of $N_b$ tuples:
$\qquad\qquad$定义$a^{\max}(s’; \theta) = \arg \max_{a’} Q(s’, a’;\theta)$
$\qquad\qquad y_j = \begin{cases}r&\qquad if\ \ s’\ \ is\ \ terminal\\ r+\gamma Q(s’, a^{\max}(s’;\theta);\theta^{-}, &\qquad otherwise\end{cases}$
$\qquad\qquad$利用loss $||y_j − Q(s, a; \theta)||^2$的梯度更新
$\qquad\qquad$每隔$N^{-}$个步骤更新一下target network 参数$\theta^{-}$
$\qquad$end
end

Experiments

Settings

Tunned Double DQN，update frequency从$10000$改成了$30000$，训练时$\epsilon$在$1$ millon内从$0.1$退火到$0.01$。Evaluation时是$0.001$。

Evaluation

和Gorila DQN一样，用了两种：no-op和human starts。

Training

在每个游戏上，网络都训练了$200$M frames，也就是$50$M steps。每隔$1$M step进行一次evaluation，从evaluations中选出最好的policy作为输出。

Metric

提出了一个指标，normalized score，计算公式如下：
$$score_{normalized} = \frac{score_{agent}- score_{random}}{score_{human}-score_{random}}\tag{8}$$
分母是human和random之差，对应$100%$。

Prioritized DQN(PER)

contributions

本文提出一种了proritizing experience的框架，在训练过程中多次使用重要的transtions replay进行更新，让训练变得的更有效率。
使用TD-errors作为prioritization mechanism，给出了两种protitization计算方式，提出了一种stochastic prioritization以及importance sampling方法。

Prioritized replay

可以从两个维度上考虑replay memeory的改进，一个是存哪些experiences，一个是使用哪些experiences进行回放。本文是从第二个方向上进行的考虑。

从buffer中随机抽样的方法中，update steps和memory size是线性关系，作者想找一个update steps和memory size是log关系的oracle，但是很遗憾，这是不现实的，所以作者想要找一种比uniform random replay好尽量接近oracle的方法。

Prioritizion with TD-error

prioritized replay最重要的部分是如何评价每一个transition的重要程度。一个理想的criterion是agent在当前的state可以从某个transition中学到多少。这个measure metric是不确定的，一个替代方案是使用TD error $\delta$，表示how ‘suprising’ 或者upexpected the transition：就是当前的value离next-step bootstrap得到的value相差多少，booststrap就是基于其他估计值进行计算。。这中方法对于incremental,online RL方法，例如SARSA以及Q-learning来说都是很合适的，因为他们会计算TD-error，然后给TD-error一个比例系数用来更新参数。然后当reward是noisy的时候，TD-error效果可能很差。
作者在一个人工设计的环境中使用了greedy TD-error prioritization算法，算法在每次存transition到replay buffer的时候，同时还会存一下该transition最新的TD-error，然后在更新的时候从memory中选择TD-error最大的transition。最新的transition TD-error没有算出来，就给它一个最大的priority，保证所有的experience都至少被看到过一次。
采用二叉堆用实现优先队列，查找复杂度是$O(1)$，更新priorities的复杂度是$O(logN)$。

Stochastic prioritization

上述方法有很多问题。第一，每次都sweep整个replay memory的计算量很大，所以只有被replayed的experiences的TD-errors才会被更新。开始时一个TD error很小的transition可能很长一段事件不会被replayed，这就导致了replay buffer的sliding window不起作用了。第二，TD-error对于noise spike很敏感，还会被bootstrap加剧，估计误差可能会是另一个noise。第三，greedy prioritization集中在experiences的一个subset：errors减小的很慢，尤其是使用function appriximation时，这就意味着初始的高error的transitions会被replayed的很频繁，然后会over-fitting因为缺乏diversity。
为了解决这些问题，引入了一个介于pure greedy prioritizaiton以及uniform random sampling之间的stochastic采样方法，priority高的transition有更大的概率被采样，而lowest-priority的transition也有概率被选中，具体的定义transition $i$的概率如下：
$$P(i) = \frac{p_i^{\alpha}}{\sum_kp_k^{\alpha}}\tag{9}$$
$\alpha$确定prioritizaiton的比重，如果$\alpha=0$就是unifrom。

有两种$p_i$的计算方法，一种是直接的proportional prioritization，$p_i = |\delta_i| + \varepsilon$，其中$\varepsilon$是一个小的正整数，确定当$p_i=0$时，该transition仍能被replay；第二种是间接的，$p_i = \frac{1}{rank(i)}$，其中$rank(i)$是所有replay memory中的experiences根据$|\delta_i|$排序后的rank。第二种方法的鲁棒性更好。
在实现上，两种方法都有相应的trick，让复杂度不依赖于memory 大小$N$。Proportional prioritization采用了’sum-tree’数据结构，每一个节点都是它的子节点的children，priorities是leaf nodes。而rank-based方法，使用线性函数估计累计密度函数，具体怎么实现没有细看。

annealing the bias

因为random sample方法，samples之间没有一点联系，选择每一个sample的概率都是相等的，但是如果加上了priority，就有一个bias toward高priority的samples。IS和prioritized replay的组合在non-learn FA中有一个用处，large steps可能会产生不好的影响，因为梯度信息可能是局部reliable，所以需要使用一个小点的step-size。
在本文中，high-error的样本可能会观测到很多次，使用IS减小gradient的大小，对应于高priority的samples的weight被微调了一下，而对应于低priority的样本基本不变。
weigth的计算公式如下：
$$w_i = (\frac{1}{N}\cdot \frac{1}{P(i)})^{\beta}\tag{10}$$
OK,这里IS的作用有些不明白。。。。

算法

算法 4
输入: minibatch $k$, 学习率（步长）$\eta$, replay period $K$ and size $N$ , exponents $\alpha$ and $\beta$, budget $T$.
初始化replay memory $H = \emptyset, \Delta = 0, p_1 = 1$
根据$S_0$选择 $A_0 \sim \pi_{\theta}|(S_0)$
for $t = 1,\cdots, T$ do
$\qquad$观测$S_t, R_t, \gamma_t$
$\qquad$存储transition $(S_{t−1}, A_{t−1}, R_t , \gamma_t, S_t)$ 到replay memory，以及$p_t$的最大priority $p_t = \max {i\lt t} p_i$
$\qquad$if $t ≡ 0$ mod $K$ then
$\qquad\qquad$for j = 1 to k do
$\qquad\qquad\qquad$Sample transition $j \sim P(j) = \frac{p_j^{\alpha}}{\sum_i p_i^{\alpha}}$
$\qquad\qquad\qquad$计算importance-sampling weight $w_j = \frac{(N \cdot P(j))^{\beta}}{\max_i w_i}$
$\qquad\qquad\qquad$计算TD-error $\delta_j = R_j + \gamma_j Q_{target} (S_j$, $arg \max_a Q(S_j, a)) − Q(S_{j−1} , A\ {j−1})$
$\qquad\qquad\qquad$更新transition的priority $p_j \leftarrow |\delta_j|$
$\qquad\qquad\qquad$累计weight-change $\Delta \leftarrow \Delta + w_j \cdot \delta_j \cdot \nabla_{\theta} Q(S_{j−1}, A_{j−1})$
$\qquad\qquad$end for
$\qquad\qquad$更新weights $\theta\leftarrow \theta+ \eta\cdot\Delta$, 重置$\Delta = 0$
$\qquad\qquad$每隔一段时间更新target network $\theta_{target} \leftarrow \theta$
$\qquad$end if
$\qquad$选择action $A_t \sim \pi_{\theta}(S_t)$
end for

Experiments

两组实验，
一组是DQN和proportional prioritization作比较。
一组是tuned Double DQN和rank-based以及proportional prioritizaiton。

Metrics

用的是double dqn提出来的nomalized score，这里在分母上加了绝对值。
主要用的median scores和mean scores。

Dueling DQN

介绍

本文作者提出来将dueling网络框架应用在model-free算法上。The dueling architecture能用一个deep model同时表示$V(s)$和优势函数$A(s,a)$，网络的输出将$V$和$A$结合产生$Q(s,a)$。和advantage不一样的是，这种方式在构建时就将他们进行了解耦，因此，dueling architecture可以应用在各种各样的model free RL算法上。
本文的架构是对算法创新的补充，它可以对之前已有的各种DQN算法进行结合。

dueling network architecture

这个新的architecture的核心想法是，没有必要估计所有states的action value。在一些states，需要action value去确定执行哪个action，但是在许多其他states，action values并没有什么用。当然，对于bootstrap算法来说，每一个state的value estimation都很重要。
dueling-dqn
作者给出了一个single Q-network的architecture，如图所示。
网络结构和nature-dqn一样，但是这里加了两个fully connected layers，一个用于输出$V$，一个用于输出$A$。然后$A$和$V$结合在一起，产生$Q$，网络的输出和nature dqn一样，对应于某个state的一系列action value。
从$Q$函数的定义$Q^{\pi}(s,a) = V^{\pi}(s)+A^{\pi}(s,a)$以及$Q$和$V$之间的关系$V^{\pi}(s) = \mathbb{E}_{a\sim\pi(s)}\left[Q^{\pi}(s,a)\right] = \pi(a|s)Q^{\pi}(s,a)$，所以有$\mathbb{E}_{a\sim\pi(s)}\left[A^{\pi}(s,a)\right]=0$。此外，对于deterministic policy，$a^{*} = \arg \max_{a’\in A}Q(s,a’)$，有$V(s) = Q(s,a^{*})$，即$A(s,a^{*}) = 0$。
如图所示的network中，一个网络输出scalar $V(s;\theta, \beta)$，一个网络输出一个$|A|$维的vector $A(s,a;\theta, \alpha)$，其中$\theta$是网络参数，$\alpha$和$\beta$是两个全连接层的参数。
根据advantage的定义，可以直接将他们加起来，即：
$$Q(s,a;\theta, \alpha, \beta) = V(s;\theta, \beta) + A(s,a;\theta, \alpha) \tag{11}$$
但是，我们需要知道的一点是，$Q(s, a;\theta, \alpha, \beta)$仅仅是$Q$的一个参数化估计。它由两部分组成，一部分是$V$，一部分是$A$，但是需要注意的是，这里的$V$和$Q$只是我们叫它$V$和$A$，它的实际意义并不是$V$和$A$。给了$Q$，我们可以得到任意的$Q(s, a) = V(s) + A(s,a)$，而$V$和$Q$并不代表value function和advantage functino。
为了解决这个问题，作者提出了选择让advantage为$0$的action，即：
$$Q(s, a; \theta,\alpha, \beta) = V(s; \theta, \beta) + \left(A(s,a;\theta,\alpha) - \max_{a’\in |A|}A(s, a’; \theta, \alpha)\right)\tag{12}$$
选择$a^{*} = \arg \max_{a’\in A} Q(s, a’; \theta, \alpha, \beta) = \arg \max_{a’\in A}A(s, a’;\theta, \alpha)$，我们得到$Q(s,a^{*}; \theta, \alpha,\beta) = V(s;\theta, \beta)$。这个时候，输出$V$的网络给出的真的是state value的估计$V(s;\theta, \beta)$，另一个网络真的给出的是advantage的估计。
另一种方法是用mean取代max操作：
$$Q(s, a; \theta,\alpha, \beta) = V(s; \theta, \beta) + \left(A(s,a;\theta,\alpha)- \frac{1}{|A|}\sum_{a’}A(s, a’; \theta, \alpha)\right)\tag{13}$$
一方面这种方法失去了$V$和$A$的原始语义，因为它们有一个常数的off-target，但是另一方面它增加了优化的稳定性，因为上式中advantage的改变只需要和mean保持一致即可，不需要optimal action’s advantange一有变化就要改变。

算法

Distributed DQN

Noisy DQN

介绍

已有方法的exploration都是通过agent policy的random perturbations，比如常见的$\varepsilon$-greedy等方法。这些方法不能找出环境中efficient exploration的behavioural patterns。常见的方法有以下几种:
第一种方法是optimism in the face of uncertainty，理论上证明可行，但是通常应用在state-action spaces很小的情况下或者linear FA，很难处理non-linearn FA，而且non-linear情况下收敛性没有保证。
另一种方法是添加额外的intrinsic motivation term，该方法的问题是将算法的generalisation mechanism和exploration分割开，即有instrinsic reward和environment reward，它们的比例如何去设置，需要认为指定。如果不仔细调整，optimal policy可能会受intrinsic reward影响很大。此外为了增加exploration的鲁邦性，扰动项仍然是需要的。这些算法很具体也能应用在参数化policy上，但是很低效，而且需要很多次policy evaluation。
本文提出NoisyNet学习网络参数的perturbations，主要想法是参数的一点改变可能会导致policy在很多个timsteps上的consistent，complex, state-dependent的变化，而如$\varepsilon$-greedy的dithering算法中，每一步添加到policy上的noise都是不相关的。pertubations从一个noise分布中进行采样，它的variance可以看成noise的energy，variance的参数和网络参数都是通过loss的梯度进行更新。网络参数中仅仅加入了噪音，没有distribution，可以自动学习。
在高维度上，本文的算法是一个randomised value function，这个函数是neural network，网络的参数并没有加倍，linear 的参数加倍，而参数是noise的一个简单变换。
还有人添加constant Gaussian niose到网络参数，而文本的算法添加的noise并不是限制在Gaussion noise distributions。添加noise辅助训练在监督学习等任务中一直都有，但是这些噪音都是不能训练的，而NoisyNet中的噪音是可以梯度下降更新的。

NoisyNets

noisy_linear_layer
用$\theta$表示noisy net的参数，输入是$x$，输出是$y$，即$y=f_{\theta}(x)$。$\theta$定义为$\theta=\mu+\Sigma\odot\varepsilon$，其中$\zeta=(\mu,\Sigma)$表示可以学习的参数，$\varepsilon$表示服从固定分布的均值为$0$的噪音,$\varepsilon$是random variable。$\odot$表示element-wise乘法。最后的loss函数是关于$\varepsilon$的期望：$\bar{L}(\zeta)=\mathbb{E}\left[L(\theta)\right]$，然后优化相应的$\zeta$，$\varepsilon$不能被优化，因为它是random variable。
一个有$p$个输入单元，$q$个输出单元的fully-connected layer表示如下：
$$y=wx+b \tag{14}$$
其中$w\in \mathbb{R}^{q\times p}$，$x\in \mathbb{R}^{p}$,$b\in \mathbb{R}^{q}$，对应的noisy linear layer定义如下：
$$y=(\mu^w+\sigma^w\odot\varepsilon^w)x + \mu^b+\sigma^b\odot\varepsilon^b \tag{15}$$
就是用$\mu^w+\sigma^w\odot\varepsilon^w$取代$w$，用$\mu^b+\sigma^b\odot\varepsilon^b$取代$b$。其中$\mu^w,\sigma^w\in \mathbb{R}^{q\times p} $，而$\mu^b,\sigma^b\in\mathbb{R}^{q}$是可以学习的参数，而$\varepsilon^w\in \mathbb{R}^{p\times q},\varepsilon^b \in \mathbb{R}^{q}$是random variable。
作者提出了两种添加noise的方式，一种是Independent Gaussian noise，一种是Factorised Gaussion noise。使用Factorised的原因是减少随机变量的计算时间，这些时间对于单线程的任务来说还是很多的。

Independent Gaussian noise

应用到每一个weight和bias的noise都是independent的，对于$\varepsilon^w$的每一项$\varepsilon_{i,j}^w$来说，它们的值都是从一个unit Gaussion distribution中采样得到的；$varepsilon^b$同理。所以对于一个$p$个输入,$q$个输出的noisy linear layer总共有$pq+q$个noise 变量。

Factorised Gaussian noise

通过对$\varepsilon_{i,j}^w$来说，可以将其分解成$p$个$\varepsilon_i$用于$p$个输入和$q$个$\varepsilon_j$用于$q$个输出，总共有$p+q$个noiss变量。每一个$\varepsilon_{i,j}^w$和$\varepsilon_{j}^b$可以写成：
$$\varepsilon_{i,j}^w = f(\varepsilon_i)f(\varepsilon_j) \tag{16}$$
$$\varepsilon_{j}^b = f(\varepsilon_j)\tag{17}$$
其中$f$是一个实函数，在第一个式子中$f(x) = sng(x)\sqrt{|x|}$，在第二个式子中可以取$f(x)=x$，这里选择了和第一个式子中一致。
因为noisy network的loss函数是$\bar{L}(\zeta)=\mathbb{E}\left[L(\theta)\right]$，是关于noise的一个期望，梯度如下：
$$\nabla\bar{L}(\zeta)=\nabla\mathbb{E}\left[L(\theta)\right]=\mathbb{E}\left[\nabla_{\mu,\Sigma}L(\mu+\Sigma\odot\varepsilon)\right] \tag{18}$$
使用Monte Carlo估计上述梯度，在每一个step采样一个sample进行optimization:
$$\nabla\bar{L}(\zeta)\approx\nabla_{\mu,\Sigma}L(\mu+\Sigma\odot\varepsilon) \tag{19}$$

Noisy DQN and dueling

相对于DQN和dueling DQN来说，noisy DQN and dueling主要做了两方面的改进：

不再使用$\varepsilon$-greedy behaviour policy了，而是使用greedy behaviour policy采样优化randomised action-value function。
网络中的fully connected layers全都换成了参数化的noisy network，noisy network的参数在每一次replay之后从noise服从的distribution中进行采样。这里使用的nose是factorised Gaussian noise。

在replay 整个batch的过程中，noisy network parameter sample保持不变。因为DQN和Dueling每执行一个action step之后都会执行一次optimization，每次采样action之前都要重新采样noisy network parameters。

Loss

$Q(s,a,\epsilon;\zeta)$可以看成$\zeta$的一个random variable，NoisyNet-DQN loss如下：
$$\bar{L}(\zeta) = \mathbb{E}\left[\mathbb{E}_{(x,a,r,y)}\sim D\left[r + \gamma \max_{b\in A}Q(y, b, \varepsilon’;\zeta^{-}) - Q(x,a,\varepsilon;\zeta)\right]^2\right]\tag{20}$$
其中外层的期望是$\varepsilon$相对于noisy value function $Q(x,a, \varepsilon;\zeta)$和$\varepsilon’$相对于noisy target value function $Q(x,a, \varepsilon’;\zeta^{-}$。对于buffer中的每一个transition，计算loss的无偏估计，只需要计算target value和true value即可，为了让target value和true之间没有关联，target network和online network采用independent noises。
就double dqn中的action选择来说，采样一个新的independent sample $\varepsilon^{’’}$计算action value，然后使用greedy操作，NoisyNet-Dueling的loss如下：
$$\bar{L}(\zeta) = \mathbb{E}\left[\mathbb{E}_{(x,a,r,y)}\sim D\left[r + \gamma Q(y, b^{*}(y), \varepsilon’;\zeta^{-} - Q(x,a,\varepsilon;\zeta)\right]^2\right]\tag{21}$$
$$b^{*}(y) = \arg \max_{b\in A} Q(y, b(y), \varepsilon^{’’};\zeta)\tag{22}$$

Noisy-A3C

Noisy-A3C相对于A3C有以下的改进：

entropy项被去掉了;
fully-connected layer被替换成了noisy network。

A3C算法中没有像$\epsilon$-greedy这样进行action exploration，选中的action通常是从current policy中选的，加入entropy是为了鼓励exploration，而不是选择一个deterministic policy。当添加了noisy weights时，对参数进行采样就表示选择不同的current policy，就已经代表了exploration。NoisyNet相当于直接在policy space中进行exploration，而entropy项就可以去掉了。

Noisy Networks的初始化

在unfactorised noisy networks中，每个$\mu_{i,j}$从独立的均匀分布$U\left[-\sqrt{\frac{3}{p}}, \sqrt{\frac{3}{p}}\right]$中采样初始化，其中$p$是对应linear layer的输入个数，$\sigma_{i,j}$设置为一个常数$0.0017$，这是从监督学习的任务中借鉴的。
在factorised noisy netowrks中，每个$\mu_{i,j}$从独立的均匀分布$U\left[-\sqrt{\frac{1}{p}}, \sqrt{\frac{1}{p}}\right]$中进行采样，$\sigma_{i,j}$设置为$\frac{\sigma_0}{p}$，超参数$\sigma_0$设置为$0.5$。

算法

算法5 NoisyNet-DQN / NoisyNet-Dueling
输入: Env Environment; $\varepsilon$ random variables of the network的集合
输入: DUELING Boolean; "true"代表NoisyNet-Dueling and "false"代表 NoisyNet-DQN
输入: $B$空replay buffer; $\zeta$初始的network parameters; $\zeta^{-}$初始的target network parameters
输入: replay buffer大小$N_B$; batch size $N_T$; target network更新频率$N^{-}$
输出: $Q(\cdot, \varepsilon; \zeta)$ action-value function
for episode $e\in {1,\cdots , M}$ do
$\qquad$初始化state sequence $x_0 \sim Env$
$\qquad$for $t \in {1,\cdots }$ do
$\qquad\qquad$设置$x \leftarrow x_0$
$\qquad\qquad$采样 a noisy network $\xi\sim \varepsilon$
$\qquad\qquad$选择an action $a \leftarrow \arg \max_{b\in A} Q(x, b, \xi; \zeta)$
$\qquad\qquad$采样 next state $y \sim P (\cdot|x, a)$, 接收 reward $r \leftarrow R(x, a) $以及$x_0 \leftarrow y$
$\qquad\qquad$将transition (x, a, r, y)添加到replay buffer
$\qquad\qquad$if $|B| \gt N_B$ then
$\qquad\qquad\qquad$删掉最老的transition
$\qquad\qquad$end if
$\qquad$采样一个大小为$N_T$的batch, transitions $((x_j, a_j, r_j, y_j) \sim D)_{j=1}^{N_T}$
$\qquad\qquad$采样noisy variables用于online network $\xi \sim\varepsilon$
$\qquad\qquad$采样noisy variables用于target network $\xi’\sim\varepsilon$
$\qquad\qquad\qquad$if DUELING then
$\qquad\qquad\qquad$采样noisy variables用于选择action的network $\xi\sim\varepsilon$
$\qquad\qquad$end if
$\qquad\qquad$for $j \in {1,\cdots, N_T}$ do
$\qquad\qquad\qquad$if $y_j$ is a terminal state then
$\qquad\qquad\qquad\qquad$$\hat{Q}\leftarrow r_j$
$\qquad\qquad\qquad$end if
$\qquad\qquad\qquad$if DUELING then
$\qquad\qquad\qquad\qquad b^{*}(y_j) = \arg \max_{b\in A} Q(y_j, b, \xi^{’’}; \zeta)$
$\qquad\qquad\qquad\qquad\qquad \hat{Q}\leftarrow r_j + \gamma Q(y_j, b^{*}(y_j), \xi’;\zeta^{-})$
$\qquad\qquad\qquad$else
$\qquad\qquad\qquad\qquad$$\hat{Q}\leftarrow r_j + \gamma \max_{b\in A} Q(y_j, b, \xi’;\zeta^{-})$
$\qquad\qquad$end if
$\qquad\qquad\qquad$利用loss $(\hat{Q}-Q(x_j,a_j, \xi;\zeta))^2$的梯度更新$\zeta$
$\qquad\qquad$end for
$\qquad\qquad$每隔$N^{-}$步更新target network:$ \zeta^{−}\leftarrow \zeta$
$\qquad$end for
end for

算法6 NoisyNet-A3C for each actor-learner thread
输入: Environment Env, 全局共享参数$(\zeta_{\pi},\zeta_{V})$ , 全局共享counter $T$和maximal time $T_{max}$
输入: 每个线程的参数 $(\zeta’_{\pi},\zeta’_{V})$, random variables $\varepsilon$的集合, 每个线程的counter $t$和TD-$\gamma$的长度$t_{max}$
输出: policy $\pi(\cdot; \zeta_{\pi}, \varepsilon)$和value $V(\cdot; \zeta_{V}, \varepsilon)$
初始化线程counter $t \leftarrow 1$
repeat
$\qquad$重置acumulative gradients: $d\zeta_{\pi}\leftarrow 0$和$d\zeta_V \leftarrow 0$
$\qquad$Synchronise每个线程的parameters: $\zeta’_{\pi}\leftarrow \zeta_{\pi}$和$\zeta_V\leftarrow \zeta_V$
$\qquad$counter $\leftarrow 0$
$\qquad$从Env中得到state $x_t$
$\qquad$采样noise: $\xi\sim\varepsilon$
$\qquad r \leftarrow []$
$\qquad a \leftarrow []$
$\qquad x \leftarrow []$和$x[0] \leftarrow x_t$
$\qquad$repeat
$\qquad\qquad$采样action: $a_t \sim\pi(\cdot|x_t;\zeta’_{\pi};\xi)$
$\qquad\qquad$$a[−1]\leftarrow a_t$
$\qquad\qquad$接收reward $r_t$和next state $x_{t+1}$
$\qquad\qquad$$r[−1]\leftarrow r_t$和$x[−1]\leftarrow x_t+1$
$\qquad\qquad$$t\leftarrow t + 1$和 $T\leftarrow T + 1$
$\qquad\qquad$$counter = counter + 1$
$\qquad\qquad$until $x_t\ \ terminal\ \ or\ \ counter == t_{max} + 1$
$\qquad$if $x_t$ is a terminal state then
$\qquad\qquad$$Q = 0$
$\qquad$else
$\qquad\qquad$$Q = V(x_t; \zeta’_{V}, \xi)$
$\qquad$end if
$\qquad$for $i \in {counter − 1, \cdots, 0}$ do
$\qquad\qquad$更新Q: $Q\leftarrow r[i] + \gamma Q$
$\qquad\qquad$累积policy-gradient: $d\zeta_{\pi} \leftarrow d\zeta_{\pi} + \nabla \zeta’_{\pi}log(\pi(a[i]|x[i]; \zeta’_{\pi}, \xi))[Q − V(x[i]; \zeta’_{\pi}V, \xi)]$
$\qquad\qquad$累积 value-gradient: $d\zeta_V \leftarrow ← d\zeta_V+ \nabla \zeta’_{V}[Q − V(x[i]; \zeta’_{V}, \xi)]^2$
$\qquad$end for
$\qquad$执行$\zeta_{\pi}$的asynchronous update: $\zeta_{\pi}\leftarrow \zeta_{\pi} + \alpha_{\pi}d\zeta_{\pi}$
$\qquad$执行$\zeta_{V}$的asynchronous update: $\zeta_{V}\leftarrow \zeta_{V} − \alpha_VdV\zeta_{V}$
until $T \gt T_{max}$

Rainbow

参考文献

1.https://blog.csdn.net/yangshaokangrushi/article/details/79774031
2.https://link.springer.com/article/10.1007%2FBF00992698
3.https://www.jianshu.com/p/b92dac7a4225
4.https://datascience.stackexchange.com/questions/20535/what-is-experience-replay-and-what-are-its-benefits/20542#20542
5.https://stats.stackexchange.com/questions/897/online-vs-offline-learning
6.https://www.freecodecamp.org/news/improvements-in-deep-q-learning-dueling-double-dqn-prioritized-experience-replay-and-fixed-58b130cc5682/
7.https://jaromiru.com/2016/11/07/lets-make-a-dqn-double-learning-and-prioritized-experience-replay/
8.https://datascience.stackexchange.com/questions/32873/prioritized-replay-what-does-importance-sampling-really-do
9.https://papers.nips.cc/paper/5249-weighted-importance-sampling-for-off-policy-learning-with-linear-function-approximation.pdf

Modeling Others using Oneself in Multi-Agent Reinforcement Learning

发表于 2019-01-29 | 更新于 2019-12-17 | 分类于强化学习

摘要

我们考虑使用不完全信息的多智能体强化学习问题，每个智能体的目标是最大化自身的效用。奖励函数取决于两个智能体的隐藏状态（或者目标），每一个智能体必须从它观察到的行为中推断出其他玩家的隐藏目标从而完成任务。我们提出了一种新的方法在这些领域中进行学习：自我其他建模（SOM），智能体使用自己的策略来预测其他智能体的动作并实时更新其他智能体隐藏状态的置信度。我们在三个不同的任务上对该方法进行了评估，结果表明智能体无论在合作还是对抗环境中都能使用他们对其他玩家隐藏状态的估计来学习到更好的策略。

引言

在多智能体系统中推理其他智能体的意图并预测它们的行为是很重要的，这些智能体可能有不同的甚至是竞争的目标集。由于多智能体系统的不稳定性，这仍然是一个非常具有挑战性的问题。
在本文中，我们介绍了一种从其他智能体的行为中估计对应的未知的目标和并利用这些估计的目标选择动作的新方法。我们证明了在本文提到的任务中，在游戏中显式的对其他玩家进行建模比将其他智能体看做环境的一部分会有更好的性能。我们将问题定义为双人随机游戏，也叫双人马尔可夫游戏，其中环境对于智能体是完全可见的，但是没有关于其他智能体目标的明确知识而且没有沟通信道。每个智能体在回合结束时收到的奖励取决于两个智能体的目标，因此是每个智能体最优的策略都必须考虑到所有智能体的目标。
认知科学研究表明，人类维持与他们联系的其他人的模型，这些模型用来捕捉那些人的目标，信仰或偏好。在某些情况下，人类利用自己的心理过程来模拟他人的行为。这使他们能够理解其他人的意图或动机，并能在社交场合采取相应的行动。受这些研究的启发，关键想法是要理解游戏中其他玩家正在做什么，智能体应该问自己“如果我扮演另一个玩家的角色，我的目标是什么？”。我们通过使用一个多层循环神经网络参数化智能体的动作和值函数来实现这个想法，该神经网络将状态和目标作为输入。当智能体玩游戏时，它通过直接使用自己的动作函数优化目标来最大化对方行动的可能性，从而推断出其他智能体的未知目标。

方法

背景两个智能体的马尔可夫游戏由描述所有智能体的可能配置的一组状态集合$S$，两组动作集合$A_1$，$A_2$和两个智能体的观察$O_1$，$O_2$以及转换函数$\Tau$：$S\times A_1 \times A_2 \rightarrow S$作为当前状态和动作的函数给出下一个状态的概率分布。每个智能体$i$通过从随机策略$\pi_{\theta_i}:S\times A_i\rightarrow [0,1]$中采样选择动作。每个智能体都有一个奖励函数，它取决于智能体的状态和动作：$r_i：S\times A_i\rightarrow R$。每个智能体$i$试图最大化自己的总预期收益$R_i = \sum_{t =0}^T\gammatr_i^t$，其中$\gamma$是折扣因子，$T$是时间范围。在本文中，我们考虑了合作以及竞争环境。
接下来介绍自我其他模型（SOM），这是一种在一个回合内以实时方式推断其他智能体的目标并使用这些估计来选择动作的新方法。为了决定一个动作并估计一个状态的值，我们使用一个神经网络$f$将它自己的目标$z_{self}$，另一个玩家的估计目标$\hat{z}{self}$，并且他自己的角度的观察状态$s{self}$作为输入，输出动作$\pi$的一个概率分布和值估计$V$，即对于每个玩游戏的智能体，有：
$$\begin{bmatrix}\pi^i\Vi\end{bmatrix}=f^i(s_{self}i,z_{self}^{i,\hat{z}_{other}}i;\theta^i)$$
其中$\theta_i$是智能体$i$的神经网络$f$的参数，包括一个softmax层输出策略，一个线性层输出值函数，所有非输出层是共享的。动作是从策略$\pi$中采样得到的。观察状态$s_{self}^i$包含$fi$智能体的位置，以及其他智能体的位置。每个智能体都有两个网络（为了简洁，省略了智能体上标$i$），一个计算它自己的动作和值函数，一个计算其他智能体的估计值，如下：
\begin{equation}
f_{self}(s_{self},z_{self},\hat{z}{other};\theta{self})
\end{equation}
\begin{equation}
f_{other}(s_{other},\hat{z}{other},z{self};\theta_{self})
\end{equation}
这两个网络使用的方式不同：$f_{self}$用于计算智能体自己的行为和价值，并以前馈方式运行。给出其他智能体观察到的动作，智能体使用$f_{other}$通过优化$\hat{z}{other}$推断其他智能体的目标。
我们建议每个智能体使用自己的策略模拟其他玩家的行为，这样$f{other}$的参数与$f_{self}$的参数是相同的。但请注意，两个网络的输入$z_{self}$和$\hat{z}{other}$的相对位置不同。另外，由于环境是完全可观测的，两个智能体的观察状态的不同仅通过地图上智能体的身份体现出来（即，每个智能体将能够区分其自己的位置和另一个智能体的位置）。因此，在acting模式下，$f{self}$网络将$s_{self}$作为输入；在推理模式下，$f_{other}$网络将$s_{other}$作为输入。在游戏的每一步，智能体需要推理$\hat{z}{other}$将其作为(1)的输入并选择其动作。为了实现这个目的，在每一步中，智能体观察另一个智能体采取的行动，并且在下一步中，智能体使用先前观察到的另一个智能体的动作作为监督信号，使用式子(2)反向传播并优化其$\hat{z}{other}$，如图1所示。
推理过程优化器中采取的步数是一个可根据游戏的不同而变化的超参数。因此，在游戏的每一步中其他智能体的目标估计$\hat{z}{other}$会被更新多次。参数$\theta{self}$在每个回合结束时使用和带有智能体获得的奖励信号的Asynchronous Advantage Actor-Critic（A3C）进行更新。
算法1给出了一个回合内训练SOM智能体的伪代码。这里考虑的所有任务的目标都是离散的，智能体的目标$\hat{z}{self}$被
表示独热向量，维度是智能体目标所有可能的情况数。另一个玩家的目标嵌入$\hat{z}{other}$有相同的维度。为了估计经过离散而不可微的变量$\hat{z}{other}$的梯度，我们用Gumbel-Softmax分布上的一个可微样本$\hat{z}{other}^G$代替它。这种重新参数化技巧被证明可以有效地产生低方差偏置的梯度。使用该方法在每一步优化过$\hat{z}{other}$之后，$\hat{z}{other}$通常偏离独热向量。在下一步中，$f_{self}$将对应于先前更新的$z_{other}$ argmax的一个独热向量量$\hat{z}{other}^OH$作为输入。
智能体的策略由长短期记忆（LSTM）单元参数化，以及两个全连接的线性层和指数线性单元（ELU）激活函数。神经网络的权重用半正交矩阵初始化。
由于$f{other}$的循环性，当推理步数$\gt 1$时必须特别小心。在这种情况下，在游戏的每一步中，我们在推理模式中的第一次前向传播之前保存$f_{other}$的循环状态，并且在每个推理步骤将循环状态初始化为此值。这个过程可以确保在动作和推理模式下$f_{other}$可以展开相同数量的步骤。

实验

在本节中，我们在三个任务上评估SOM模型：

硬币游戏，这是一个完全合作的任务，智能体的角色是对称的。
配方游戏，它是对抗的，但具有对称角色。
门禁游戏，它是完全合作的，但是两个玩家拥有不对称的角色。

我们将SOM与其他三个baselines以及一个可以访问其他智能体目标的ground truth的模型进行比较。所有任务都是在Mazebase gridworld环境中创建的。

Baselines

TRUE-OTHER-GOAL（TOG）：我们提供了一个给出的模型性能上限的策略网络，该网络将其他智能体的真正目标$z_{other}$，以及状态特征$s_{self}$和自己的目标$z_{self}$作为输入。因为这个模型可以直接访问其他智能体真正的目标，因此不需要单独的网络来模拟其他智能体的行为。 TOG的结构与SOM的一个策略网络$f_{self}$相同。 NO-OTHER-MODEL（NOM）：我们使用的第一个baseline仅使用观察状态$s_{self}$和自身目标$z_{self}$作为输入。NOM与SOM的一个策略网络$f_{self}$有相同的架构。该baseline没有对其他智能体的显式建模或估计它们的目标。
集成-策略-预测器（IPP）：从NOM的体系结构和输入开始，我们构建了一个更强的baseline IPP，它有一个额外的最终线性层输出另一个智能体下一个动作的概率分布。除了用于训练该网络策略的A3C损失函数，我们还添加交叉熵损失项训练其他智能体的行为的预测。
分离-策略-预测器（SPP）：He等人提出了一个基于DQN的对手建模框架。在他们的方法中，给定对手特有的人工提取的状态信息，训练一个神经网络预测对手的动作。该网络的中间隐藏表示用作Q网络的输入。
我们修改了He等人的模型应用到本文的场景中。特别的，我们使用A3C而不是DQN，我们不使用特定领域的特征表示对手的隐藏状态。
最后产生的SPP模型由两个独立的网络组成，一个策略网络用于决定智能体的动作，一个对手网络用于预测其他智能体的动作。对手网络将世界状态$s$和自己的目标$z_{self}$作为输入，并输出其他智能体在下一步采取动作的概率分布，以及其隐藏状态（由网络的循环给出）。与IPP一样，我们使用其他智能体的真实动作训练对手策略预测器的交叉熵损失。在每一步中，该网络输出的隐藏状态以及智能体观察状态和智能体自身的目标被作为智能体的策略网络的输入。策略网络和对手策略预测器都是与SOM结构相同的LSTM网络。
与SOM作对比，SPP没有显式推断出其他智能体的目标。相反，它通过预测智能体在每个时间步的动作来隐式的构建对手模型。在SOM中，一个参考的目标作为策略网络的附加输入。而在SPP，类似的参考目标是从对手策略预测器得到的隐藏表示，把它作为策略网络的附加输入。
训练细节。在我们的所有实验中，我们使用系数为$0.01$的熵，价值损失系数为$0.5$，折扣系数为$0.99$的A3C训练智能体的策略。使用Adam优化智能体商策略的参数，其中$\beta_1= 0.9,\beta_2= 0.999,\epsilonn =1\times 10^{-8}$，权重衰减为$0$。学习率为$0.1$的SGD用于推断另一个智能体的目标，$\hat{z}_{other}$。
硬币和食谱游戏中策略网络的隐藏层维度为$64$，门游戏中为$128$。所有游戏和模型的学习率都是$1\times 10^{-4}$。
观测状态$s$用一些独热向量表示，包括环境中所有物体的位置，以及智能体和另一个智能体的位置。这个输入状态的维度是$1\times n$特征，其中Coin，Recipe和Door游戏的特征数分别为$384$,$192$和$900$。对于每个实验，我们使用5个不同的随机种子训练模型。除非特殊说明，否则论文中展示的所有游戏结果都是每步进行的10次优化更新的结果。

讨论

在本文中，我们介绍了一种新方法，用于从其他智能体的行为中推断他们的隐藏状态，并使用这些估计来选择动作。我们证明了智能体能够在合作和竞争环境中估计其他参与者的隐藏目标，这使他们能够收敛到更好的政策并获得更高的回报。在本文提出的任务中，对其他智能体的显式建模比仅仅考虑其他代理成为环境的一部分更好的性能。 SOM的一个限制是它比其他baseline需要更长的训练时间，因为我们在每一步都进行了反向传播。但是，它的online更新方式对于适应环境中其他智能体的动作变化至关重要。SOM的一些主要优点是简单性和灵活性，它不需要任何额外参数来模拟环境中的其他代理，可以使用任何强化学习算法进行训练，并且可以轻松地与任何策略参数化或网络结构集成。SOM可以适应具有两个以上智能体的环境，因为智能体可以使用自己的策略来模拟任意数量的智能体的动作并推断其目标。而且，它可以很容易地推广到许多不同的环境和任务。
我们计划通过评估更复杂环境中的模型来扩展这项工作，包括两个以上的参与者，混合策略，更多样化的智能体类型（例如具有不同动作空间的智能体，奖励函数，角色或策略），以及假设其他玩家和自己一样的模型偏差。
未来研究的其他重要途径是设计能够适应环境中其他智能体非平稳策略的模型，处理具有分层目标的任务，并在测试时遇到新智能体时表现良好。
最后，许多研究领域可以从拥有其他智能体的模型中受益，这些智能体能够推理其他智能体的意图并预测他们的动作。这些模型可能对人机或师生互动，以及价值对齐问题有恒大帮助。此外，这些方法可用于多智能体任务中基于模型的强化学习，因为前向模型的准确性很大程度上取决于预测其他智能体动作的能力。

Policy Gradient With Value Function Approximation For Collective Multiagent Planning

发表于 2019-01-26 | 更新于 2019-12-17 | 分类于强化学习

摘要

分布式的部分可观测马尔科夫决策过程(Dec POMDP)为解决多智能体系统中的序列决策问题提供了一个框架。考虑到POMDP的计算复杂度，最近的研究主要集中在Dec-POMDP中一些易于处理但是比较实用的子问题。本文解决的就是其中的一个子问题叫做CDec-POMDP其中一系列智能体的共同行为影响了它们公共的reward和环境变化。本文的主要贡献是提出了一个actor-critic(AC)强化学习算法优化CDec-POMDP问题的policy。普通的AC算法对于大型问题收敛的很慢，为了解决这个问题，本文展示了如何将智能体的估计动作值函数进行分解从而产生有效的更新以及推导出一个基于局部奖励信号的新的critic训练方式。通过在一个合成的benchmark以及真实的出租车车队优化问题上和其他方法进行对比，结果表明本文的AC方法提供了比之前最好的方法还要高质量的方法。

引言

近些年来，分布式的部分可观测马尔科夫决策过程已经发展成了解决多智能体协作的序列决策问题的一个很有前景的(promising)方法。Dec-POMDP对智能体基于环境和其他智能体的不同部分观测最大化一个全局的目标进行建模。Dec-POMDP的具体应用包括协调行星探测，多机器人协调控制以及无线网络的吞吐量优化。然而，解决分布式的部分马尔科夫决策过程是相当困难的，即使对于只有$2$个智能体的问题呢是NP难的。
为了增大规模和提高真实问题中的应用，过去的研究已经探索了智能体之间严格的交互，如状态转换和观测独立，事件驱动的的交互以及智能体之间的弱耦合性。最近，一系列工作开始关注于智能体的身份不影响它们之间的交互上，环境的变化主要受到智能体的共同影响，和著名的阻塞游戏很像。一些城市交通中的问题如出租车调度可以用这样的协同规划模型进行建模。
在本文中，作者着重于集中的Dec-POMDP框架将一类不确定情况下的集中多智能体序列决策问题形式化。Nguyen等人提出了一个采样方法优化CDec-POMDP模型中的policy。之前方法的一个主要缺点是policy是用表格形式展现的，随着智能体的observation spaces改变时，表格形式的policy不能很好的进行扩展。受到最近一些强化学习工作的启发，本文的贡献是一个AC框架的强化学习算法用来优化CDec-POMDP的policy。Policy用函数如神经网络来表示可以避免表格形式的policy的扩展性问题。我们推导出了策略梯度并且基于CDec-POMDP中智能体的交互提出了一个估计的因子动作值函数。普通的AC算法因为学习全局reward的原因，在解决大型多智能体系统问题时收敛的很慢。为了解决这个问题，本文提出了一种新的方式去训练critic，高效利用智能体的局部值函数的估计动作值函数。
我们在一个合成的多机器人导航领域和现实世界中一个亚洲城市的出租车调度问题上测试了本文的方法，结果展示了本文的方法可以扩展到大型多智能体系统上。根据经验，我们因式AC方法比以前最好的方法给出的解决方案都要好。因式AC方法收敛的也比普通的AC方法快很多，验证了我们提出的critic训练方法的有效性。
相关工作 我们的工作基于具有近似值函数的策略梯度框架。然而，根据以往的经验显示，直接应用原始的策略梯度到多智能体任务中，尤其是CDec-POMDP模型中会产生较高方差。在本文中，我们展示了一个和CDec-POMDP兼容的估计值函数，它能产生高效且低方差的策略梯度更新。Peshkin很早之前就研究过了应用于分布式policy的强化学习，Guestrin还提出使用REINFORCE从协调图中训练一个因子值函数的softmax策略。然而，这些以前的工作中，策略梯度都是从全局的经验回报而不是分解后的critic中估计的。我们在第四章中展示了一个分解后ciritc和基于训练这个critic得到的一个单个值函数对于高效的采样学习是很重要的。我们的实验结果表明了我们提出的critic训练方式比用全局经验回报训练收敛的还要快。

集中分布式POMDP模型

我们首先介绍一下Nguyen提出的CDec-POMDP模型。一个对应于这个模型的$T$步的动态贝叶斯网络如图所示。它由以下几个部分组成：

一个有限的计划范围$H$
智能的数量$M$，一个智能体m可能处在state space $S$中的任意一个状态，联合state space是$\times_{m=1}^MS$，我们用$i\in S$表示一个state。
每一个智能体m都有一个action spaceA，我们用$j\in A$表示一个action。
用$(s_{1:H},a_{1:H})^m=(s_1m,a_1^m,\cdots,s_Hm,a_H^{m)$表示一个智能体m完整的state-action轨迹。用随机变量$s_t}m,a_t^m$表示智能体$m$在$t$时刻的state和action。不同的指示函数$I_t(\cdot)$如表$1$所示。给定每一个智能体$m\in M$的轨迹，定义以下的计数方式：
$$n_t(i,j,i’) = \sum_{m=1}^M I_t^m(i,j,i’),\forall i,i’\in S,j\in A.$$
如表$1$所示，计数器$n_t(i,j,i’)$表示在$t$时刻处于state $i$，采取action $j$，转换到state $i’$的智能体数量。其他计数器$n_t(i)$和$n_t(i,j)$的定义类似。使用这些计数器，我们可以定义$t$时刻的计数表$\bf{n}{s_t}$和$\bf{n}{s_ta_t}$如表$1$所示。
我们假设一个普遍的部分观测环境，其中智能体基于其他智能体的总体影响可以有不同的ovservation。一个智能体观测到它的局部state $s_t^{m$。此外在$t$时刻基于它的局部状态$s_t}m$和计数表$\bf{n}_{s_t}$观测到$o_t^m$。例如，一个智能体m在$t$时刻处于state $i$，可以观测到其他也处在state $i(=n_t(i))$的智能体或者其他处在state $i$临近状态$j$的智能体，即$n_t(j),\forall j\in Nb(i)$。
状态转换函数是$\Phi_t(s_{t+1}^m=i’|s_tm=i,a_t^m=j,\bf{N}{s_t})$。所有智能体的状态转换函数是一样的，注意它会受到$\bf{n}{s_t}$的影响，而$\bf{n}_{s_t}$依赖于智能体的共同行为。
每一个智能体m有一个不平稳的policy $\pi_t^m(j|i,o_tm(i,\bf{n}{s_t}))$，表示在$t$时刻给定智能体m的observation $(i,o_t^m(i,\bf{n}{s_t})$之后，智能体采取action $j$的概率。我们用$\pi^m=(\pi_1,\cdots,\pi_H)$表示智能体m水平范围的policy。
一个智能体接收到的reward $r_t^m=r_t(i,j,\bf{n}{s_t}$取决于它的局部state和action，以及计数表$\bf{n}{s_t}$。
初始的state分布，$b_o=(P(i)\forall i \in S)$，对于所有的智能体都是相同的。

我们在这里展示了最简单的版本，所有的智能体的类型都相同，并且有相似的state transition，observation和reward模型。模型也可以处理多种类型的智能体，不同类型的智能体有不同的变化。我们还可以引入一个不受智能体action影响的external state，如交通领域的出租车需求。我们的结果也可以扩展到解决类似的问题。
像CDec-POMDP之类的模型对于解决智能体数量很大或者智能体的身份不影响reward或者transition function之类的问题是很有用的。其中一个应用是出租车车队优化问题，这个问题是计算出出租车调度的policy使得车队的利润最大化。一个出租车的决策过程如下。在时刻$t$时，每个出租车观测到它当前的城市空间$z$，不同的空间构成了state space $S$，以及当前空间和它的相邻空间的其他出租车的计数和当前局部请求的一个估计。这构成了出租车基于计数的observation $o(\cdot)$。基于这个observation，出租车必须决定待在当前空间$z$寻找乘客还是移动到下一个空间。这些决策选择取决于不同的因子，如请求比率和当前空间其他出租车的计数。类似的，环境是随机的，在不同时间出租车请求是变化的。使用出租车车队的的GPS记录可以得到这些历史的请求数据。
基于计数的统计数据用于规划 CDec-POMDP模型的一个关键属性是模型的变换取决于智能体的集中交互而不是智能体的身份。在出租车车队优化问题中，智能体数量可以相当大（大约有$8000$个智能体在现实世界的实验中）。给出这么大数量的智能体个数，为每一个智能体计算出独一无二的policy是不可能的。因此，和之前的工作类似，我们的目标是对所有智能体计算出一个相同的policy $\pi$。因为policy $\pi$取决于计数，它代表了一种富有表现力的policy。
对于一个固定的数量M来说，用${(s_{1:T},a_{1:T})^m\forall m}$表示从图$1$的DBN网络中采样得到的不同智能体的state-action轨迹。用$\mathbf{n}{1:T}={(\mathbf{n}{s_t},\mathbf{n}{s_ta_t},\mathbf{n}{s_ta_ts_{t+1}})\forall t=1:T}$表示每一个时间步$t$的结果计数表的组合向量。Nguyen等人展示了计数器$\mathbf{n}$中拥有足够的统计数据用来规划。也就是说，一个policy $\pi$在水平范围H内的联合值函数可以通过计数器的期望进行计算：
$$V(\pi) = \sum_{m=1}^M\sum_{T=1}H E[r_T^m] = \sum_{\mathbf{n}\in \Omega_{1:H}}P(\mathbf{n};\pi) \left[\sum_{T=1}^H\sum_{i\in S,j\in A} n_T(i,j)r_T(i,j,\mathbf{n}T)\right]$$
集合$\Omega{1:H}$是所有允许的一致计数表的集合，如下所示：
$$\sum_{i\in S}n_T(i) = M \forall T;$$
$$\sum_{j\in A}n_T(i,j) = n_T(i) = \forall j \forall T;$$
$$\sum_{i’\in S}n_T(i,j,i’) = n_T(i,j)\forall i\in S,\forall j \in A, \forall T;$$
$P(\mathbf{n},\pi)$是计数器的分布。这个结果的一个关键好处是我们可以直接从分布$P(\mathbf{n})$中对计数器$\mathbf{n}$采样而不是对单个不同智能体的轨迹$(s_{1:H},a_{1:H})进行采样来$评估policy $\pi$，这显著节省了计算开销。我们的目标是计算最优的policy $\pi$来最大化$V(\pi)$。我们假设一个集中式学习，分布式执行的强化学习设置。我们假设有一个模拟器可以从$P(\mathbf{n};\pi)$中提供计数器样本。

CDec-POMDP的策略梯度

之前的工作提出了一个基于采样的EM算法来优化policy $\pi$。这个policy被表示成计数器$\mathbf{n}$空间中的一个线性分段表policy，其中每一个线性片段指定了下一个action的分布。然而，这种表格形式的表示限制了它的表达能力，因为片段的数量是固定的先验，并且每个范围都必须手动定义，这可能会对性能产生不利影响。此外，当observation o是多维的时候，即，一个智能体观测到它位置相邻区域的计数器时，需要指数多个片段。为了解决这个问题，我们的目标是优化函数形式（如神经网络）的policy。
我们首先扩展策略梯度理论到CDec-POMDP上，用$\theta$表示policy参数的向量。我们接下来展示如何计算$\Delta_\theta V(\pi)$。用$\mathbf{s}_t,\mathbf{a}t$表示$t$时刻所有智能体的联合state和联合action。给定一个policy $\pi$，值函数表示形式如下：
$$V_t(\pi)=\sum{\mathbf{s}_t,\mathbf{a}t}P^{{\pi}(\mathbf{s}_t,\mathbf{a}_t|b_o,\pi)Q_t}{\pi}(\mathbf{s}t,\mathbf{a}T)$$
其中$P^{{\pi}(\mathbf{s}_t,\mathbf{a}_t|b_o)=\sum_{\mathbf{s}_{1:t-1},\mathbf{a}_{1:t-1}}P}{\pi}(\mathbf{s}{1:t},\mathbf{a}{1:t}|b_o)$是policy $\pi$下联合state $\mathbf{s}t$，和联合action $\mathbf{a}t$的分布。值函数$Q_t^{\pi}(\mathbf{s}t,\mathbf{a}t)$的计算过程如下：
$$Q_t^{\pi}(\mathbf{s}t,\mathbf{a}t) = r_t(\mathbf{s}t,\mathbf{a}t)+\sum{\mathbf{s}{t+1},\mathbf{a}{t+1})}P^{{\pi}(\mathbf{s}_{t+1},\mathbf{a}_{t+1}|\mathbf{s}_t,\mathbf{a}_t))Q_{t+1}}{\pi}(\mathbf{s}{t+1},\mathbf{a}{t+1})$$
接下来介绍以下CDec-POMDP的策略梯度理论：
定理1. 对于任何CDec-POMDP，策略梯度计算公式如下：
$$\Delta{\theta}V_1(\pi)=\sum{t=1}^{HE_{\mathbf{s}_t,\mathbf{a}_t)|b_o,\pi}\left[Q_t}{\pi}(\mathbf{s}t,\mathbf{a}t)\sum{i\in S,j\in A}n_t(i,j)\Delta{\theta}log\pi{t}(j|i,o(i,\mathbf{n}{s_t}))\right]$$
这个定理的证明和其他后续结果在附录中。
注意由于许多原因利用上述结果计算策略梯度是不切实际的。联合state-action $\mathbf{a}_t,\mathbf{s}_t$空间是组合的。考虑到智能体的个数可能有很多个，对每一个智能体的轨迹进行采样是计算上不可行的。为了补救，我们接下来会展示类似policy评估直接对计数器$\mathbf{n}~P(\mathbf{n};\pi)$进行采样计算梯度。类似的，也可以使用经验回报作为动作值函数$Q_t^{{\pi}(\mathbf{s}_t,\mathbf{a}_t)$的一个近似估计。这是标准的REINFORCE算法在CDec-POMDP上的应用。众所周知，REINFORCE可能比其他使用学习的动作值函数的方法学习的慢。因此，我们提出了一个$Q_t}{\pi}$的近似函数，展示了直接采样计数器$\mathbf{n}$来计算策略梯度。

使用估计动作值函数的策略梯度

估计动作值函数$Q_t^{\pi}(\mathbf{s}t,\mathbf{a}t)$有几种不同的方式。我们考虑下列特征形式的近似值函数$f_w$：
$$Q_t^{\pi}(\mathbf{s}t,\mathbf{a}t)\approx f_w(\mathbf{s}t,\mathbf{a}t)=\sum{m=1}^Mf_wm(s_t^m,o(s_tm,\mathbf{n{s_t}}),s_t^m)$$
每一个智能体m都定义了一个$f_w^{m$，它的输入是智能体的局部state，action和observation。注意不同的$f_w}m$是相关的，因为它们依赖于公共的计数器表$\mathbf{n}{s_t}$。这样的一种分解方式是很有用的，因为它产生了有效的策略梯度计算方式。此外，CDec-POMDP中一类很重要的这种形式的估计值函数是兼容值函数最后会产生一个无偏的策略梯度。
命题1 CDec-POMDP中的兼容值函数可以分解成：
$$f_w(\mathbf{s}t\mathbf{a}t) = \sum_mf_w^m(s_tm,o(s_t^{m,\mathbf{n}_{s_t}),a}m)$$
我们可以直接用估计值函数$f_w$取代$Q^{\pi}(\cdot)$。经验上来说，我们发现使用这个估计的方差很大。我们利用$f_w$的结构进一步分解策略梯度会有更好的效果。
定理2 对于任何具有如下的分解的值函数：
$$f_w(\mathbf{s}t\mathbf{a}t) = \sum_mf_w^m(s_tm,o(s_t^{m,\mathbf{n}_{s_t}),a}m)$$
策略梯度可以写成：
$$\Delta{\theta}V_1(\pi)=\sum{t=1}^{HE_{\mathbf{s}_t,\mathbf{a}_t)|b_o,\pi}\left[\sum_m\Delta_{\theta}log\pi(a_t}m|s_t^m,o(s_tm,\mathbf{n}{s_t}))f_w^m(s_tm,o(s_t^{m,\mathbf{n}_{s_t}),a_t}m)\right]$$
上述结果展示了如果估计值函数被分解了，那么得到的策略梯度也是分解的。上述结果也可以应用到多种类型的智能体上，只要我们假设不同的智能体有不同的函数$f_t^m$。最简单的情况下，所有的智能体都是相同类型的，每一个智能体都有相同的函数$f_w$，推断出下式：
$$f_w(\mathbf{s}t,\mathbf{a}t) = \sum{i,j}n_t(i,j)f_w(i,j,o(i,\mathbf{n}{s_t}))$$
使用上式，我们可以将策略梯度简化成：
$$\Delta{\theta}V_1(\pi) = \sum_tE{\mathbf{s}t,\mathbf{a}t}\left[\sum{i,j}n_t(i,j)\Delta{\theta}log\pi (j|i,o(i,\mathbf{n}{s_t}))f_w(i,j,o(i,\mathbf{n}{s_t}))\right]$$

基于计数器的策略梯度计算

注意在上式中，期望仍然和联合state，action，$(\mathbf{s}t,\mathbf{a}t)$相关，当智能体的个数很大时效率很低。为了解决这个问题是
定理3 对于任何拥有形式$f_w(\mathbf{s}t,\mathbf{a}t) = \sum{i,j}n_t(i,j)f_w(i,j,o(i,\mathbf{n}{s_t}))$的值函数，策略梯度都可以用下式计算：
\begin{equation}
E{\mathbf{n}{1:H}\in \Omega_{1:H}} \left[\sum_{t=1}^H\sum_{i\in S,j\in A}n_t(i,j) \Delta_{\theta}log\pi (j|i,o(i,\mathbf{n}_t)) f_w(i,j,o(i,\mathbf{n}t))\right]
\end{equation}
上述结果展示了策略梯度可以类似于计算policy的值函数一样通过从底层分布$P(\cdot)$中采样计数表向量$\mathbf{n}{1:H}$来计算策略梯度，在智能体数量很大的情况下也是可行的。

训练动作值函数

在我们的方法中，在计数器样本$\mathbf{n}{1:H}$生成用来计算策略梯度后，我们还需要调整critic $f_w$的参数。注意对于每一个动作值函数$f_w(\mathbf{s}t,\mathbf{a}t)$只取决于联合state，action $(\mathbf{s}t,\mathbf{a}t)$生成的计数器。训练$f_w$可以通过一个梯度步最下化下列loss函数实现：
\begin{equation}
min_w\sum{\xi=1}^K\sum_{t=1}H\left(f_w(\mathbf{n}t^{\xi})-R_t{\xi}\right)^2
\end{equation}
其中$\mathbf{n}{1:H}^{{\xi}$是从分布$P(\mathbf{n};\pi)$中生成的一个计数器样本；$f_w(\mathbf{n}_t}{\xi})$是动作值函数，$R_t^{\xi}$是用式子$(1)$计算的$t$时刻的所有经验回报：
\begin{equation}
f_w(\mathbf{n}t^{\xi}) = \sum{i,j}n_t^{{\xi}(i,j)f_w(i,j,o(i,\mathbf{n}_t}{\xi});R_t^{{\xi}=\sum_{T=t}}H]\sum{i\in S,j\in A}n_T{\xi}(i,j)r_T(i,j,\mathbf{n}T^{\xi})
\end{equation}
然而，我们发现公式$(11)$中的loss函数在训练较大问题的critic时表现并不好。需要一定数量的计数器样本可靠的训练$f_w$，这对于拥有较多数量智能体的大问题的扩展有不利影响。已知在多智能体强化学习中单独利用全局reward信号的算法要比利用局部reward信号的方法多用一些样本。受到这些现象的启发，接下来我们提出了一个基于策略的局部reward信号去训练critic $f_w$。
单个值函数 用$\mathbf{n}{1:H}^{{\xi}$表示一个计数器样本。给定计数器样本$\mathbf{n}_{1:H}}{\xi}$，用$V_t^{{\xi}(i,j)=E\left[\sum_{t’=t}}Hr{t’}^m|s_tm=i,a_m^t=j,n_{1:H}{\xi}\right]$表示一个智能体在时刻$t$处于state $i$，采取action $j$，所能得到的所有期望reward。这个单个的值函数可以用动态规划算法来计算。基于这个值函数，我们接下来展示了式子$(12)$中全局经验reward的重新参数化：
引理(Lemma)1 给定计数器样本$\mathbf{1:H}^{{\xi}$，$t$时刻的经验回报$R_t}{\xi}$可以被重新参数化为：
$$R_t^{\xi} = \sum{i\in S,j\in A}n_t^{{\xi}(i,j)V_t}{\xi}(i,j).$$
基于单个值函数的loss 给出引理$1$，我们推导出式子$11$中真实loss的上界，它有效利用了单个值函数：
\begin{align*}
&\sum{\xi}\sum_t\left(f_w(\mathbf{n}^{\xi})-R_t{\xi}\right)^2 \
= &\sum_{\xi}\sum_t\left(\sum_{i,j}n_t^{{\xi}(i,j)f_w(i,j,o(i,\mathbf{n}_t}{\xi}))-\sum_{i,j}n_t^{{\xi}(i,j)V_t}{\xi}(i,h)\right)^2\
= &\sum_{\xi}\sum_t\left( \sum_{i,j}n_t^{{\xi}(i,j)(f_w(i,j,o(i,\mathbf{n}_t}{\xi}))-V_t^{{\xi}(i,h))\right)}2\
\le &M\sum_{\xi}\sum_{t,i,j}n_t(i,j)\left(f_w(i,j,o(i,\mathbf{n}_t^{\xi}))-V_t{\xi}(i,j)\right)^2
\end{align*}
其中最后一部用了柯西施瓦茨不等式。我们用式子(14)中修改过的loss训练critic。按照经验来说，对于较大的问题，式子(14)中的新loss比式子(13)中的原始loss要收敛的快很多。直观上来说，这是因为式子(14)中的新loss尝试调整每一个critic组件$f_w(i,j,o(i,\mathbf{n}_t^{{\xi}))$更接近它的经验回报$V_t}{\xi}(i,j)$。然而，原始的式子(13)中的loss着重于最小化全局loss，而不是调整每一个单个的critic因子$f_w(\cdot)$到相对应的每一个经验回报。
算法$1$展示了CDec-POMDP中AC算法的大纲。第$7$行和第$8$行展示了两种不同的方式训练critic。第$7$行代表基于局部值函数的critic更新，也可以称为factored cirtic更新(fC)。第$8$行展示了基于全局reward或者全局critic的更新©。第$10$行展示了使用定理$2$(fA)计算的策略梯度。第$11$行展示了直接使用$f_w$计算的梯度。

实验

这一节中比较了我们的AC算法和另外两个解决CDec-POMDP问题的算法，Soft-Max based flow update(SMFU)，和期望最大化方法。SMFU只能优化智能体的action依赖于局部state的policy，$\pi(a_t^m|s_tm)$，因为它通过计算在规划阶段单个最有可能的计数器向量来估计计数器$\mathbf{n}$的作用。EM方法优化基于计数器的分段线性policy，其中$\pi(a_t^m|s_tm,\cdot)$是所有可能的计数器observation $o_t$空间上的一个分段函数。
算法$1$展示了更新critic的两种方式（第$7$行和第$8$行）和更新actor的两种方式（第$10$行和第$11$行），所以就有四种可能的AC方法－fAfC,AC,FfC,fAC。我们也研究了不同actor-critic方法的属性。在附录中有神经网络的结构和其他一些实验设置。
为了和之前方法公平的进行比较，我们使用了三种不同的模型用于基于计数的observation $o_t$。在$o0$设置中，policy只取决于智能体的局部state $s_t^m$并不需要计数器。在$o1$设置中，policy取决于局部state $s_t^m$和单个计数器observation $n_t(s_t^{m)$。也就是说，智能体只能观测到其他也在当前状态$s_t}m$的智能体的计数器。在$oN$设置中，智能体能观测到它的局部state $s_t^{m$和当前状态$s_t}m$的局部相邻状态内其他智能体的计数器。$oN$ observation模型提供给智能体最多的信息。然而，它也是最难优化的因为policy有更多的参数。SMFU方法只能在$o0$设置中起作用，EM方法和本文中的AC方法在所有设置中都能起作用。

EM(Expectation Maximization)算法

发表于 2019-01-21 | 更新于 2019-12-17 | 分类于机器学习

引言(Introduction)

什么是期望最大化算法

期望最大化算法(Expectation Maximization,EM)，是利用参数估计的迭代法求解最大似然估计的一种方法。

EM和MLE关系

MLE的目标是求解已知分布类型的单个分布的参数。
EM的目标是求解已知分布类型的多个混合分布的参数。
一般我们用到的极大似然估计都是求某种已知分布类型的单个分布的参数，如求高斯分布的均值和方差；而EM算法是用来求解已知分布类型，多个该已知类型分布的混合分布的参数，这句话听起来可能有些拗口，举个最常见的例子，高斯混合分布参数的求解，这个混合分布都是高斯分布，只是每个分布的参数不同而已。如果一个高斯分布，一个卡方分布是没有办法求解的。

为什么叫它EM算法

因为这个算法总共有两个迭代步骤，E步和M步。第一步是对多个分布求期望，固定每一个分布的参数，计算出混合分布的参数，即E步，第二步是对这个混合分布利用最大似然估计方法进行参数估计，即M步。

推理过程

假设我们要求一个混合分布p的参数$\theta$，比如校园内男生和女生的身高参数，显然，男生和女生的身高服从的分布类型是相同的，但是参数是不一样的。这里通过引入一个隐变量$z$，求解出对应不同$z$取值的参数$\theta$的值。
\begin{align*}
p(x|\theta) &= \sum_zp(x,z|\theta)\\
&=\sum_zp(z|\theta)p(x|\theta, z) \tag{0}
\end{align*}
如果我们假设男女生的身高分布是一个高斯混合模型，现在要求它的参数$\theta$。混合模型的表达式可以写为：
\begin{align*}
p(x|\theta) &= \sum_zw(z)N(x|\mu_z,\sigma_z)\\
&=\sum_zp(z|\theta)p(x|\theta,z)
\end{align*}
其中$\sum_zw(z) = 1,\theta={w, \mu, \sigma}$，如果用最大似然估计来解该问题的话，log函数内有和式，不好优化，所以就要换种方法。
观测数据：$x=(x_1,\cdots, x_N)$
对应的隐变量：$z=(z_1,\cdots, z_N)$，$z_i$有$c$种取值。

\begin{align*}
l(\theta;x) &= log p(x|\theta) \tag{1}\\
&= log\prod_{i=1}^N\ p(x_i|\theta) \tag{2}\\
&= \sum_{i=1}^Nlog\ p(x_i|\theta) \tag{3}\\
&= \sum_{i=1}^Nlog\sum_zp(x_i,z|\theta) \tag{4}\\
\end{align*}
这里式子(4)中$\sum_zp(x,z|\theta)$该怎么变形，因为现在解不出来了。
最开始我想的是使用条件概率进行展开，即：
$$\sum_zp(x_i, z|\theta) = \sum_zp(x_i|z, \theta)p(z|\theta)$$
但是如果展开成这样子，就变成了文章开头给出的式子(0)，并没有什么用，不能继续化简了。
所以就对式子(4)做个变形
\begin{align*}
&\ \ \ \ \sum_{i=1}^Nlog\sum_zp(x_i,z|\theta) \tag{4}\\
&= \sum_{i=1}^Nlog\sum_zq(z|x_i)\frac{p(x_i,z|\theta)}{q(z|x_i)}, \ \ s.t.\sum_zq(z|x_i)=1 \tag{5}\\
&\ge \sum_{i=1}^N \underbrace{\sum_zq(z|x_i)log\frac{p(x_i,z|\theta)}{q(z|x_i)}}_{L(q,\theta)},\ \ s.t. \sum_zq(z|x_i)=1 \tag{6}\\
\end{align*}
第(4)步到第(5)步引入了一个分布$q(z|x)$，就是给定一个观测数据$x$，隐变量$z$取值的概率分布。注意，$q(z)$是一个函数，但是给定$x$之后，$q(z|x)$是一个变量。然后因为变形之后还是没有求解，就利用杰森不等式做了缩放，将$log(sum())$变成了$sum(log())$，就变成了(6)式。
这里使用Jensen不等式的目的是使得缩放后的值还能取得和原式相等的值，重要的是等号能够取到。

Jensen不等式

对于随机变量的Jensen不等式，当函数$f(x)$是凸函数的时候可以用下式表示：
$$f(E(x)) \le E(f(x))$$
当$f(x)$是凹函数的时候，有
$$f(E(x)) \ge E(f(x))$$

接下来我们就要求解使得式子(6)中杰森不等式等号成立的$q$分布的取值。这里有两种方法可以求解。

拉格朗日乘子法

令
$$L(q,\theta) = \sum_z q(z|x_i)log{\frac{p(x_i,z|\theta)}{q(z|x_i)}}, s.t.\sum q(z|x_i) = 1 \tag{7}$$
构建拉格朗日目标函数：
\begin{align*}
L &= L(q, \theta) + \lambda(\sum_zq(z|x)- 1) \tag{8}\\
&= \sum_z q(z|x_i)log{\frac{p(x_i,z|\theta)}{q(z|x_i)}} + \lambda(\sum_z q(z|x_i) - 1) \tag{9}
\end{align*}

对$L$求导，得到：
$$\frac{\partial L}{\partial q(z|x_i)} = log\frac{p(x_i, z|\theta)}{q(z|x_i)} + q(z|x_i)(-\frac{1}{q(z|x_i)}) + \lambda \tag{10}$$
令$\frac{\partial L}{\partial q(z|x_i)}$等于$0$，得到：$$log\frac{p(x_i, z|\theta)}{q(z|x_i)} = 1 - \lambda$$
两边同取$e$的对数：
$$\frac{p(x_i, z|\theta)}{q(z|x_i)} = e^{1-\lambda} \tag{11}$$
$$q(z|x_i) = e^{\lambda - 1}p(x_i, z|\theta) \tag{12}$$
两边同时求和得：
$$1 = e^{\lambda - 1}\sum_z p(x_i, z|\theta) \tag{13}$$
用$p$表示$e^{\lambda-1}$得到：
$$e^{\lambda-1} = \frac{1}{\sum_z p(x_i, z|\theta)}$$
将其代入式子(12)得：
\begin{align*}
q(z|x_i) &= \frac{p(x_i, z|\theta)}{\sum_z p(x_i, z|\theta)}\\
&= \frac{p(z, x_i|\theta)}{p(x_i|\theta)}\\
&= p(z|x_i, \theta) \tag{14}
\end{align*}

最后求出来$q(z|x_i) = p(z|x_i, \theta)$。

杰森不等式成立条件

杰森不等式成立条件是常数，即：
$$\frac{p(x_i, z|\theta)}{q(z|x_i)} = c, s.t. \sum q(z|x_i)=1 \tag{15}$$
则有:
$$p(x, z_i|\theta) = cq(z_i|x) \tag{16}$$
同时对式子左右两边求和，得到：
$$\sum p(x_i, z|\theta) = \sum cq(z|x_i) = c \tag{17}$$
将$c = \sum p(x_i, z|\theta)$代入式子(14)得：
\begin{align*}
q(z|x_i) &= \frac{p(x_i, z|\theta)}{\sum p(x_i,z|\theta)}\\
&= \frac{p(x_i, z)|\theta}{p(x_i|\theta)}\\
&= p(z|x_i, \theta) \tag{18}
\end{align*}

等号成立证明

上面两个方法都算出来在$q(z|x_i) = p(z|x_i, \theta)$时$L$能取得最大值。接下来证明这个这个$L$的最大值和$l$相等。
将$q = p(z|x_i, \theta)$代入$L(q, \theta)$得：
\begin{align*}
L(q, \theta) &= L(p(z|x_i, \theta^t), \theta^t)\\
&= \sum_z p(z|x_i, \theta^t) log\frac{p(z, x_i|\theta^t)}{p(z|x_i, \theta^t)} \\
&= \sum_z p(z|x_i, \theta^t) log p(x_i|\theta^t)\\
&= 1\cdot log p(x_i|\theta^t)\\
&= log p(x_i|\theta^t)\\
&= l(\theta^t; x_i)
\end{align*}

另一种等号成立推导

\begin{align*}
l(\theta; x) - L(q, \theta) &= l(\theta; x_i) - \sum_z q(z|x_i) log{\frac{p(z, x_i|\theta)}{q(z|x_i)}}\\
&= \sum_z q(z|x_i) log p(x_i|\theta) - \sum_z q(z|x_i) log{\frac{p(z, x_i|\theta)}{q(z|x_i)}}\\
&= \sum_z q(z|x_i)log {\frac{p(x_i|\theta)q(z|x_i)}{p(z, x_i|\theta)}}\\
&= \sum_z q(z|x_i)log {\frac{q(z|x_i)}{p(z|x_i, \theta)}}\\
&= KL(q(z|x_i)||p(z|x_i,\theta))
\end{align*}
最后算出来两个函数之差是一个KL散度，是从$p$到$q$的KL散度。当前仅当$p=q$时取等，否则就非负。

M步

\begin{align*}
L(q, \theta) & = \sum_z q(z|x_i) log\frac{p(z, x_i|\theta)}{q(z|x_i)} \\
& = \underbrace{\sum_z q(z|x_i)log{p(z, x_i|\theta)}}_{Expected\ complete\ log-likelyhood} - \underbrace{\sum_z q(z|x_i)l{q(z|x_i)}}_{Entropy}
\end{align*}

EM流程

计算流程

（１）首先随机初始化模型的不同隐变量对应的参数，
（２）对于每一个观测，首先判断它对应的隐变量的分布。
（３）求期望
（４）最大似然估计求参数
用公式来表示如下：
E步：$q^{t+1} = arg\ max_q L(q, \theta^t)$
M步：$\theta^{t+1} = arg max_{\theta}L(q^{t+1}, \theta)$
E步就是根据$t$时刻的$\theta^t$利用概率$q$求出$L$的期望，然后M步使用最大似然估计计算出新的$\theta$，就这样迭代下去。

EM收敛性分析

EM算法的收敛性就是要证明$L(q=p(z|x_i, \theta^t) , \theta)$的值一直在增大。
\begin{align*}
L(p(z|x_i, \theta^{t+1}) , \theta^{t+1}) - L(p(z|x_i, \theta^{t}) , \theta^{t}) &= log p(x_i|\theta^{t+1}) - log p(x_i|\theta^t)\\
& \ge 0
\end{align*}

例子

假如有两个硬币A和B，假设随机从A,B中选一个硬币，掷$10$次，重复$5$次实验，分别求出两个硬币正面向上的概率。假设硬币服从二项分布
$5$次实验结果如下：
5H5T
9H1T
8H2T
4H6T
7H3T

这个时候有两种情况

知道每次选的是A还是B

这个时候就变成了极大似然估计。

不知道每次选的是A还是B

这个时候就用EM算法了。
首先随机初始化$\theta_A = 0.5, \theta_B = 0.5$，
对于每一个观测，首先判断它对应的隐变量的分布。
$i={1,2,3,4,5}$，分别代表$5$个实验。
首先求出$\theta_A$的参数。
$$P(z = A|x_i, \theta_A, \theta_B) = \frac{P(z = A|x_i, \theta_A)}{P(z = A|x_i, \theta_A) + P(z = B|x_i, \theta_B)}$$
$$P(z = B|x_i, \theta_A, \theta_B) = 1 - P(z = A|x_i,\theta_A,\theta_B)$$
然后计算下式：
\begin{align*}
L(q,\theta_A) &= \sum_{i=1}^5 \sum_zp(z|x_i, \theta_A, \theta_B)log p(x_i|\theta)\\
&= \sum_{i=1}^5 (p(z=A|x_i, \theta_A)log p(x_i|\theta_A) + p(z=B|x_i, \theta_B)log p(x_i|\theta_B))
\end{align*}
然后利用极大既然估计计算$\theta_A$和$\theta_B$的值。

参考文献

1.https://www.zhihu.com/question/27976634/answer/153567695
2.https://en.wikipedia.org/wiki/Jensen's_inequality

maximum likelyhood estimation

发表于 2019-01-20 | 更新于 2019-12-17 | 分类于机器学习

最大似然估计

前提：数据的分布已知（如服从高斯分布或者指数分布），但是分布参数未知。例如某学校学生的身高服从高斯分布$N(\mu, \sigma^2 )$，$\mu$和$\sigma^2 $未知。现随机抽取$200$个学生的身高，估计该学校学生身高的均值$\mu$和方差$\sigma^2 $。
概率密度：$p(x|\theta)$
样本集：$X=(x_1, x_2,\cdots, x_N), N=200$，$x_i$为第$i$个人的身高。
假设分布：$N(\mu, \sigma^2 )$，$\mu, \sigma^2 $未知

联合概率密度函数

样本集中的$N$个样本是独立同分布的，它们的联合概率可以表示为：
$$L(\theta; X) = L(x_1, \cdots, x_n;\theta) = \prod_{i=1}^{N} p(x_i|\theta), \theta \in \Theta$$

取对数

因为$L$中包含乘法，不方便求导，可以对其取log，不改变函数的单调性，并且方便计算：
$$ln\ L(\theta; X) = ln\ L(\theta; x_1, \cdots, x_n) = \sum_{i=1}^N ln\ p(x_i|\theta), \theta \in \Theta$$

求极值

计算偏导数，令其等于$0$，取函数的极值点，因为只有一个极值点，所以一定是最大值点，即
$$\hat{\theta} = arg\ max_{\theta} ln\ L(\theta; x)$$
求偏导等于$0$即：
$$\frac{\partial ln\ L(\theta; X)}{\partial \theta} = \sum_{i=1}^N \frac{\partial ln\ p(x_i;\theta)}{\partial \theta}, \theta={\mu, \sigma^2}$$

最大似然估计求解高斯分布

若$p$为高斯分布，即$p(x; \theta) = \frac{1}{\sqrt{2\pi} \sigma} e^{-\frac{(x-\mu)^2 }{2 \sigma^2 } } $，$ln\ p(x;\theta) = -ln\ \sqrt{2\pi } - ln\ \sigma - \frac{(x_i-\mu)^2 }{ 2\sigma^2 }$
则：
$$ln\ L(\theta; X) = ln\ L(\theta; x_1, \cdots, x_n) = \sum_{i=1}^N ln\ p(x_i;\theta) = \sum_{i=1}^N \left(-ln\ \sqrt{2\pi} - ln\ \sigma - \frac{(x_i-\mu)^2 }{ 2\sigma^2 }\right), \theta \in \Theta$$

求解$\mu$

对$ \mu $ 求偏导得：
$$\frac{\partial ln\ L(\mu; X)}{\partial \mu} = \sum_{i=1}^N \frac{\partial ln\ p(x_i;\mu)}{\partial \mu} = \sum_{i=1}^N -\frac{2(x_i-\mu) }{2\sigma^2 } $$
令其等于$0$，即
$$\sum_{i=1}^N(x_i-\mu) = 0$$
解得$$ \mu = \frac{\sum_{i=1}^N x_i }{N}$$

求解$\sigma^2 $

对$ \sigma^2 $求偏导，令$t=\sigma^2 $，则$ln\ p(x;\sigma^2 ) = ln\ p(x;t) = -ln\ \sqrt{2\pi} - ln\ \sqrt{t} - \frac{(x_i-\mu)^2 }{2t}$，有：
$$ \frac{\partial ln\ L(t; X)}{\partial t} = \sum_{i=1}^N \frac{\partial ln\ p(x_i;t)}{\partial t} = \sum_{i=1}^N\left( -\frac{1}{2t} + \frac{(x-\mu)^2 }{2t^2 }\right) $$
令其等于$0$，即
$$ \sum_{i=1}^N\left( -\frac{1}{2t} + \frac{(x-\mu)^2 }{2t^2 }\right) = 0 $$
解得：
$$ \sigma^2 = t = \frac{\sum_{i=1}^N (x_i) -\mu)^2 }{N} $$

参考文献

bayesian classifier bayesian networks

发表于 2019-01-06 | 更新于 2019-12-17 | 分类于机器学习

概述

朴素贝叶斯分类假设属性之间是条件独立的，而在实际中，属性与属性之间难免会有一定的相关性，而不是完全独立，这时就不能用朴素贝叶斯分类器了。
贝叶斯网络不要求给定类的所有属性都条件独立，而是允许一些属性条件独立，一些属性不条件独立，它可以表现独立和条件独立之间的关系，然后根据这些关系计算后验概率。
概括的来说，本文的内容可以分为以下几个部分：

首先定义了贝叶斯网络的语法(syntax)和语义(semantics)，并且展示了如何用贝叶斯网络表示不确定知识。贝叶斯网络是一个有向无环图，节点代表随机变量，边代表因果关系。贝叶斯公式的两个意义，一个是数值意义，用贝叶斯网络表示全概率分布，另一个是拓扑意义，给定某个节点的父节点，这个节点条件独立于所有它的非后裔节点，或者给定某个节点的马尔科夫毯，这个节点条件独立于所有其他节点。
条件独立的有效表示，噪音或模型表示离散型父节点和离散型子节点之间的关系，用参数化模型表示连续型父节点和连续型子节点之间的关系，用probit模型或者logit模型表示连续型父节点和离散型子节点之间的关系。
贝叶斯精确推理计算后验分布的两种方法，一种是枚举推理，一种是消元法。
因为精确推理的复杂度太高了，没有实际应用价值，给出了估计推理的方法，直接采样，拒绝采样，以及可能性加权，还有另一类采样方法，蒙特卡洛算法，主要介绍了吉布森采样，大概就是这些。

贝叶斯网络的定义

我们可以根据联合概率分布(full joint probability distribution)算出任何想要的概率值（边缘分布什么的），但是随着随机变量个数的增加，联合概率分布可能会变得特别大。此外，一个一个的指定可能世界中的概率是不可行的。
如果在联合概率中引入独立和条件独立，会显著的减少定义联合概率分布所需要的概率。所以这节就介绍了贝叶斯网络来表示变量之间的依赖关系。本质上贝叶斯网络可以表示任何联合概率分布，而且在很多情况下是非常精确地表示。一个贝叶斯网络是一个有向图，图中的节点包含量化后的概率信息。具体的说明如下：

每一个节点对应一个随机变量，这个随机变量可以是离散的也可以是连续的。
有向边或者箭头连接一对节点。如果箭头是从节点$X$到节点$Y$，那么节点$X$称为节点$Y$的父节点。图中不能有环，因此贝叶斯网络是一个有向无环图(directed acyclic graph,DAG)。
每一个节点$X_i$有一个条件概率分布$P(x_i|Parents(X_i))$量化(quantifiy)父节点对其影响。

网络的拓扑，即节点和边的集合，指定了条件概率分布之间的关系。箭头的直观意义是节点$X$对节点$Y$有直接的影响，$Y$发生的原因是其父节点的影响。通常对于一个领域(domain)的专家来说，指出该域受哪些因素的直接影响要比直接给出它的概率值简单的多。一旦贝叶斯网络的拓扑结构定了，给出一个变量的父节点，我们仅仅需要给出每个节点的条件概率分布。我们能看出，拓扑和条件概率的组合能计算出所有变量的联合概率分布。

贝叶斯网络的示例

牙疼和天气

给定一组随机变量牙疼(Toothache)，蛀牙(Cavity)，拔牙(Catch)和天气(Weather)。Weather是独立于另外三个随机变量的，此外，给定Cavity，Catch和Toothache是条件独立的，即给定Cavity，Catch和Toothache是相互不受影响的，如下图所示。正式的：给定Cavity，Toochache和Catch是条件独立的，图中Toothache和Catch之间缺失的边体现出了条件独立。直观上，网络表现出Cavity是Toothache和Catch发生的直接原因，然而在Toothache和Catch之间没有直接的因果关系。
figure 14.1

警报和打电话

我家里有一个新安装的防盗警报(burglar alarm)，这个警报对于小偷的检测是相当可靠的，但是也会对偶然发生的微小的地震响应。我有两个邻居(Mary和John)，他们听到警报后会打电话给我。John有时会把电话铃和警报弄混了，也会打电话。Mary听音乐很大声，经常会错过警报。现在给出John或者Mary谁是否打电话，估计警报响了的概率。
figure 14.2
该例子的贝叶斯网络如上图所示。该网络体现了小偷和地震两个因素会直接影响警报响的概率，但是John和Mary会不会打电话只取决于警报有没有响。贝叶斯网络展示出了我们的假设，即John和Mary不直接观察小偷有没有来，也不直接观察小的地震是否，也不受之前是否打过电话的影响。上图中的条件概率分布以一个条件概率分布表(conditional probability table,CPT)的形式展现了出来。这个表适合离散型的随机变量，但是不适合连续性随机变量。没有父节点的节点只有一行，用来表示随机变量的可能取值的先验概率(prior probabilities)。
注意到这个网络中没有节点对应Mary听音乐很大时，也没有节点对应John把电话铃声当成了警报。事实上这些因素都被包含在和边Alarm到JohnCalls和MaryCall相关的不确定性中了，概率包含了无数种情况可能让警报失灵（停电，老鼠咬坏了，等等）或者John和Mary没有打电话的原因（吃饭去了，午睡了，休假了等等），这些不确定性都包含在了概率中了。

贝叶斯网络的意义

有两种方式可以理解贝叶斯网络的意义。第一个是一种数值化的意义(numerical semantics)，把它当成联合概率分布的一种表示形式。第二个是一种拓扑的意义(topological semantics)，将它看成条件独立的一种编码方式。事实上，这两种方式是等价的，只不过第一种方式有助于理解如何构建贝叶斯网络，而第二种方式更有助于设计推理过程。

贝叶斯网络表示联合分布(Representing the full joint distribution)

定义

一个贝叶斯网络是一个有向无环图，并且每个节点都有一个数值参数。数值方式给出这个网络的意义是，它代表了所有变量的联合概率分布。之前说过节点上的值代表的是条件概率分布$P(X_i|Parents(X_i)$，这是对的，但是当赋予整个网络意义以后，这里我们认为它们只是一些数字$\theta(X_i|Parents(X_i)$。
联合概率中的一个具体项(entry)表示的是每一个随机变量取某个值的联合概率，如$P(X_1=x_1 \wedge \cdots\wedge X_n = x_n)$，缩写为$P(x_1,\cdots,x_n)$。这个项的值可以通过以下公式进行计算：
$$P(x_1,\cdots,x_n) = \prod_{i=1}^n \theta(x_i|parents(X_i)),$$
其中$parents(X_i)$表示节点$X_i$在$x_1,\cdots,x_n$中的父节点。因此，联合概率分布中的每一项都可以用贝叶斯网络中某些条件概率的乘积表示。从定义中可以看出，很容易证明$\theta(x_i|parents(X_i))$就是条件概率$P(x_i|parents(X_i))$，因此，我们可以把上式写成：
$$P(x_1,\cdots,x_n) = \prod_{i=1}^n P(x_i|parents(X_i)),$$
换句话说：根据上上个式子定义的贝叶斯网络的意义，我们之前叫的条件概率表真的是条件概率表。（这句话。。。）

示例

我们可以计算出警报响了，但是没有小偷或者地震发生，John和Mary都打电话了的概率。即计算联合分布$P(j,m,a,\neg b, \neg e)$（使用小写字母表示变量的值）：
\begin{align*}
P(j,m,a,\neg b, \neg e) &=P(j|a)P(m|a)P(a|\neg b \wedge \neg e)P(\neg b)P(\neg e)\\
&=0.90\times 0.70\times 0.001 \times 0.999 \times 0.998\\
&=0.000628
\end{align*}

构建贝叶斯网络(Constructing Bayesian networks)

上面给出了贝叶斯网络的一种意义，接下来给出如何根据这种意义去构建一个贝叶斯网络。确定的条件独立可以用来指导网络拓扑的构建。首先，我们把联合概率的项用乘法公式写成条件概率表示：
$$P(x_1,\cdots,x_n) = P(x_n|x_{n-1},\cdots,x_1)P(x_{n-1},\cdots,x_1)$$
接下来重复这个过程，将联合概率(conjunctive probability)分解成一个条件概率和一个更小的联合概率。最后得到下式：
\begin{align*}
P(x_1,\cdots,x_n) &= P(x_n|x_{n-1},\cdots,x_1)P(x_{n-1}|,x_{n-2}\cdots,x_1)\cdots P(x_2|x_1)P(x_1)\\
&= \prod_{i=1}^nP(x_i|x_{i-1},\cdots,x_1)
\end{align*}
这个公式被称为链式法则，它对于任意的随机变量集都成立。对于贝叶斯网络中的每一个变量$X_i$，如果给定$Parents(X_i) \subset {X_{i-1},\cdots,X_1}$（每一个节点的序号应该和图结构的偏序结构一致），那么有：
$$P(x_1,\cdots,x_n) = \prod_{i=1}^n P(x_i|parents(X_i)),$$
将它和上式对比，得出：
$$P(X_i|X_{i-1},\cdots,X_1) = P(X_i|Parents(X_i).$$
这个公式成立的条件是给定每个节点的父节点，它条件独立于所有它的非父前置节点。这里给出一个生成贝叶斯网络的方式：

节点：首先，确定需要对领域建模所需要的随机变量集合。对它们进行排序：${X_1,\cdots,X_n}$，任意顺序都行，但是如果随机变量的因(causes)在果(effects)之前，最终的结果会更加紧凑。
边：从$i = 1$到$n$，

从$X_1,\cdots,X_{i-1}$中选出$X_i$的最小父节点集合。
对于每一个父节点，插入一条从父节点到$X_i$的边。
写下条件概率表，$P(X_i| Parents(X_i))$。

直观上，$X_i$的父节点应该包含$X_1,\cdots,X_{i-1}$中所有直接影响$X_i$的节点。因为每一个节点都只和它前面的节点相连，这就保证了每个网络都是无环的(acyclic)。此外，贝叶斯网络还不包含冗余的概率值，如果有冗余值，就会产生不一致：不可能生成一个违反概率论公理的贝叶斯网络。

紧凑性和节点顺序(Compactness and node ordering)

紧凑性(compactness)

因为不包含冗余信息，贝叶斯网络会比联合概率分布更加紧凑，这让它能够处理拥有很多变量的任务。贝叶斯网络的紧凑性是稀疏(sparse)系统或者局部结构化(local structured)系统普遍拥有的稀疏性的一个例子。在一个局部结构化系统中，每一个子部件仅仅和有限数量的其他部件进行交互，而不用管整个系统。局部结构化的复杂度通常是线性增加的而不是指数增加的。在贝叶斯网络中，一个随机变量往往最多受$k$个其他随机变量直接影响，这里的$k$是一个常数。为了简化问题，我们假设有$n$个布尔变量，指定一个条件概率表所需要的数字最多是$2^{k$个，整个网络则需要$n2}k$个值；作为对比，联合概率分布需要$2^n$个值。举个例子，如果我们有$n=30$个节点，每一个节点至多有五个父节点(k=5)，那么贝叶斯网络只需要$960$个值，而联合概率分布需要超过十亿个值。
但是在某些领域，可能每一个节点都会被所有其他节点直接影响，这时候网络就成了全连接的网络(fully connected)，它和联合概率分布需要同样多的信息。有时候，增加一条边，也就是一个依赖关系，可能会对结果产生影响，但是如果这个依赖很弱(tenuous)，添加这条边的花费比获得的收益还要大，那么就没有必要加这条边了。比如，警报的那个例子，如果John和Mary感受到了地震，他们认为警报是地震引起的，所以就不打电话了。是否添加Earthquake到JohnCalls和MaryCalls这两条边取决于额外的花费和得到更高的警报率之间的关系。

节点顺序(node ordering)

即使在一个局部结构化的领域，只有当我们选择好的节点顺序的时候，我们才能得到一个紧凑的贝叶斯网络。考虑警报的例子，我们给出下图：

Figure 14.2和Figure 14.3两张图中的三个贝叶斯网络表达的都是同一个联合分布，但是Figure 14.3中的两张图没有表现出来条件独立，尤其是Figure 14.3(b)中的贝叶斯网络，它需要用和联合分布差不多相同个数的值才能表现出来。可以看出来，节点的顺序会影响紧凑性。

贝叶斯网络中的条件独立 (Conditional independence relations in Bayesion networks)

贝叶斯网络的一个数值意义(“numerical” semantics)是用来表示联合概率分布。根据这个意义，给定每个节点的父节点，使得每一个节点条件独立于它的父节点之外的节点，我们能构建一个贝叶斯网络。此外我们也可以从用图结构编码整个条件独立关系的拓扑意义出发，然后推导出贝叶斯网络的数值意义。拓扑语义说的是给定每个节点的父节点，则该节点条件独立于所有它的非后裔(non-descendants)节点。举例来说，Figure 14.2的警报例子中，给定alarm后，JohnCalls独立于Burglary,Eqrthquake和MaryCalls。如图Figuree 14.4(a)中所示。从条件独立断言(assertions)和网络参数$\theta(x_i|parents(X_i))$就是条件概率$P(x_i|parents(X_i))$的解释中，联合概率可以计算出来。在这种情况下，数值意义和拓扑语义是相同的。
另一个拓扑意义的重要属性是：给定某个节点的马尔科夫毯(Markov blanket)，即节点的父节点，子节点，子节点的父节点，这个节点条件独立于所有其他的节点。如图Figure 14.4(b)所示。

条件分布的高效表示(Efficient representation of conditional distributions)

即使每个节点有$k$个父节点，一个节点的CPT还需要$O(2^k )$，最坏的情况下父节点和子节点是任意连接的。一般情况下，这种关系可以用符合一些标准模式(standard pattern)的规范分布(canonical distribution)表示，这样子就可以仅仅提供分布的一些参数就能生成整个CPT。
最简单的例子是确定性节点(deterministic node)。一个确定性节点的值被它的父节点的值精确确定。这个确定性关系可以是逻辑关系：父节点是加拿大，美国和墨西哥，子节点是北美洲，它们之间的关系是子节点是父节点所在的洲。这个关系也可以是数值型的，一条河的流量是流入它的流量减去流出它的流量。
不确定关系通常称为噪音逻辑关系(noisy logical relationships)。一个例子是噪音或(noisy-OR)，它是逻辑或的推广。在命题逻辑中，当且仅当感冒(Cold)，流感(Flu)或者疟疾(Malaria)是真的时候，发烧(Fever)才是真的。噪音或模型允许不确定性，即每一个父节点都有可能让子节点为真，可能父节点和子节点之间的关系被抑制了(inhibited)，可能一个人感冒了，但是没有表现出发烧。这个模型做了两个假设。第一个，它假设所有的原因都被列了出来，有时候会加一个节点(leak node)包含所有的其他原因(miscellaneous causes)。第二个，抑制每一个父节点和子节点之间的原因是独立的，比如抑制疟疾产生发烧和抑制感冒产生发烧的原因是独立的。所以，当且仅当所有的父节点都是假的时候，发烧才一定不会发生。给出以下的假设：
$q_{cold} = P(\neg fever| cold,\neg flu, \neg malaria) = 0.6$
$q_{flu} = P(\neg fever|\neg cold, flu, \neg malaria) = 0.2$
$q_{malaria} = P(\neg fever|\neg cold,\neg flu, malaria) = 0.1$
根据这些信息，以及噪音或的假设，整个CPT可以被创建。一般的规则是：
$P(x_i|parents(X_i)) = 1 - \prod_{j:X_j=ture} q_j.$
最后生成如下的表：

Cold	Flu	Malaria	P(Fever)	P($\neg$Fever)
F	F	F	$0.0$	$1.0$
F	F	T	$0.9$	$0.1$
F	T	F	$0.8$	$0.2$
F	T	T	$0.98$	$0.1\times 0.2=0.02$
T	F	F	$0.4$	$0.6$
T	F	T	$0.94$	$0.6\times 0.1 = 0.06 $
T	T	F	$0.88$	$0.5\times 0.2 = 0.12 $
T	T	T	$0.988$	$0.6\times 0.2\times 0.1 = 0.012$

对于这个表，感觉自己一直有点转不过来圈。就是有症状不一定发烧，也可能不发烧，没有症状一定不发烧。什么时候不发烧呢，只有某个症状表现出来不发烧，如果多个症状的话，直接把有症状表现但不发烧的概率相乘。
一般情况下，噪声逻辑模型中，有$k$个父节点的变量可以用$O(k)$个参数表示而不是$O(2^k )$去表示整个CPT。这让访问(assessment)和学习(learning)更容易了。

连续性随机变量的贝叶斯网络(Bayesian nets with continuous variables)

常用方法

现实中很多问题都是连续型的随机变量，它们有无数可能的取值，所以显式的指定每一个条件概率行不通。常用的总共有三种方法，第一个可能的方法是离散(discretization)连续型随机变量，将随机变量的可能取值划分成固定的区间。比如，温度可以分成，小于$0$度的，$0$度到$100$度之间的，大于$100$度的。离散有时候是可行的，但是通常会造成精度的缺失和非常大的CPT。第二个方法也是最常用的方法是通过指定标准概率密度函数的参数，比如指定高斯分布的均值和方差。第三种方法是非参数化(nonparametric)表示，用隐式的距离去定义条件分布。

示例

一个同时拥有离散型和随机性变量的网络被称为混合贝叶斯网络(hybrid Bayesian network)。为了创建这样一个网络，我们需要两种新的分布。一种是给定离散或者连续的父节点，子节点是连续型随机变量的条件概率，另一种是给定连续的父节点，子节点是离散型随机变量的条件概率。

连续型子节点

考虑Figure 14.5的例子，一个顾客买了一些水果，买水果的量取决取水果的价格(Cost)，水果的价格取决于收成(Harvest)和政府是否有补助(Subsidy)。其中，Cost是连续型随机变量，他有连续的父节点Harvest和离散的父节点Subsidy，Buys是离散的，有一个连续型的父节点Cost。
对于变量Cost，我们需要指定条件概率$P(Cost|Subsidy,Harvest)$。离散的父节点通过枚举(enumeration)来表示，指定$P(Cost|subsidy,Harvest)$和$P(Cost|\neg subsidy,Harvest)$。为了表示Harvest，可以指定一个分布来表示变量Cost的值$c$取决于连续性随机变量Harvest的值$h$。换句话说，将$c$看做一个$h$的函数，然后给出这个函数的参数即可，最常用的是线性高斯分布。比如这里，我们可以用两个不同参数的高斯分布来表示有补贴和没补贴时Harvest对Cost的影响：
$$P(c|h, subsidy) = N(a_th+b_t,\sigma_t^2 )© = \frac{1}{\sigma_t \sqrt{2\pi} }e^{- \frac{1}{2}(\frac{c-(a_th+b_t)}{\sigma_t})^2 } $$
$$P(c|h,\neg subsidy) = N(a_fh+b_t,\sigma_f^2 )© = \frac{1}{\sigma_f \sqrt{2\pi} }e^{- \frac{1}{2}(\frac{c-(a_fh+b_f)}{\sigma_f})^2 } $$
所以，只需要给出$a_t,b_t,\sigma_t,a_f,b_f,\sigma_f$这几个参数就行了，Figure 14.6(a)和(b)就是一个示例图。注意到坡度(slope)是负的，因为随着供应的增加，cost在下降，当然，这个线性模型只有在harvest在很小的一个区间内才成立，而且cost有可能为负。假设有补贴和没补贴的两种可能性相等，是$0.5$，那么就有了Figure 14.6©的图$P(c|h)$。

连续型父节点

当离散型随机变量有连续型父节点时，如Figure 14.5中的Buys节点。我们有一个合理的假设是：当cost高的时候，不买，cost底的时候，买，在中间区域买不买是一个变化很平滑的概率。我们可以把条件分布当成一个软阈值函数(soft-threshold)，一种方式是用标准正态分布的积分(intergral)。
$$\Phi(x) = \int_{-\infty}^{x} N(0,1)(x)dx$$
给定Cost买的概率可能是:
$$P(buys|Cost = c) = \Phi((-x+\nu)/ \sigma))$$
其中cost的阈值在$\nu$附近，阈值的区域和正比于$\sigma$，当价格升高的时候，买的概率会下降。这个probit distribution模型如Figure 14.7(a)所示。
另一个可选择的模型是logit distribution，使用logistic function $1/(1+e^{-x} )$来生成一个软阈值：
$$P(buys|Cost = c) = \frac{1}{1+exp(-2\frac{-c+u}{\sigma})}.$$
如Figure 14.7(b)所示，这两个分布很像，但是logit有更长的尾巴。probit更符合实际情况，但是logit数学上更好算。它们都可以通过对父节点进行线性组合推广到多个连续性父节点的情况。

贝叶斯网络的精确推理(Exact inference in bayesian networks)

概率推理系统的基本任务就是给出一些观察到的事件，即给证据变量(evidence variable)赋值，然后计算一系列查询变量(query variable)的后验概率。我们用$X$表示查询变量，用$\mathbf{E}$表示证据变量$E_1,\cdots,E_m$的集合，$\mathbf{e}$是一个特定的观测事件，$\mathbf{Y}$表示既不是证据变量，也不是查询变量的变量$Y_1,\cdots,Y_l$的集合（隐变量,hidden variables)。变量的所有集合是$\mathbf{X}={X}\cup \mathbf{E}\cup \mathbf{Y}$。一个典型的查询是求后验概率$P(X|\mathbf{e})$。
在这一节中主要讨论的是计算后验概率的精确算法以及这些算法的复杂度。事实上，在一般情况下精确推理的复杂度都是很高的，为了降低复杂度，就只能进行估计推理(approximate inference)了，这个会在下一节中介绍到。

枚举实现精确推理(Inference by enumeration)

任何条件概率都可以用联合概率分布的项相加得到，即：
$$P(X|\mathbf{e}) = \alpha P(X,\mathbf{e}) = \alpha \sum_{\mathbf{y}}P(X,\mathbf{e},\mathbf{y})$$
贝叶斯网络给出了所有的联合概率分布，任何项$P(x,\mathbf{e},\mathbf{y})$都可以用贝叶斯网络中的条件概率的乘积表示出来。比如警报例子中的查询$P(Burglary|JohnCalls=true,MaryCalls=true)$。隐变量是Earthquake和Alarm，我们可以算出：
$$P(B|j,m) = \alpha P(B,j,m) = \alpha \sum_{e}\sum_{a}P(B,j,m,e,a).$$
贝叶斯网络已经给出了所有CPT项的表达式，比如当Burglary = true时：
$$P(b|j,m) = \alpha \sum_e\sum_aP(b,j,m,e,a) = \alpha \sum_e\sum_aP(b)P(e)P(a|b,e)P(j|a)P(m|a).$$
为了计算这个表达式，我们得计算一个四项的加法，分别是e为true和false,a为true和false对应的$P(b,j,m)$的值，每一项都是五个数的乘法。最坏的情况下，所有的变量都用到了，那么拥有$n$个布尔变量的贝叶斯网络的时间复杂度是$O(n2^n )$。我们可以做一些简化，将一些重复的计算保存下来，比如将上面的式子变成：
$$P(b|j,m) = \alpha \sum_e\sum_aP(b,j,m,e,a) = \alpha P(b) \sum_eP(e)\sum_aP(a|b,e)P(j|a)P(m|a).$$
这样子可以按照顺序进行计算，具体的计算过程如Figure 14.8所示。这种算法叫做ENUMERATION-ASK，它的空间复杂度是线性的，但是它的事件复杂度是$O(2^n )$比$O(n2^n )$要好，却仍然是实际上不可行的。（这里我理解的是$O(2^n )$而不是$O(n2^n )$的原因是，总共有$n$个布尔变量，所以总共有$2^n $个可能的取值，每次算一个，存一个，而原来的是算完之后不存。）
事实上，Figure 14.8中的计算过程还有很多重复计算，比如$P(j|a)P(m|a)$和$P(j|\neg a)P(m|\neg a)$这两项被计算了两次。我原来在想这里是不是和上面一段说的冲突了，事实上是没有的，这$2^n $个值，其中可能会有$P(b,j,m,e,a)$和$P(b,j,m,e,\neg a)$，这两个概率中都用到了$P(j|a)P(m|a)$，但是这里就会计算两次，事实上有很多值都会被重复计算很多次。下面就介绍一个避免这种运算的方法。

消元法(The variable elimination algorithm)

上面问题的解决思路就是保存已经计算过的值，实际上这是一种动态规划。还有很多其他方法可以解决这个问题，这里介绍了最简单的消元算法。消元法对表达式进行从右至左的计算，而枚举法是自底向上的。所有的中间值被报存起来，最对和每个变量有关的表达式进行求和。例如对于下列表达式：
$$P(B|j,m) = \alpha \underbrace{P(B)}_{f_1(B)} \sum_e\underbrace{P(e)}_{f_2(E)} \sum_a\underbrace{P(a|B,e)}_{f_3(A,B,E)} \underbrace{P(j|a)}_{f_4(A)} \underbrace{P(m|a)}_{f_5(A)}.$$
表达式的每一部分都是一个新的因子，每一个因子都是由它的参数变量(argument variables)决定的矩阵，参数变量指定的取值是没有固定的变量。比如因子$f_4(A)$和$f_5(A)$对应$P(j|a)$和$P(m|a)$的表达式只取决于$A$的值因为$J$和$M$在这个查询中都是固定的。它们都是两个元素的向量：
$$f_4(A) = \begin{pmatrix}P(j|a)\\P(j|\neg a)\end{pmatrix} = \begin{pmatrix}0.90\\0.05\end{pmatrix}$$
$$f_5(A) = \begin{pmatrix}P(m|a)\\P(m|\neg a)\end{pmatrix} = \begin{pmatrix}0.70\\0.01\end{pmatrix}$$
$f_3(A,B,E)$是一个$2\times 2\times 2$的矩阵。用因子表达的话，查询的表达式变成了：
$$P(B|j,m) = \alpha f_1(B)\times \sum_ef_2(E)\times \sum_af_3(A,B,E)\times f_4(A)\times f_5(A)$$
其中$\times$不是普通的矩阵乘法，而是对应元素相乘(pointwise product)。整个表达式的计算过程可以看成从右到左变量相加的过程，将现有的因子消去产生新的因子，最后只剩下一个因子的过程。具体的步骤如下：
首先先利用$f_3,f_4,f_5$把变量$A$消掉，产生一个新的$2\times 2$的只含有变量$B$和$E$的新因子$f_6(B,E)$：
\begin{align*}
f_6(B,E) &= \sum_af_3(A,B,E)\times f_4(A) \times f_5(A)\\
&= (f_3(a,B,E)\times f_4(a) \times f_5(a)) + (f_3(\neg a,B,E)\times f_4(\neg a)\times f_5(\neg a)
\end{align*}
这样目标变成了：
$$P(B|j,m) = \alpha f_1(B)\times \sum_ef_2(E)\times \sum_af_6(B,E)$$
利用$f_2,f_6$消去$E$：
\begin{align*}
f_7(B) &= \sum_ef_2(E)\times \sum_af_6(B,E)\\
& = f_2(e)\times f_6(B,e) + f_2(\neg e)\times f_6(B,\neg e)
\end{align*}
将表达式化成：
$$P(B|j,m) = \alpha f_1(B)\times f_7(B)$$
显然，根据这个表达式就可以计算出我们想要的结果了。上面的过程可以总结成两步，第一步是point-wise的因子乘法，第二步是利用因子的乘法进行消元。

因子运算(Operations on factors)

两个因子$f_1$和$f_2$进行point-wise乘法运算产生新的因子(factor)$f$的变量是$f_1$和$f_2$变量的并，新的因子中的元素的值是$f_1$和$f_2$中对应项的积。假设两个因子有公共变量$Y_1,\cdots,Y_k$，那么就有：
$$f(X_1,\cdots,X_j,Y_1,\cdots,Y_k,Z_1,\cdots,Z_l)=f_1(X_1,\cdots,X_j,Y_1,\cdots,Y_k)f_2(Y_1,\cdots,Y_k,Z_1,\cdots,Z_l).$$
如果所有的变量都是二值化的，那么$f_1$和$f_2$各有$2^{j+l} $和$2^{l+k} $项，$f$有$2^{j+l+k} $项。比如，$f_1(A,B),f_2(B,C)$，那么point-wise乘法产生的$f_3(A,B,C)=f_1\times f_2$有$8$项，如Figure 14.10所示。

根据图中给出的值，消去$f_3(A,B,C)$中的$A$：
\begin{align*}
f(B,C) &= \sum_af_3(A,B,C)\\
&= f_3(a,B,C) + f_3(\neg a,B,C)\\
&= \begin{pmatrix} 0.06&0.24\\0.42&0.28\end{pmatrix} + \begin{pmatrix}0.18&0.72\\0.06&0.04\end{pmatrix}\\
&= \begin{pmatrix}0.24&0.96\\048&0.32\end{pmatrix}
\end{align*}
产生新的因子用的是pointwise乘法，消元用的是累乘。给定pointwise乘法和消元函数，消元算法就变得很简单，一个消元算法如Figure 14.11所示。

变量顺序和变量相关性(Variable ordering and variable relevance)

Figure 14.11中的算法包含一个没有给出具体实现的排序函数Order()对要消去的变量进行排序，每一种排序选择都会产生一组有效的算法，但是不同的消元顺序会产生不同的中间因子。一般情况下，消元法的时间和空间复杂度是由算法产生的最大因子决定的，这个最大因子是由消元的顺序和贝叶斯网络的结构决定的，选取最优的消元顺序是很困难的，但是有一些小的技巧：总是消去让新产生的因子最小的变量。
另一个属性是：每一个不是查询变量或者证据变量的祖先变量都和这次查询无关，在实现消元算法的时候可以把这些变量都去掉。（具体的示例可以看第十四章，在$528$页）。

精确推理的复杂度(The complexity of exact inference)

贝叶斯网络的精确推理跟网络的结构有很大的关系。
Figure 14.2中警报贝叶斯网络中的复杂度是线性的。该网络中任意两个节点只有一条路径，这种网络称为单连接的(singly-connected)或者多树(polytrees)，这种结构有一个很好的属性就是：多树结构中精确推理的时间，空间复杂度对于网络大小来说都是线性关系，这里网络大小指的是CPT项的个数。如果每一个节点的父节点都是一个有界的常数，那么复杂度和节点数之间也是线性关系。
对于多连接(multiply connected)的网络，如Figure 14.12(a)所示，最坏情况下，即使每一个节点的父节点个数都是有界常数，消元法的时间和空间复杂度也都是指数级别的。因为贝叶斯网络的推理也是NP难问题。

聚类算法(clustering algorithms)

用消元法来计算单个的后验概率是简单而高效的，但是如果要计算网络中所有变量的后验概率是很低效的。例如：在单连接的网络中，每一个查询都是$O(n)$，总共有$O(n)$个查询，所以总共的代价是$O(n^2 )$。使用聚类算法(clustering algorithms)，代价可以降到$O(n)$，因此贝叶斯网络中的聚类算法已经被广泛商用。（这里不明白为什么？）。
聚类算法的基本思想是将网络中的一些节点连接成聚点(cluster nodes)，最后形成一个多树(polytree)结构。例如Figure 14.12(a)中的多连接网络可以转换成Figure 14.12(b)所示的多树，Sprinkler和Rain节点形成了SPrinkler+Rain聚点，这两个布尔变量被一个大节点(meganode)取代，这个大节点有四个可能的取值：$tt,tf,ft,ff$。一旦一个多树形式的网络生成了以后，就需要特殊的推理算法进行推理了，因为普通的推理算法不能处理共享变量的大节点，有了这样一个特殊的算法，后验概率的时间复杂度就是线性于聚类网络的大小。但是，NP问题并没有消失，如果消元需要指数级别的时间和空间复杂度，聚类网络中的CPT也是指数级别大小。

贝叶斯网络的估计推理(Approximate inference in bayesian networks)

因为多连接网络中的推理是不可行的，所以用估计推理取代精确推理是很有用的。这一节会介绍随机采样算法，也叫蒙特卡洛算法(Monte Carlo)，它的精确度取决于生成的样本数量。我们的目的是采样用于计算后验概率。这里给出了两类算法，直接采样(direct sampling)和马尔科夫链采样(Markov chain sampling)。变分法(variational methods)和循环传播(loopy propagation)将会在本章的最后进行介绍。

直接采样(Direct sampling methods)

任何采样算法都是通过一个已知的先验概率分布生成样本。比如一个公平的硬币，服从一个先验分布$P(coin) = <0.5,0.5 > $，从这个分布中采样就像抛硬币。
一个最简单的从贝叶斯网络中进行随机采样的方法就是：从没有证据和它相关的网络中生成事件，即按照拓扑顺序对每一个变量进行采样。如Figure 14.13所示的算法，每一个变量的采样都取决于前之前已经采样过了的父节点变量的值。按照Figure 14.13中的算法对Figure 14.12(a)中的网络进行采样，假设一个采样顺序是[Cloudy,Sprinkler,Rain,WetGrass]：

从$P(Cloudy)=<0.5,0.5>$中采样，采样值是true；
从$P(Sprinkler|Cloudy=true) = <0.1,0.9>$中采样，采样值是false；
从$P(Rain|Cloudy=true)=<0.8,0.2>$中采样，采样值是true；
从$P(WetGrass|Sprinkler=false,Rain=true)=<0.9,0.1>$中采样，采样值是true；

这个例子中，PRIOR-SAMPLE算法返回事件[true,false,true,true]。可以看出来，PRIOR-SAMPLE算法根据贝叶斯网络指定的先验联合分布生成样本。假设$S_{PS}(x_1,\cdot,x_n)$是PRIOR-SAMPLE算法生成的一个样本事件，从采样过程中我们可以得出：
$$S_{PS}(x_1,\cdots,x_n) = \prod_{i=1}^n P(x_i|parents(X_i))$$
即每一步采样都只取决于父节点的值。这个式子和贝叶斯网络的联合概率分布是一样的，所以，我们可以得到：
$$S_{PS} = P(x_1,\cdots,x_n).$$
通过采样让这个联合分布的求解很简单。
figure 14.13
事实上在任何采样算法中，结果都是通过对产生的样本进行计数得到的。假设生成了$N$个样本，$N_{PS}(x_1,\cdots,x_n)$是样本集中的一个具体事件$(x_1,\cdots,x_n)$发生的次数。我们希望这个值比上样本总数取极限和采样概率$S_{PS}$是一样的，即：
$$ lim_{N\rightarrow \infty}\frac{N_{PS}(x_1,\cdots,x_n)}{N} = S_{PS}(x_1,\cdots,x_n) = P(x_1,\cdots,x_n).$$
例如之前利用PRIOR-SAMPLE算法产生的事件[true,false,true,true]，这个事件的采样概率是：
$$S_{PS}(true,false,true,true) = 0.5 \times 0.9 \times 0.8 \times 0.9 = 0.324.$$
即当$N$取极限时，我们希望有$32.4%$的样本都是这个事件。(这里为什么要用采样进行计算呢，我的想法是因为实际情况中，采样概率$S_{PS}$是很难计算的，就通过不断的采样，计算出某个样本出现的概率。)
我们用$\approx$表示估计概率(estimated probability)在样本数量$N$取极限时和真实概率一样的估计，这叫一致(consistent)估计。比如，对于任意的含有隐变量的事件(partially spefified event)，$x_1,\cdots,x_m,m\le n$，会产生一个一致估计：
$$P(x_1,\cdots,x_m)\approx N_{PS}(x_1,\cdots,x_m)/N.$$
这个事件的概率可以看成所有满足观测变量条件的样本事件（隐变量所有值都可以取）比上所有样本事件的比值。比如在Spinkler网络中，生成$1000$个样本，其中有$511$个样本的Rain=true，那么rain的估计概率就是$\hat{P}(Rain=true) = 0.511.$

贝叶斯网络的拒绝采样(Rejection sampling in Bayesian networks)

算法

拒绝采样(rejection sampling)利用容易采样的分布来生成难采样分布的样本，计算后验概率$P(X|\mathbf{e})$，算法流程如Figure 14.14所示，首先根据贝叶斯网络的先验分布生成样本，接下来拒绝(reject)那些和证据变量不匹配的结果，最后在剩下的样本中统计每个$X=x$出现的概率，估计$\hat{P}(X|\mathbf{e}).$
用$\hat{P}(X|\mathbf{e})$表示估计概率分布，利用拒绝采样算法的定义计算：
$$\hat{P}(X|\mathbf{e}) = \alpha N_{PS}(X,\mathbf{e}) = \frac{N_{PS}(X,\mathbf{e})}{N_{PS}(\mathbf{e})}.$$
而根据$P(x_1,\cdots,x_m)\approx N_{PS}(x_1,\cdots,x_m)/N$，就有：
$$\hat{P}(X|\mathbf{e}) = \alpha N_{PS}(X,\mathbf{e}) = \frac{N_{PS}(X,\mathbf{e})}{N_{PS}(\mathbf{e})} = \frac {P(X,\mathbf{e})}{P(\mathbf{e})} = P(X|\mathbf{e}).$$
所以，拒绝采样产生了真实概率的一个一致估计(consistent estimate)，但是这个一致估计和无偏估计还不一样。

示例

举一个例子来说明，假设我们要估计概率$P(Rain|Sprinkler=true)$，生成了$100$个样本，其中$73$个是$Sprinkler=false$，$27$是$Sprinkler=true$，这$27$个中有$8$个$Rain=true$，有$19$个$Rain=false$，因此：
$$P(Rain|Sprinkler=true)\approx NORMALIZE \lt\lt 8,19>> = <0.296,0.704>.$$
正确答案是$<0.3,0.7>$，可以看出来，估计值和真实值差的不多。生成的样本越多，估计值就会和正确值越接近，概率的估计误差和$1/\sqrt{n}$成比例，$n$是用来估计概率的样本数量。

不足

拒绝采样最大的问题是它拒绝了很多样本，随着证据变量的增加，和证据$\mathbf{e}$一致的样本指数速度减少，所以这个方法对于复杂的问题是不可行的。拒绝假设和现实生活中条件概率是很像的，比如估计观测到晚上天空是红的，第二天下雨的概率$P(Rain|RedSkyAtNight=ture)$，这个条件概率的估计就是根据日常生活的观察实现的。但是如果天空很少是红的，就需要很长时间才能估计它的值，这就是拒绝假设的缺点。

可能性加权(Likelihood weighting)

算法

可能性加权(Likelihood weighting)只产生和证据$\mathbf{e}$一致的事件，因此避免了拒绝采样的低效。它是统计学中重要性采样的一个例子，专门为贝叶斯推理设计的。
如Figure 14.15所示，加权似然固定证据变量$\mathbf{E}$的值，只对非证据变量进行采样，这就保证了每一个事件都是和证据一致的。但是，不是所有的事件权重都是一样的。给定每一个证据变量的父节点，它的可能性(likelihood)是证据变量的条件概率的乘积，每一个事件都根据证据的可能性进行加权。

示例

对于Figure 14.12(a)中的例子，计算后验概率$P(Rain|Cloudy=true,WetGrass=true)$，采样顺序是Cloudy,Sprinkler,Rain,WetGrass。过程如下，首先，权重$w$设为$1$，一个事件生成过程如下：

Cloudy是一个证据变量，它的值是true,因此，令：
$$w\leftarrow w\times P(cloudy=true) = 0.5.$$
Sprinkler是隐变量，所以从$P(Sprinkler|Cloudy=true)=<0.1,0.9>$中采样，假设采样结果是false；
Rain是隐变量，从$P(Rain|Cloudy=true)=<0.8,0.2>$中采样，假设采样结果是true；
WetGrass是证据变量，值是true,令：
$$w\leftarrow w\times P(WetGrass=true|Sprinkler=false,Rain=true) = 0.45.$$

所以WEIGHTED-SAMPLE算法生成事件[true,false,true,true]，相应的权重是$0.45$。

原理

用$S_{WS}$表示WEIGHTED-SAMPLE算法中事件的采样概率，证据变量$\mathbf{E}$的取值$\mathbf{e}$是固定的，用$\mathbf{Z}$表示非证据变量，包括隐变量$\mathbf{Y}$和查询变量$\mathbf{X}$。给定变量$\mathbf{Z}$的父节点，算法对变量$\mathbf{Z}$进行采样：
$$S_{WS}(\mathbf{z},\mathbf{e}) = \prod_{i=1}^l P(z_i|parents(Z_i)).$$
其中$Parents(Z_i)$可能同时包含证据变量和非证据变量。
和先验分布$P(\mathbf{z})$不同的是，每一个变量$Z_i$的取值会受到$Z_i$的祖先(ancestor)变量的影响。比如，对Sprinkler进行采样的时候，算法会受到它的父节点中的证据变量Cloudy=true的影响，而先验分布不会。另一方面，$S_{WS}$比后验分布$P(\mathbf{z}|\mathbf{e})$受证据的影响更小，因为对$Z_i$的采样忽略了$Z_i$的非祖先(non-ancestor)变量中的证据。比如，对Sprinkler和Rain进行采样的时候，算法忽略了子节点中的证据变量WetGrass=true，事实上这个证据已经排除了(rule out)Sprinkler=false和Rain=false的情况，但是WEIGHTED-SAMPLE还会产生很多这样的样本事件。
理想情况下，我们想要一个采样分布和真实的后验概率$P(\mathbf{z}|\mathbf{e})$相等，不幸的是不存在这样的多项式时间的算法。如果有这样的算法的话，我们可以用多项式数量的样本以任意精度逼近想要求的概率值。
可能性权重$w$弥补了实际的分布和我们想要的分布之间的差距。一个由$\mathbf{z}$和$\mathbf{e}$组成的样本$\mathbf{x}$的权重是给定了父节点的证据变量的可能性乘积：
$$w(\mathbf{z},\mathbf{e}) = \prod_{i=1}^m P(e_i|parents(E_i)).$$
将上面的两个式子乘起来，可以得到一个样本的加权概率(weighted probability)是：
$$S_{WS}(\mathbf{z},\mathbf{e})w(\mathbf{z},\mathbf{e}) = \prod_{i=1}^l P(z_i|parents(Z_i))\prod_{i=1}^m P(e_i|parents(E_i)) = P(\mathbf{z},\mathbf{e}).$$
可能性加权估计是一致估计。对于任意的$x$，估计的后验概率按下式计算：
\begin{align*}
\hat{P}(x|\mathbf{e}) &= \alpha \sum_{\mathbf{y}} N_{WS}(x,\mathbf{y},\mathbf{e})w(x,\mathbf{y},\mathbf{e})\\
&\approx \alpha’\sum_{\mathbf{y}}S_{WS}(x,\mathbf{y},\mathbf{e})w(x,\mathbf{y},\mathbf{e})\\
&=\alpha’\sum_{\mathbf{y}}P(x,\mathbf{y},\mathbf{e})\\
&=\alpha’\sum_{\mathbf{y}}P(x,\mathbf{y},\mathbf{e})\\
&=P(x|\mathbf{e})
\end{align*}
算法中真实实现的是第一行，即统计出用WEIGHTED-SAMPLE产生的样本$(x,\mathbf{y},\mathbf{e})$数量$N_{WS}$，以及对应的权重$w(x,\mathbf{y},\mathbf{e})$，后面的都是理论推导，当$N$取极限的时候$lim_{N\rightarrow \infty}\frac{N_{WS}(x_1,\cdots,x_n)}{N} = S_{WS}(x_1,\cdots,x_n)$，后面的都是为了证明算法是一致估计。

不足

可能性加权算法使用了所有生成的样本，它比拒绝假设算法更高效。然而，随着证据变量的增加，算法性能会退化(degradation)，这是因为很多样本的权重都会很小，因此加权估计可能会受一小部分权重很大的样本的影响(dominated)。如果证据变量在非证据变量的后边，这个问题会加剧，因为它们的父节点或者祖先节点没有证据变量来指导样本的生成。这就意味着生成的样本和证据变量支撑的真实情况可能差距很大(bear little resemblance)。

马尔科夫链仿真推理(Inference by Markov chain simulation)

马尔科夫链蒙特卡洛(Markov chain Monte Carlo,MCMC)算法和拒绝采样以及可能性加权很不一样。那两个方法每次都从头开始生成样本，而MCMC算法在之前的样本上做一些随机的变化。可以将MCMC算法看成指定了每一个变量值的特殊当前状态(current state)，通过对当前状态(current state)做任意的改变生成下一个状态(next state)。这一节要介绍的一种MCMC算法是吉布森采样(Gibbs sampling)。

贝叶斯网络中的吉布森采样(Gibbs sampling in Bayesian networks)

算法

贝叶斯网络中的吉布森采样从任意一个状态开始，其中证据变量的取值固定为观测值，通过随机选取非证据变量$X_i$的值生成下一个状态。变量$X_i$的采样取决于变量$X_i$的马尔科夫毯的当前值。算法在状态空间（所有非证据变量的全部可能取值空间）中随机采样，每次采样都保持证据变量不变，一次改变一个非证据变量的值。完整的算法如Figure 14.16所示。

示例

Figure 14.12(a)中的查询(query)$P(Rain|Sprinkler=true, WetGrass=true)$，证据变量Spinkler和WetGrass取它们的观测值不变，非证据变量Cloudy和Rain随机初始化，假设取的是true和false。那么初始状态就是[true,true,false,true]，接下来对非证据变量进行重复的随机采样。
比如第一次对Cloudy采样（也可以对Rain采样），给定它的马尔科夫毯变量，然后从$P(Cloudy|Sprinkler=true,Rain=false)$中进行采样，假设采样结果是false，新的状态就是[false,true,false,true]。接下来随机可以对Rain采样（也可以对Cloudy采样），给定Rain的马尔科夫毯变量的取值，从$P(Rain|Cloudy=false,Sprinkler=true,WetGrass=true)$中进行采样，假设采样值是true,那么新的状态是[true,true,false,false]。接下来可以一直进行采样。。最终利用生成的样本计算出相应的概率。

为什么吉布森采样有用(Why Gibbs sampling works)

接下来给出为什么吉布森采样计算后验概率是一致估计。基本的解释是直截了当的：采样过程建立了一个动态平衡，每个状态花费的时间长期来说和它的后验概率是成比例的。
具体的，不想看了。。。就随缘吧

参考文献

1.《Aritifical Intergence》

PRML chapter 8 Graphical Models

发表于 2019-01-06 | 更新于 2019-12-17 | 分类于机器学习

$\newcommand{\mmm}{\mathbf}$
概率图模型
概率论在现代模式识别中有很重要的地位。第一章中介绍了概率论可以被表示成两个简单的加法和乘法公式。事实上在这本书中讨论的所有概率推理和学习的计算（无论有多复杂）都可以看成这两个公式的重复应用。我们可以只用代数计算(algebraic manipulation)来形式化并解决复杂的概率问题。但是，使用概率分布(probability distributions)的图表示(diagrammatic representations)，即概率图模型(graphical models)会更有优势。概率图模型有以下几个有用的属性：

概率图模型提供了一个简单的方式可视化(visualize)概率模型的结构，并且能够用来设计和产生新的模型。
通过观察概率图模型，可以看到模型的一些属性，包括条件独立性(conditional independence)等等。
在复杂模型上进行的需要推理和学习的复杂计算，可以被表示为图计算，底层的数据表达式隐式的被执行。

一个图由节点(nodes)，有时也叫顶点(vertices)，连接顶点的连接(links)，也叫边(edges)。在一个概率图模型中，每一个节点代表一个随机变量，或者一组随机变量，边代表着变量之间的概率关系(probabilistic relationships)。所有随机变量的联合分布可以被分解成一系列部分随机变量的乘积。
本章从有向图(directed graphical models)中的贝叶斯网络(Beyesian networks)开始介绍，有向图中的边通过箭头表示方向。另一个主要的图模型是马尔科夫随机场(Markov random fields)，它是一个无向图模型(undirected graphical models)，没有明显的方向性。有向图用来描述随机变量之间的因果关系(causal relationships)，而无向图用来描述随机变量之间的一些软约束(soft constraints)。为了解决推理问题，将无向图和有向图转化成另一种因子图(factor graph)表示是很方便的。

贝叶斯网络(Bayesian Networks)

图的一个很强大的特点就是一个具体的图可以用来解释一类概率分布。给定随机变量$a,b,c$的联合概率分布$p(a,b,c)$，通过利用乘法公式，我们可以把它写成以下形式：
$$p(a,b,c) = p(c|a,b)p(b|a)p(a).$$
这个公式对于任意的联合分布都成立，我们用节点$a,b,c$表示随机变量，按照上式的右边找出每个节点对应的条件分布，在图中添加一个有向箭头从依赖变量指向该变量。如Figure 8.1所示,$a$到$b$的边表示$a$是$b$的父节点。上式中左边是$a,b,c$是对称的，但是右边不是，事实上，在做分解的时候，一个隐式的顺序$a,b,c$已经被确定了，当然也可以选其他顺序，这样会得到一个新的分解和一个新的图。

如果把三个变量可以扩展到$K$个变量，则对应的联合概率为$p(x_1,\cdots,x_k)$，写成如下形式：
$$p(x_1,\cdots,x_K) = p(x_K|x_{K-1},\cdots,x_1)\cdots p(x_2|x_1)p(x_1).$$
这个式子也叫链式法则，微积分中也有链式法则，这个是概率论中的链式法则。给定$K$值，我们也能生成一个含有$K$个节点的有向图，每一个节点都对应一个条件分布，每一个节点都和比它序号小的节点全部直接相连，所以这个图也叫全连接图，因为任意两个节点都直接相连，但是只有一条有向边由小号节点指向大号节点，所以没有环。

目前为止，所有的操作都是在完全的联合概率分布，相应的分解以及全连接网络上进行的，它们可以应用到任何分布。但是图中也可能有缺失的边，如Figure 8.2所示，它不是一个全连接的图。我们可以直接根据这个图将联合分布表示为很多条件分布的乘积。每一个条件分布的取值只跟图中对应的父节点。比如，$x_5$只取决于$x_1$和$x_3$，$7$个变量的联合概率分布可以写成：
$$p(x_1)p(x_2)p(x_3)p(x_4|x_1,x_2,x_3)p(x_5|x_1,x_3)p(x_6|x_4)p(x_7|x_4,x_5)$$
从上面我们可以看出有向图和变量的条件概率之间的关系，图中定义的联合概率分布是图中所有节点给定其父节点的条件概率的乘积，即:
$$p(\mathbf{x}) = \prod_{k=1}^Kp(x_k|pa_k).$$
其中$pa_k$是$x_k$节点的父节点的集合，$\mathbf{x} = {x_1,\cdots,x_k}$，这个式子给出了一个有向图的联合概率具有因式分解属性。
贝叶斯网络中不能存在有向的圈，即不能存在闭路，所以这种图也叫有向无环图。另一种说法是如果图中的节点有顺序的话，不能存在大号节点到小号节点的有向边。

示例：多项式回归(Example: Polynomial regression)

figure 8.3
这里给出了一个用有向图描述概率分布的例子，贝叶斯多项式回归模型。模型中的随机变量是多项式系数向量$\mathbf{w}$以及观测值$\mathbf{t}=(t_1,\cdots,t_N)^{T$，此外，还有一些模型中确定的参数，它们不是随机变量，如输入数据$\mathbf{x}=(x_1,\cdots,x_N)}T$，噪音方差$\sigma^2$，还有$\mathbf{w}$上高斯分布精度的超参数$\alpha$。如果只关注随机变量，联合分布可以看成先验分布$p(\mathbf{w})$和$N$个条件分布$p(t_n|\mathbf{w}),n=1,\cdots,N$的乘积：
$$p(\mathbf{t},\mathbf{w}) = p(\mathbf{w})\prod_{n=1}^Np(t_n|\mathbf{w}).$$
这个模型可以用Figure 8.3表示。

为了方便表示，我们把$t_1,\cdots,t_N$用一个单独的节点，外面用一个盒子包着，叫做盘子(plate)，盘子上写上$N$代表有$N$个这样的节点，得到Figure 8.4中的图。如果把模型确定的参数写出来，我们可以得到下式：
$$p(\mathbf{t},\mathbf{w}|\mathbf{x},\alpha,\sigma^2) = p(\mathbf{w}|\alpha)\prod_{n=1}^{Np(t_n|\mathbf{w},x_n,\sigma}2).$$

如果在图中把模型参数和随机变量都表示出来，用空心圆圈代表随机变量，用实心圆点代表确定性参数(deterministic parameters)，用图形表示如Figure 8.5。

当用图模型去解决机器学习或者模型时，有时候会固定一些随机变量的值，比如在多项式拟合问题中训练集的变量${t_n}$，在图模型中，将对应节点加上阴影，表示观测变量(observed variables)。如Figure 8.6所示，变量${t_n}$是观测变量。$\mathbf{w}$没有被观测到，所以是一个隐变量(latent variable)或者是(hidden variable)。
利用观测到的${t_n}$的值，我们可以估计多项式系数$\mathbf{w}$，利用贝叶斯公式：
$$p(\mathbf{w}|\mathbf{T}) \propto p(\mathbf{w}) \prod_{n=1}^Np(t_n|\mathbf{w})$$
为了整洁(uncluttered)，模型的确定性参数被略去了。

一般来说，我们对于如$\mathbf{w}$之类的模型参数不感兴趣，因为我们的目标是用模型对新的输入进行预测。即在给定观测数据之后，我们给出一个新的输入$\hat{x}$，要找到对应的$\hat{t}$的概率分布，如Figure 8.7所示。给定确定性参数之后，图中所有随机变量的联合分布如下所示：
$$p(\hat{t},\mathbf{t},\mathbf{w}|\hat{x},\mathbf{x},\alpha,\sigma^2) = \left[\prod_{n=1}^{Np(t_n|x_n,\mathbf{w},\sigma}2)\right] p(\mathbf{w}|\alpha)p(\hat{t}|\hat{x},\mathbf{w},\sigma^2).$$
刚开始有一些不理解，但是实际上就是这样一个公式$p(a,b,c) = p(a)p(b|a)p(c|a)$，把$\hat{t}$和$\mathbf{t}$当成两个变量看就行了。
利用概率论的加法公式$p(X) = \sum\limits_Yp(X,Y)$，对模型的参数$\mathbf{w}$积分就得到了$\hat{t}$的预测分布：
\begin{align*}
p(\hat{t}|\hat{x},\mathbf{x},\mathbf{t},\alpha,\sigma^2) = \int p(\hat{t},\mathbf{w}|\hat{x},\mathbf{x},\mathbf{t},\alpha,\sigma^2)d\mathbf{w}
\propto \int p(\hat{t},\mathbf{t},\mathbf{w}|\hat{x},\mathbf{x},\alpha,\sigma^2)d\mathbf{w}
\end{align*}
其中随机变量$\mathbf{t}$被隐式的赋值为数据集中的观测值，即是一个$p(t)$是一个定值。这里刚开始有些不理解,实际上是当$p(b)$为定值的时候，$p(a|b) \propto p(ab)$。

生成模型(Generative models)

这里实际上介绍的是采样方法，叫祖先采样，实际上就是直接采样，AI的第十四章有讲很多采样，可以直接看那个。
很多时候我们需要从一个给定的分布中进行采样，十一章还会更详细的讲采样，这里要介绍一种采样分布叫祖先采样(ancestral sampling)，是一种和概率图模型相关的采样方法。给定$K$个变量的联合分布$p(x_1,\cdots,x_K)$对应的有向无环图，假设所有变量的父节点的序号都比它本身小。我们的目标是从联合分布中采样$\hat{x_1},\cdots,\hat{x_k}$。
首先从最小的序号根据$p(x_1)$开始采样，采样结果称为$\hat{x_1}$，接下来按顺序对第$n$个节点按照条件分布$p(x_n|pa_n)$进行采样，每个节点的父节点都取采样值，因为每个父节点都已经采完样了，所以这里不用担心。一直到第$K$个节点采样完成，就生成了一个样本。为了对某些边缘分布进行采样，对需要的节点进行采样，忽略其他节点即可，比如为了对边缘分布$p(x_2,x_4)$进行采样，从联合分布中进行采样，保留$\hat{x_2},\hat{x_4}$的值，其他的值不用管即可。
在概率图的实际应用中，通常小节点对应的是隐变量，大节点对应的图上的最终节点代表着一些观测变量。隐变量的目的是让观测变量的复杂概率分布可以表示成多个简单的条件概率分布的乘积。
我们可以把这样的模型解释为观测变量产生的过程，比如，一个模式识别任务中，每一个观测数据对应一张图片。隐变量解释为物体的位置和方向，给定一个观测图像，我们的目标是找到物体的一个后验分布，在后验分布中对所有可能的位置和方向进行积分，如Figure 8.8所示。
图模型通过观测数据的生成过程描述了一种因果关系过程，因为这个原因，这样的模型也叫做生成式模型(generative model)。相反，Figure 8.5中的模型不是生成式模型，因为多项式回归模型中的输入变量$x$没有概率分布，所以不能用来合成数据。通过引入一个合适的先验分布$p(x)$，我们可以把它变成一个生成式模型。
事实上，概率图模型中的隐变量不是必须要有显式的物理意义，它的引入只是为了方便从简单的条件概率生成复杂的联合分布。在任何一种情况下，应用到生成式模型的祖先采样模拟了观测数据的生成过程，因此产生了和观测数据分布相同（如果模型完美的表现了现实）的美好(fantasy)数据。实际应用中国，利用生成模型产生合成的观测数据，对于理解模型表达的概率分布很有帮助。

离散型随机变量(Discrete variables)

指数分布是很重要的一类分布，它们虽然很简单，但是可以形成更复杂的概率分布，概率图的框架对于表达这些概率分布是如何连接的很有用。
如果我们有向图中亲本和子节点对之间的关系选择为conjugate，会发现这些模型有很好的属性。这里主要探讨两种情况，父节点和子节点都是离散的以及父节点和子节点都对应高斯变量，因为这两种关系可以分层扩展(extended hierarchically)构建任何复杂的有向无环图。首先从离散变量开始：
有$K$个可能状态的单个离散变量$\mathbf{x}$的概率分布是：
$$p(\mathbf{x}|\nu) = \prod_{k=1}^k\nu_k{x_k}$$
由参数$\nu = (\nu_1,\cdots,\nu_K)^T$控制，由于有约束条件$\sum_k\nu_k=1$，为了定义分布有$K-1$个$\nu_k$的值需要指定。
假设有两个离散型随机变量$\mathbf{x}1,\mathbf{x}2$，每个变量都有$K$个可能的取值。用$\nu{kl}$表示同时观测到$x{1k}=1$和$x_{2l}=1$，其中$x_{1k}$表示$\mathbf{x}1$的第$k$个分量，$x{2l}$类似。联合分布可以写成：
$$p(\mathbf{x}1,\mathbf{x}2|\nu) = \prod{k=1}^K\prod_{l=1}K\nu{kl}^{x_{1k}x_{2l}}.$$
$\nu_{kl}$满足约束条件$\sum_k\sum_l\nu_{kl} =1$，被$K^{2-1$个参数控制，任意$M$个具有$K$个取值的随机变量的联合分布需要$K}M-1$个参数，随着随机变量$M$个数的增加，参数的个数以指数速度增加。
使用乘法公式，联合分布$p(\mathbf{x}_1,\mathbf{x}_2)$可以分解成$p(\mathbf{x}_2|\mathbf{x}_1)p(\mathbf{x}_1)$，对应的图如Figure 8.9(a)所示，边缘分布$p(\mathbf{x}_1)$的分布需要$K-1$个参数，$p(\mathbf{x}_2|\mathbf{x}_1)$对于$K$个可能的$\mathbf{x}_1$，每个都需要$K-1$个参数。所以，和联合分布一样，总共需要的参数为$K-1+K(K-1) = K^2-1$个。
假设$\mmm{x}_1$和$\mmm{x}_2$是独立的，如Figure 8.9(b)所示，每一个变量可以用分开的多峰分布(multinomial distribution)表示，所需的参数量为$2(K-1)$个。类似的，$M$个独立变量需要$M(K-1)$个参数，和变量个数之间是线性关系。从概率图的角度来看，通过在图中去掉边减少了参数的数量，同时代价是这个图只能代表有限类别的分布。
更普通的是，如果我们有$M$个离散型随机变量$\mmm{x}_1,\cdots,\mmm{x}_M$，我们可以用一个节点代表一个随机变量，建立一个有向图表示联合概率分布。每个节点处的条件概率由一组非负参数给出，并且需要满足归一化条件。如果图是全连接的，那么这个分布需要$K^M-1$个参数，如果图中没有连接，那么联合分布可以分解成边缘分布的乘积，需要的所有参数是$M(K-1)$个。拥有中间水平连接性的图比分解成单个边缘分布的乘积能解释更多的分布同时比普遍的联合概率分布需要更少的参数。如Figure 8.10中的节点链，边缘分布$p(\mmm{x}_1)$需要$K-1$个参数，其余的$M-1$个条件分布$p(\mmm{x}i|\mmm{x}{i-1}),i = 2,\cdots,M$，需要$K(K-1)$个参数，总共需要的参数是$K-1+(M-1)K(K-1)$个，是$K$的二次函数(quadratic)，随着链的长度$M$增加，参数个数线性增加。
另一个减少模型中独立参数个数的方法是共享参数。例如，Figure 8.10中的链，我们可以用共享的$K(K-1)$个参数去控制条件概率$p(\mmm{x}i|\mmm{x{i-1}),i=2,\cdots,M$，用$K-1$个变量去控制$\mmm{x}_1$的概率分布，总共需要$K^2-1$个参数需要被指定去定义联合概率分布。
通过引入每个参数对应的Dirichlet先验，我们可以把一个随机变量图转换成贝叶斯模型。从概率图的角度来看，每一个节点需要一个额外的父节点代表和这个离散节点相关的Dirichlet分布，如Figure 8.11所示。将控制条件分布$p(\mmm{x}i|\mmm{x}{i-1}),i=2,\cdots,M$的参数共享，得到如Figure 8.12所示的图。
另一种控制模型中离散变量参数指数速度增加的方法是用参数模型而不是条件概率表来表示条件分布。如Figure 8.13中，所有的节点都是一个二值变量，用参数$\nu_i$表示每一个父节点$\mmm{x}i$取值为$1$的概率$p(x_i=1)$，总共有$M$个父节点，所以总共需要$2^{M$个参数表示条件概率$p(y|x_1,\cdots,x_M)$的$2}M$可能取值，如$p(y=1)$。所以指定这个条件分布所需要的参数随着$M$指数级增长。我们可以通过使用logistic sigmoid函数作用于父节点的线性组合上，得到一个更简洁的条件概率分布：
$$p(y=1|x_1,\cdots,x_M) = \sigma\left(w_0+\sum{i=1}^{Mw_ix_i\right)=\sigma(\mmm{w}}T\mmm{x})$$
其中$\sigma(a) = \frac{1}{1+exp(-a)}$是logistic sigmoid，$\mmm{x}=(x_0,x_1,\cdots,x_M)^{T$是由$M$个父节点状态和一个$x_0=1$构成的$M+1$维向量，$\mmm{w}=(w_0,w_1,\cdots,w_M)}T$是$M+1$维参数项。与一般情况相比，这是一个更加严格的条件概率分布形式，但是它的参数个数随着$M$的增加线性增加。在这种情况下，类似于选择多元高斯分布的协方差矩阵的限制形式（如对角矩阵等）。

线性高斯模型(Linear-Gaussian models)

条件独立性(Conditional Independence)

马尔科夫随机场(Markov Random Fields)

概率图模型中的推理(Inference in Graphical Models)

参考文献(references)

ESL chapter 1 Introduction

发表于 2019-01-05 | 更新于 2019-12-17 | 分类于机器学习

引言

这本书主要介绍的是统计学习。一些典型的学习问题如下：

基于一个心脏病患者的饮食，临床检测等等，去预测一个这个因为心脏病住院的人会不会第二次患心脏病。
基于一个公司的运行状况或一些经济数据，去预测未来六个月股票的价格。
从一个手写字母图像中识别出来其中的字母。
从一个糖尿病(diabetic)患者血液的红外吸收频谱去预测他的血糖(glucose)含量。
基于人口统计(demographic)和临床检测，分析前列腺癌的致病因素。

在一个典型的学习场景下，我们通常有一些定量的结果(outcome measurement)，如上面例子中的股票价格或者分类问题中问题的类别，我们希望基于一系列的特征进行预测。
接下来给了几个真实的学习问题的示例。下面就简要介绍一下这几个例子。

邮件分类

给定一封邮件，邮件分类的目标就是根据邮件的特征去判断这封邮件是正常邮件还是垃圾邮件。这是监督学习中的二分类问题，因为该问题有ouputs，且只有两个类别。

前列腺癌(prostate cancer)

该问题的目标是给定一系列临床检测，如记录癌症量(log cancer volume)，去预测前列腺特异性抗原(prstate specific antigen)的数量。该问题是监督学习中的回归问题，因为结果(outcome measurement)是定量的(quatitative)。

手写数字识别

给定一个手写数字的图片，该问题的目标是识别出图片中的数字。

DNA Expression Microarrays

这个问题是通过基因数组去学习基因和不同基因样本之间的关系，一些典型的问题是：

哪些样本之间是相似的？在不同的基因之间都相似。
哪些基因是相似的？在不同的样本之间都相似。
一些特定的基因对于特定的癌症患者表达是达不是很明显？

这个问题可以看成回归问题，或者更有可能是无监督问题。

本书结构

第一章就是本章。第二章讲监督学习的介绍。第三章和第四章介绍线性回归和分类。第五章介绍仿样(splines)，小波(wavelets)，正则化(regularization)和惩罚(penalization)。第六章介绍核方法(kernel methods)和局部回归(local regression)。第七章将模型估计和选择(model assessment and selection)，涉及到偏置(bias)和方差(variance)，过拟合(overfitting)以及交叉验证(cross-validation)等等。第八章讲模型推理。第十章讲boosting。
第九到十三章讲监督学习的一系列结构化方法。十四章介绍非监督学习。十五和十六章分别介绍随机森林(random forests)和集成学习(ensemble learning)。第十七章介绍无向图(undirected graphical models)。第十八章介绍高维问题。
第一到四章是基础最好按顺序阅读，第七章也是。其他的可以不按顺序。

ESL chapter 2 Overview of supervised learning

发表于 2019-01-05 | 更新于 2019-12-17 | 分类于机器学习

引言

在机器学习领域，监督学习(supervised learning)的每一个样本都由输入(inputs)和输出(outputs)组成。监督学习的目标就是根据inputs的值去预测outpus的值。
在统计学(statistical)中，inputs通常被称为预测器？？(predictors)，或者叫自变量(independent variables)。
在模式识别(pattern recognition)领域，inputs通常称为特征(features)，或者叫因变量(dependent variables)

变量类型和一些术语(terminology)

不同的问题中，输出也不一样。血糖预测问题中，输出是一个定量的(quantitative)测量。手写数字识别问题中，输出是十个不同的类，是定性的(qualitative)，定性的输出也通常被称为类别(catrgorical)，这里的类别是无序的。通常，预测定量的输出被称为回归问题(regression)，预测定性的输出被称为分类问题。这两个问题很相像，多可以看成函数拟合。第三种输出是有序类别，像小，中，大，没有合适的度量表示，因为中和小之间的差别和中和大之间的差别是不同的。
定性分析在代码实现中进行二值化数值表示。即如果只有两类的话，用一个二进制位$0$或者$1$表示，或者$1$和$-1$。当超过两类的时候，通常用虚拟变量(dummy variables)来表示，一个$K$级变量是一个长度为$K$的二进制位，每一个时刻只有一位被置一。
一些常用的表示，$X$表示inputs，$Y$表示定量outputs，$G$表示定性outputs。大写字母表示通用的表示，观测值用小写字母表示，inputs $X$的第$i$个观测值用$x_i$表示，其中$x_i$是一个标量或者向量。矩阵用粗体的大写字母表示，如具有$N$个$p$维向量$x_i, j= 1,\cdots,N$的$N\times p$矩阵$\mathbf{X}$。所有的向量都用的是列向量表示，$\mathbf{A}$的第$i$行是$x_i^T$，第$i$列的转置。

马晓鑫爱马荟荟

记录硕士三年自己的积累

RSS

GitHub E-Mail