概率估计与Loss function

决策与概率估计

假设一个很简单的实验来看我们做决策的过程，有一个有偏的硬币，投了100次，60次出现人头，40次出现字，下一次抛硬币，如果我们猜对我们得1元，猜错我们交出1元。很直觉的，我们每次都会猜人头的啦。细细一想，我们在做这样一件事，max(60%*1+40%*(-1), 60%*(-1)+40%*1) 就是在估计了概率的情况下，我们看两个决策（人头/字）哪个带来的收益大。

换一个想法，我们不先估计出经验分布P(Y|X)，而是直接看有没有可能做出一个决策，使得它在历史的数据上得到收益是最大的，即max(60*1+40*(-1), 60*(-1)+40%*1)。这两种方法可以称为概率估计法及决策函数法。

这个问题的两种解决方案都有明确的最优解，但实际情况是决策函数无法直接优化（上面的例子的事件服从伯努利分布，优化比较简单，但常常我们遇到的问题是包含像sign()这样的函数）， 所以我们会用代理损失函数来做. 所谓代理损失函数就是用一个高阶可微的损失函数去接近决策的损失函数。但是有多接近呢，谁也不知道。于是很多人转向另一个阵营：更精确的估计出条件概率分布，然后做决策。如果能精确估计出条件概率，当然很好，决策的问题也迎刃而解，但事实又不是，不考虑决策损失直接做概率分布，即使是适定的，也会有问题，这个后面再说。

这里介绍的是第二种方案，即条件概率分布的估计。对于决策而言，这是个两步的过程，1）忽视Risk估概率 2）有了概率算对应的决策对应的Risk。2）比较简单，我们又focus在第一步，即“估概率”。

概率估计常见的loss

估算概率常见的思路是最大似然或是最小均方误差，这两者形式上非常的不相同，甚至很多人认为最小均方误差不适于做概率估计。原因有二，1）概率用和0及1的距离平方和做误差有问题，P=0.4和P=0.5之间的误差和0.1与P=0.001能一样吗？ 2）使均方误差最小的参数组合，会让某些(X,Y) samples小于0或大于1。这两个直觉有道理吗？有道理，但这是直觉，能把你的疑问说清楚吗？为什么你会觉得P=0.4和P=0.5之间的误差和0.1与P=0.001之间的误差不一样？

由于一个概率分布P(Y|X)在每一点X上都可以看成伯努利分布，如果我们所假设的P(Y|X)的形式足够复杂，这个条件概率分布函数的估计就转化为每一点上的伯努利分布的估什，所以我们还是举P(Y|X)为伯努利的例子看，设真实的概率为p,估计值为\hat{p}，对最大似然（logistic loss）而言，损失函数是

-p\log{\hat{p}}-(1-p)\log(1-\hat{p})

, 对最小均方误差而言，损失函数是

p(\hat{p}-1)^2+(1-p)(1-\hat{p})^2

一个显见的事实是，这两个损失函数都在\hat{p}=p时取得最值，并且是凸的。既然能保证期望风险下经数值优化后最优的估计概率等于真实概率，二者都能达到相同的目的。这正是我们设计损失函数的目的：能正确地估计出概率分布函数？这就是最近很多人在研究的两类和多类损失函数的设计，关于这个主题可以参考"On the Design of Loss Functions for Classification"及Savage有篇老文章“Elicition of Personal Probabilities”，最近有一篇关于多类问题的引申，可以看"composite multiclass loss"。

满足上一节所说的两个条件的loss是适定的(proper)。其它的一些loss像是hinge loss就不满足，只是对Decision Risk的直接优化，即构造surrogate loss（当然SVM最初并没有从surrogate loss的角度来设计）。下面我们只讨论lse和logistic loss的区别，用一个实际的例子来说明。

两个loss有何区别

lse和logistic loss就伯努利分布的估计而言是没有任何区别。我们现在稍加扩展，假定我们想要估计的P(Y|X)符合

P(Y|X)=\frac{1}{(1+\exp{(-\lambda_1 x-\lambda_0)})}

而已知的训练数据是我们人为设计的一组略有矛盾的数据,使得上面的条件分布函数无法精确描述数据。我们这样来构造这组数据，固定\lambda_0=0 使得\lambda_1=0.5及 \lambda_1=1.0

我们有下面的表格

\lambda_1	x=0.0	x=-1.0	x=-2.0
0.5	0.5	0.3775	0.2689
1.0	0.5	0.2689	0.1192

