hangscer

scala——figaro概率编程入门

2017/09/14

《概率论与数理统计》是学生的必修课,本人隐隐约约还能想起个一二三。figaro(费加罗)是基于scala语言的概率编程软件,你可以在这里获取相关信息。虽说数学思想才是核心,但是要借助计算机得以计算大规模数据。本博文意在说明数学原理与算法的同时,还要详细讲解figaro库的使用方法(“调参侠”不可取)。
概率推理是领域模型在不确定性条件下做出决策的一种方法。这里以角球模型来讲解:有统计显示9%的角球会造成进球,攻方中锋身高1.9米、以头球著称,守方门将刚刚受伤,被仅仅第一次出场的替补门将换上,而且球场的风很大。现在来简要说明三种推理方式:

  • 预测未来事件,如下图所示,根据当前情况预测是否进球,证据一般包括中锋身高、守门员经验和风力。

  • 推断事件的原因,现在快进20秒后,发现进球了,请评估这位守门员的水平如何?用于预测未来事件的模型同样也可以用于事后推断结果的根源。证据与以前一样结合进球得分这一事实,以此来评估守门员水平。模型本身遵守时间顺序,但是推理可以向前和向后。

  • 从过去的事件中学习,更好的预测未来事件,再次快进10分钟,发现又有一次角球机会,现在与前一次角球的情况相类似,请评估此次角球进球的概率。证据来自上一次进球的情况。

说明:
推理系统得到数据越多,预测越准确。预测的质量取决于两点:原始模型的准确度和提供数据量。提供的数据量越多,原始模型越不重要,系统会去忘掉那些不那么重要的原始模型。如果提供的数据量越小,那么只有在原始模型十分准确的时候,才能作出确切的推理。

hello world

现在有这样的概率分布:

上图说明:今天有20%的可能性放晴,在晴天的时候,问候语有60%的可能性是“hello”,如果今天是晴天,那么明天有80%的可能性是晴天。
问题:

  • 如果今天的问候语是“hello”的概率是多少?
  • 如果今天的问候语是“hello”,那么今天是晴天的概率是多少?
  • 如果今天的问候语是“hello”,那么明天的问候语也是“hello”的概率是多少?

在这里我们先说明贝叶斯公式(贝叶斯方法十分复杂,三言两语说不清,这里只是说明其形式):

  • 条件概率:事件A在事件B已经发生条件下的发生概率,记为P(A|B)
  • 联合概率:两件事共同发生的概率,记为P(A∩B)
  • 边缘概率(先验概率):在联合概率中,把不需要的事件合并为它的全概率,从而消去这些无关变量。事件A的的边缘概率记为P(A)

解答以上的三个问题
记今天为晴天为事件A,今天问候语为“hello”为事件B

  • P(B)=0.2*0.6+0.8*0.2
  • 第二个问题需要贝叶斯公式,在今天的问候语是“hello”的情况下,去计算今天是晴天的概率:P(A|B)=(0.6)*0.2/(0.2*0.6+0.8*0.2)
  • 第三个问题先计算出今天晴天的概率(已由问题二给出答案),然后计算明天晴天的概率(根据链式法则)

解答完毕,接下来采用编程方式解决以上三个问题。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
import com.cra.figaro.language.{Flip,Select}
import com.cra.figaro.library.compound.If
import com.cra.figaro.algorithm.factored.VariableElimination
object Main{
val sunnyToday=Flip(0.2)
val greetingToday=If(sunnyToday,Select(0.6->"hello",0.4->"no"),Select(0.2->"hello",0.8->"no"))
val sunnyTomorrow=If(sunnyToday,Flip(0.8),Flip(0.05))
val greetingTomorrow=If(sunnyTomorrow,Select(0.6->"hello",0.4->"no"),Select(0.2->"hello",0.8->"no"))
//以上是定义模型部分
def predict()={
val result=VariableElimination.probability(greetingToday,"hello")
println("Today's greeting is hello with probability "+result)
}
def infer()={
greetingToday.observe("hello")
val result=VariableElimination.probability(sunnyToday,true)
println("If today's greeting is hello ,today's weather is sunny with probaility :"+result)
}
def learnAndPredict()={
greetingToday.observe("hello")
val result=VariableElimination.probability(greetingTomorrow,"hello")
println("If today's greeting is hello,tommorrow's greeting is hello with probability :"+result)
}
def main(args: Array[String]): Unit = {
predict() //预测
infer() //推理
learnAndPredict() //学习与预测
}
}
// output
Today's greeting is hello with probability 0.28
If today's greeting is hello ,today's weather is sunny with probaility :0.4285714285714285
If today's greeting is hello,tommorrow's greeting is hello with probability :0.3485714285714286

程序解释:

1
val sunnyToday=Flip(0.2)

,这句话定义了一个变量,代表了产生true值的概率为0.2、产生false值的概率为0.8的随机过程。因此,Flip(0.2)是一个可能取值为true和false的元素。

接下来,使用IfSelect构建更复杂的元素:

1
val greetingToday=If(sunnyToday,Select(0.6->"hello",0.4->"no"),Select(0.2->"hello",0.8->"no"))

,这个元素代表着一个随机过程,如果sunnyToday为true,那么选择”hello”的概率为0.6而且”on”的概率为0.4,如果sunnyToday为false,就选择后者。

假设今天已经观测到今天的问候语是”hello”,那么可为系统提供这一证据:

1
greetingToday.observe("hello")

然后再查询今天是晴天的概率:

1
2
3
val result=VariableElimination.probability(sunnyToday,true)
println(result)
// 0.4285714285714285

我们看到输出结果明显高于先前的0.2,因为今天的问候语是”hello”,那么今天是晴天的概率更高。这一推理是运用贝叶斯方法的简单实例。
移除观察值可以使用greetingToday.unobserver()方法。

在没有今天问候语证据情况下,得到明天的问候语为”hello”的概率为:

1
2
println(VariableElimination.probability(greetingTomorrow,"hello"))
//0.27999999999999997

若得到今天问候语为”hello”,则明天的问候语为”hello”的概率为:

1
2
3
greetingToday.observe("hello")
println(VariableElimination.probability(greetingTomorrow,"hello"))
// 0.3485714285714286

我们可以看到结果比之前的更大,今天的问候语是“hello”,则今天是晴天的可能性更高,那么明天是晴天的概率也很大,最终明天的问候语更可能是”hello”。推断过去是为了预测未来,figaro负责了所有计算。