学习参数

2018-03-16 本文已影响3人灵妍

声明：概率图模型基于R语言学习笔记

独立同分布（iid.）:每个变量都假设服从同样的概率分布，且每个观测又独立于数据集中的其它变量。

最大似然估计只能用于数据完备的情况下，如果存在数据缺失或者隐变量，就要用到EM（期望最大化）。

1、一个简单的例子引入

在这个例子中我们通过三种花的花萼长宽以及花瓣长宽判断花的类别。数据是完整的，就是我们有足够多的三种花的四个特征参数作为训练集，通过频率统计法求出概率分布。
代码：

x=read.csv("http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data",
           col.names = c("sepal_length","sepal_width","petal_length","petal_width","class"))
head(x)
mean(x$sepal_length)
var(x$sepal_length)
library(plyr)
y=daply(x,.(class),nrow)/nrow(x)
y
daply(x,.(class),function(n) mean(n$sepal_length))
daply(x,.(class),function(n) var(n$sepal_length))
q<-quantile(x$sepal_width,seq(0,1,0.33))
x$dsw[x$sepal_width<q['33%']]="small"
x$dsw[x$sepal_width>=q['33%']&x$sepal_width<q['66%']]="medium"
x$dsw[x$sepal_width>=q['66%']]="large"

p1<-daply(x,.(dsw,class),function(n) nrow(n))
p1

p1<-p1/colSums(p1)
p1

代码解释：
前面我们通过对不同类别花的行数统计得到三种花的先验概率。
中间我们假设花萼长度服从高斯分布，得出三种花花萼长度的均值与方差，从而得到后验分布密度函数。
后面我们通过对花萼宽度做离散化处理，得到分位数将样本做离散化处理，均分为三份，统计不同种类花的大小行数并归一化处理。
执行结果：

列归一化.PNG

离散化.PNG

花萼长均值及方差.PNG

类型先验分布.PNG

mean$var.PNG

head.PNG

2、作为推断学习参数

前面我们看到学习任务可以通过频率主义的计数法来完成。在多数情况下，这已经足够了，但这还只是机器学习的狭窄认识。更普遍的讲，学习的过程是把数据与领域知识整合在一起的过程，以便创建新的模型或改进现有的模型。

符号： CodeCogsEqn (5).gif 代表正比于

平板表示法：

平板表示法.png

在这个例子中，我们已知硬币是正面不公平，反面不公平，公平的先验概率，通过后面的实验序列得到后验概率更新先验概率，也就是置信度。
我们要求的是在已知序列的条件下，参数的后验概率。公式推导比较简单，这里不做赘述。
程序：

posterior<-function(prob,nh,nt,Theta=c(0.2,0.5,0.8))
{
  x=numeric(3)
  for(i in 1:3)
    x[i]=prob[i]*(Theta[i]^nh)*((1-Theta[i])^nt)
  norm=sum(x)
  return(x/norm)
}
posterior(c(0.2,0.75,0.05),2,8)
posterior(c(0.2,0.75,0.05),8,2)
posterior(c(0.2,0.75,0.05),5,5)
posterior(c(0.2,0.75,0.05),10,10)
posterior(c(0.2,0.75,0.05),50,50)

posterior(c(1/3,1/3,1/3),2,8,c(0.2,0.5,0.8))
posterior(c(1/3,1/3,1/3),8,2,c(0.2,0.5,0.8))
posterior(c(1/3,1/3,1/3),5,5,c(0.2,0.5,0.8))

代码解释：
函数中第一个参数代表先验概率，第二个参数代表正面朝上的次数，第二个采纳数代表反面朝上的次数，第三个参数硬币三种可能的分类结果，以及在此分类结果下正面朝上的概率。程序的目的就是算出后验概率更新先验概率。
运行结果：

先验是均匀分布的检测结果.PNG

不同证据的参数更新.PNG

当你需要把许多小数值相乘时，使用对数的加和，而不是原始值和乘法

3、最大似然概率

最大似然估计的目的是要找出参数的值，最大化似然率如果想更准确的刻画，我们可以采用贝叶斯方法，也给出的先验概率分布，在这个例子中找出参数值可以分解为找出的最大值。这个过程叫最大化后验概率。

Kullback-Leibler散度（也叫作相对熵）是对两种概率分布q和p差异的非对称度量，记作KL（q|p）。它给出了从q中样本到p中样本转化所需的比特数。直觉上讲，如果两个分布是一样的，那么Kullback-Leibler散度是0。
最大化似然率等价于最小化对数似然率，最小化对数似然率也会最小化经验分布q和模型分布p之间的KL散度。这只是意味着，找出p(x)最大化似然参数等价于最小化经验分布和模型分布之间的KL散度。

代码：

library(graph)
library(Rgraphviz)
library(plyr)

data0 <- data.frame(
    x=c("a","a","a","a","b","b","b","b"),
    y=c("t","t","u","u","t","t","u","u"),
    z=c("c","d","c","d","c","d","c","d"))

edges0 <- list(x=list(edges=2),y=list(edges=3),z=list())
g0 <- graphNEL(nodes=names(data0),edgeL=edges0,edgemod="directed") 
plot(g0)

data1 <- read.csv("http://archive.ics.uci.edu/ml/machine-learning-databases/nursery/nursery.data", col.names=c("parents","has_nurs","form","children","housing","finance","social","health","class"))
edges1 <- list( parents=list(), has_nurs=list(), form=list(), children=list(), 
         housing=list(), finance=list(), social=list(), health=list(),
         class=list(edges=1:8) )
g1 <- graphNEL(nodes=names(data1), edgeL=edges1,edgemod="directed")
plot(g1)

make_cpt<-function(df,pa)
{
    prob <- nrow(df)
    parents <- data.frame(df[1,pa])
    names(parents) <- pa

    data.frame(parents,prob)
}

learn <- function(g,data)
{
    rg <- reverseEdgeDirections(g)
    result <- list()

    for(var in rg@nodes)
    {
        pa <- unlist(adj(rg,var))
        if(length(pa)>0)
        {
            X <- ddply(data, c(var,pa), make_cpt, pa)
            Y <- ddply(data, pa, make_cpt, pa)
            for(i in 1:nrow(Y))
            {
                c <- sapply(1:nrow(X), function(j) all(X[j,pa] == Y[i,pa]))
                c <- which(c)
                X$prob[c] <- X$prob[c]/Y$prob[i]
            }
        }
        else
        {
            X <- ddply(data,var, function(df) c(prob=nrow(df)))
            X$prob <- X$prob/sum(X$prob)
        }

        result[[length(result)+1]] <- X
    }

    return(result)
}
learn(g0,data0)
learn(g1,data1)

代码解释：
第一个函数中计算通过输入某一节点及其父节点计算它的条件概率分布表。在程序中我们需要初始化树，然后在循环结构中计算每一个节点的概率分布表，这一过程通过判断其是否拥有父节点实现。
运行结果：