按上面表格的每一行P(Y|X)都能精确的估计参数，但我们选了加黑部分的P(Y|x)来生成训练数据。每个单元格对应的概率值生成1000个样本。生成如下表格的样本

P(Y|X)	正例数量	负例数量
0.5	500	500
0.3775	378	622
0.1192	119	881

在这个训练集合上我们可以看看用不同的loss训练出来的参数有什么特点。Logistic loss训练出的参数是0.841，LSE训练出的参数是0.795。的确有差别，在这种简单的情况之下，尝试分析一下原因。我们画出-p\log{\hat{p}}-(1-p)\log(1-\hat{p})在p=0.1192及p=0.500时的曲线，也许你能找到一些线索。这个解释可以作为一个思维练习一下。

训练的代码可以用最近比较火的theano库写，支持符号计算，不需要写梯度函数。

import re,sys,os
import time
import numpy

import theano
import theano.tensor as T

data_path = sys.argv[1]
max_feat_num = 1

X_var = T.dmatrix()
Y_var = T.ivector()

s_Wb = theano.shared(numpy.zeros(max_feat_num+1, dtype='float64'), name='s_Wb')
s_W = s_Wb[:max_feat_num].reshape((max_feat_num,1))
s_b = s_Wb[max_feat_num:]
#exp_linear = T.dot(X_var, s_W)+s_b
exp_linear = T.dot(X_var, s_W)
exp_py_given_x = T.nnet.sigmoid(exp_linear);#p
exp_npy_given_x = T.nnet.sigmoid(-exp_linear);# 1-p (softmax and sigmoid)

# loss
# 1. logistic loss
exp_py_npy = T.concatenate((exp_npy_given_x, exp_py_given_x), axis = 1)
exp_logistic = -T.mean(T.log(exp_py_npy)[T.arange(Y_var.shape[0]), Y_var])

# 2. LSEI_p = (exp_py_given_x-1)*(exp_py_given_x-1)
I_n = exp_py_given_x*exp_py_given_x
exp_lse = T.mean(T.concatenate((I_n, I_p), axis = 1)[T.arange(Y_var.shape[0]), Y_var])

# gradientexp_logistic_grad = T.grad(exp_logistic, s_Wb)
exp_lse_grad = T.grad(exp_lse, s_Wb)

X=[]
Y=[]
data_file = open(data_path, "r")
for line in data_file:
  ss = re.split("\s+", line.strip())
  Y.append((int(ss[0])+1)/2)
  x = [0.0]*max_feat_num;
for i in range(1, len(ss)):
    ss1 = re.split(":", ss[i])
    x[int(ss1[0])] = float(ss1[1])
  X.append(x)
data_file.close()

# try to use logistic loss or adaptive loss or mse loss
loss_and_grad = [exp_lse, exp_lse_grad]
s_X = theano.shared(numpy.asarray(X, dtype="float64"), name='s_X')
s_Y = theano.shared(numpy.asarray(Y, dtype="int32"), name='s_Y')
func_cost = theano.function([], loss_and_grad[0], givens={
  X_var: s_X,
  Y_var: s_Y
  }, name='func_cost');
func_grad = theano.function([], loss_and_grad[1], givens={
  X_var: s_X,
  Y_var: s_Y
  }, name='func_grad');

def cost(wb):
  s_Wb.set_value(wb, borrow=True)
  cost = func_cost()
print cost
return cost

def grad(wb):
  s_Wb.set_value(wb, borrow=True)
return func_grad()

def cb(wb):
  s_Wb.set_value(wb, borrow=True)
print wb
# try solve the optimization through lbfgs
import scipy.optimize
best_w_b = scipy.optimize.fmin_l_bfgs_b(
  func=cost,
  x0=numpy.zeros(max_feat_num+1, dtype="float64"),
  fprime=grad,
  iprint=1,
  maxfun=100)
print best_w_b

设计你所需要的loss

上面两种是最常见的两种proper的loss。在实际应用中，我们如果将loss design和decision risk完全隔裂，就会产生一种观点，只需要一个loss是proper的，估完probability再做decision.但我们往往需要将decision risk本身考虑到conditional probability estimation的过程中。因为有时把某个区域内的条件概率估对(或估错)对于最后的decision making 没有任何的帮助（或损失）。这种因素如果我们不在probability estimation时就考虑，就让训练专注于那些本不该专注的地方。如何把decision risk考虑到probability estimation中is art rather than science。所以可能需要更多的工程经验和想象力才行。