Loading [MathJax]/jax/element/mml/optable/GeneralPunctuation.js

PRML-第2章 概率分布 总结

一些记号

D={x1,...,xN}
观测数据集


2.1 二元变量-伯努利分布

伯努利概率分布为:(x只能取0或1,取1的概率是μ,p(x=1|μ)=μ)

Bern(x|μ)=μx(1μ)1x


均值
E[x]=μ


方差
var[x]=μ(1μ)


似然函数

p(D|μ)=Nn=1p(xn|μ)=Nn=1μxn(1μ)1xn


对数似然函数

lnp(D|μ)=Nn=1lnp(xn|μ)=Nn=1xnlnμ+(1xn)ln(1μ)


最大似然解

μML=1NNn=1xn

μML=mN

m为x=1的观测次数


2.1 二项分布

换一个角度,在给定数据集规模N的条件下,x=1的观测出现数量m的概率分布:二项分布。
(x只能取0或1,取1的概率是μ,p(x=1|μ)=μ)

Bin(m|N,μ)=(Nm)μm(1μ)Nm

其中 (Nm)N!(Nm)!m!


均值

E[m]Nm=0mBin(mN,μ)=Nμ


方差

var[m]Nm=0(mE[m])2Bin(mN,μ)=Nμ(1μ)

2.1.1 Beta分布

注意,这里的μ是Beta分布的随机变量,不是均值

Beta(μ|a,b)=Γ(a+b)Γ(a)Γ(b)μa1(1μ)b1

其中

Γ(x)0ux1eudu


Beta函数性质

Γ(x+1)=x!


均值:
E[μ]=aa+b


方差
 var[μ]=ab(a+b)2(a+b+1)
这里a和b是超参数。


Γ分布是二项分布关于参数μ的共轭分布族,相关概念见https://www.cnblogs.com/boyknight/p/15898454.html
公式2.13 乘以公式2.9(非\mu的常数项拿掉)得到下面的正比例关系,其中l=Nm

p(μ|m,l,a,b)μm+a1(1μ)l+b1

仿照公式2.13 做归一化,得到归一化系数,则二项分布关于参数μ的后验分布如下(前提是先验是Γ分布)

p(μ|m,l,a,b)=Γ(m+a+l+b)Γ(m+a)Γ(l+b)μm+a1(1μ)l+b1

贝叶斯顺序推断 python实现

书上 图2.3

点击查看代码

import numpy as np
import matplotlib.pyplot as plt

from prml.rv import (
    Bernoulli,
    Beta,
    Categorical,
    Dirichlet,
    Gamma,
    Gaussian,
    MultivariateGaussian,
    MultivariateGaussianMixture,
    StudentsT,
    Uniform
)

np.random.seed(1234)

x = np.linspace(0, 1, 100)
beta = Beta(2, 2)
plt.subplot(2, 1, 1)
plt.xlim(0, 1)
plt.ylim(0, 2)
plt.plot(x, beta.pdf(x))
plt.annotate("prior", (0.1, 1.5))

model = Bernoulli(mu=beta)
model.fit(np.array([1]))
plt.subplot(2, 1, 2)
plt.xlim(0, 1)
plt.ylim(0, 2)
plt.plot(x, model.mu.pdf(x))
plt.annotate("posterior", (0.1, 1.5))

plt.show()


点击查看代码

import numpy as np


class RandomVariable(object):
    """
    base class for random variables
    """

    def __init__(self):
        self.parameter = {}

    def __repr__(self):
        string = f"{self.__class__.__name__}(\n"
        for key, value in self.parameter.items():
            string += (" " * 4)
            if isinstance(value, RandomVariable):
                string += f"{key}={value:8}"
            else:
                string += f"{key}={value}"
            string += "\n"
        string += ")"
        return string

    def __format__(self, indent="4"):
        indent = int(indent)
        string = f"{self.__class__.__name__}(\n"
        for key, value in self.parameter.items():
            string += (" " * indent)
            if isinstance(value, RandomVariable):
                string += f"{key}=" + value.__format__(str(indent + 4))
            else:
                string += f"{key}={value}"
            string += "\n"
        string += (" " * (indent - 4)) + ")"
        return string

    def fit(self, X, **kwargs):
        """
        estimate parameter(s) of the distribution

        Parameters
        ----------
        X : np.ndarray
            observed data
        """
        self._check_input(X)
        if hasattr(self, "_fit"):
            self._fit(X, **kwargs)
        else:
            raise NotImplementedError

    # def ml(self, X, **kwargs):
    #     """
    #     maximum likelihood estimation of the parameter(s)
    #     of the distribution given data

    #     Parameters
    #     ----------
    #     X : (sample_size, ndim) np.ndarray
    #         observed data
    #     """
    #     self._check_input(X)
    #     if hasattr(self, "_ml"):
    #         self._ml(X, **kwargs)
    #     else:
    #         raise NotImplementedError

    # def map(self, X, **kwargs):
    #     """
    #     maximum a posteriori estimation of the parameter(s)
    #     of the distribution given data

    #     Parameters
    #     ----------
    #     X : (sample_size, ndim) np.ndarray
    #         observed data
    #     """
    #     self._check_input(X)
    #     if hasattr(self, "_map"):
    #         self._map(X, **kwargs)
    #     else:
    #         raise NotImplementedError

    # def bayes(self, X, **kwargs):
    #     """
    #     bayesian estimation of the parameter(s)
    #     of the distribution given data

    #     Parameters
    #     ----------
    #     X : (sample_size, ndim) np.ndarray
    #         observed data
    #     """
    #     self._check_input(X)
    #     if hasattr(self, "_bayes"):
    #         self._bayes(X, **kwargs)
    #     else:
    #         raise NotImplementedError

    def pdf(self, X):
        """
        compute probability density function
        p(X|parameter)

        Parameters
        ----------
        X : (sample_size, ndim) np.ndarray
            input of the function

        Returns
        -------
        p : (sample_size,) np.ndarray
            value of probability density function for each input
        """
        self._check_input(X)
        if hasattr(self, "_pdf"):
            return self._pdf(X)
        else:
            raise NotImplementedError

    def draw(self, sample_size=1):
        """
        draw samples from the distribution

        Parameters
        ----------
        sample_size : int
            sample size

        Returns
        -------
        sample : (sample_size, ndim) np.ndarray
            generated samples from the distribution
        """
        assert isinstance(sample_size, int)
        if hasattr(self, "_draw"):
            return self._draw(sample_size)
        else:
            raise NotImplementedError

    def _check_input(self, X):
        assert isinstance(X, np.ndarray)


点击查看代码

import numpy as np
from scipy.special import gamma
from prml.rv.rv import RandomVariable


np.seterr(all="ignore")


class Beta(RandomVariable):
    """
    Beta distribution
    p(mu|n_ones, n_zeros)
    = gamma(n_ones + n_zeros)
      * mu^(n_ones - 1) * (1 - mu)^(n_zeros - 1)
      / gamma(n_ones) / gamma(n_zeros)
    """

    def __init__(self, n_zeros, n_ones):
        """
        construct beta distribution

        Parameters
        ----------
        n_zeros : int, float, or np.ndarray
            pseudo count of zeros
        n_ones : int, float, or np.ndarray
            pseudo count of ones
        """
        super().__init__()
        if not isinstance(n_ones, (int, float, np.number, np.ndarray)):
            raise ValueError(
                "{} is not supported for n_ones"
                .format(type(n_ones))
            )
        if not isinstance(n_zeros, (int, float, np.number, np.ndarray)):
            raise ValueError(
                "{} is not supported for n_zeros"
                .format(type(n_zeros))
            )
        n_ones = np.asarray(n_ones)
        n_zeros = np.asarray(n_zeros)
        if n_ones.shape != n_zeros.shape:
            raise ValueError(
                "the sizes of the arrays don't match: {}, {}"
                .format(n_ones.shape, n_zeros.shape)
            )
        self.n_ones = n_ones
        self.n_zeros = n_zeros

    @property
    def ndim(self):
        return self.n_ones.ndim

    @property
    def size(self):
        return self.n_ones.size

    @property
    def shape(self):
        return self.n_ones.shape

    def _pdf(self, mu):
        return (
            gamma(self.n_ones + self.n_zeros)
            * np.power(mu, self.n_ones - 1)
            * np.power(1 - mu, self.n_zeros - 1)
            / gamma(self.n_ones)
            / gamma(self.n_zeros)
        )

    def _draw(self, sample_size=1):
        return np.random.beta(
            self.n_ones, self.n_zeros, size=(sample_size,) + self.shape
        )


点击查看代码

import numpy as np
from prml.rv.rv import RandomVariable
from prml.rv.beta import Beta


class Bernoulli(RandomVariable):
    """
    Bernoulli distribution
    p(x|mu) = mu^x (1 - mu)^(1 - x)
    """

    def __init__(self, mu=None):
        """
        construct Bernoulli distribution

        Parameters
        ----------
        mu : np.ndarray or Beta
            probability of value 1 for each elzement
        """
        super().__init__()
        self.mu = mu

    @property
    def mu(self):
        return self.parameter["mu"]

    @mu.setter
    def mu(self, mu):
        if isinstance(mu, (int, float, np.number)):
            if mu > 1 or mu < 0:
                raise ValueError(f"mu must be in [0, 1], not {mu}")
            self.parameter["mu"] = np.asarray(mu)
        elif isinstance(mu, np.ndarray):
            if (mu > 1).any() or (mu < 0).any():
                raise ValueError("mu must be in [0, 1]")
            self.parameter["mu"] = mu
        elif isinstance(mu, Beta):
            self.parameter["mu"] = mu
        else:
            if mu is not None:
                raise TypeError(f"{type(mu)} is not supported for mu")
            self.parameter["mu"] = None

    @property
    def ndim(self):
        if hasattr(self.mu, "ndim"):
            return self.mu.ndim
        else:
            return None

    @property
    def size(self):
        if hasattr(self.mu, "size"):
            return self.mu.size
        else:
            return None

    @property
    def shape(self):
        if hasattr(self.mu, "shape"):
            return self.mu.shape
        else:
            return None

    def _fit(self, X):
        if isinstance(self.mu, Beta):
            self._bayes(X)
        elif isinstance(self.mu, RandomVariable):
            raise NotImplementedError
        else:
            self._ml(X)

    def _ml(self, X):
        n_zeros = np.count_nonzero((X == 0).astype(np.int))
        n_ones = np.count_nonzero((X == 1).astype(np.int))
        assert X.size == n_zeros + n_ones, (
            "{X.size} is not equal to {n_zeros} plus {n_ones}"
        )
        self.mu = np.mean(X, axis=0)

    def _map(self, X):
        assert isinstance(self.mu, Beta)
        assert X.shape[1:] == self.mu.shape
        n_ones = (X == 1).sum(axis=0)
        n_zeros = (X == 0).sum(axis=0)
        assert X.size == n_zeros.sum() + n_ones.sum(), (
            f"{X.size} is not equal to {n_zeros} plus {n_ones}"
        )
        n_ones = n_ones + self.mu.n_ones
        n_zeros = n_zeros + self.mu.n_zeros
        self.prob = (n_ones - 1) / (n_ones + n_zeros - 2)

    def _bayes(self, X):
        assert isinstance(self.mu, Beta)
        assert X.shape[1:] == self.mu.shape
        n_ones = (X == 1).sum(axis=0)
        n_zeros = (X == 0).sum(axis=0)

        assert X.size == n_zeros.sum() + n_ones.sum(), (
            "input X must only has 0 or 1"
        )
        self.mu.n_zeros += n_zeros
        self.mu.n_ones += n_ones

    def _pdf(self, X):
        assert isinstance(self.mu, np.ndarray)  # 这里貌似漏了点
        return np.prod(
            self.mu ** X * (1 - self.mu) ** (1 - X)
        )

    def _draw(self, sample_size=1):
        if isinstance(self.mu, np.ndarray):
            return (
                    self.mu > np.random.uniform(size=(sample_size,) + self.shape)
            ).astype(np.int)
        elif isinstance(self.mu, Beta):
            return (
                    self.mu.n_ones / (self.mu.n_ones + self.mu.n_zeros)
                    > np.random.uniform(size=(sample_size,) + self.shape)
            ).astype(np.int)
        elif isinstance(self.mu, RandomVariable):
            return (
                    self.mu.draw(sample_size)
                    > np.random.uniform(size=(sample_size,) + self.shape)
            )


直接采用贝叶斯方法写出后验公式如下(跳过对自然参数的估计,直接采用加法法则和乘法法则),证明方法见公式1.68

p(x=1|D)=10p(x=1|μ)p(μ|D)dμ=10μp(μ|D)dμ=E[μ|D]


根据公式 2.18,2.13,2.15可知,因为2.19就是后验分布,服从2.18公式,所以,其均值就是2.15所示

p(x=1|D)=m+am+a+l+b


m,l趋近于无穷大,2.20(贝叶斯后验),2.8(最大似然估计)两者会趋近于统一,不仅仅在beta分布下成立,在其他分布也有这样性质)

μML=mN


a,b,var(μ)0


用频率派的角度解释,推导见 https://www.cnblogs.com/boyknight/p/16010734.html

Eθ[θ]=ED[Eθ[θ|D]]

其中
Eθ[θ]=p(θ)θdθ

ED[Eθ[θD]]{θp(θD)dθ}p(D)dD

θ的后验均值(在产生数据集的分布上的平均)等于θ的先验均值。

同样的我们可以得到:

varθ[θ]=ED[varθ[θ|D]]+varD[Eθ[θ|D]]

公式(2.24)中左边是$ \theta \theta \theta \theta $的后验方差小于先验方差。

2.2 多项式变量

二元变量:2个状态。推广到k个互斥状态。用one-hot表示。比如K=6,x3=1

x=(0,0,1,0,0,0)T

向量满足Kk=1xk=1。如果用参数μk来标记xk=1的概率,那么我们就得到x的分布:

p(x|μ)=Kk=1μxkk

其中μ=(μ1,...,μK)T,由于参数μk表示概率,所以需要满足$\mu_k \geq 0 \sum_k\mu_k = 1 $。公式(2.26)分布可以看作伯努利分布在多于两种输出时的泛化。很容易证明这个分布是标准化的。

xp(x|μ)=Kk=1μk=1

E[x|μ]=xp(x|μ)x=(μ1,...,μM)T=μ

考虑一个有N个独立观测值$x_1,...,x_N D$。其对应的似然函数的形式为

p(D|μ)=Nn=1Kk=1μxnkk=Kk=1μ(nxnk)k=Kk=1μmkk

令:

mk=nxnk
它表示观测到xk=1的次数。这些别称为这个分布的充分统计量。

求最大似然解,我们需要在μk的和等于1的约束下,关于μk最大化lnp(D|μ)。这可以通过拉格朗日乘数法得到

μk=mk/λ

把公式(2.32)代入限制条件$ \sum_k\mu_k = 1 \lambda = -N $。所以我们的最大似然解:

μMLk=mkN

就是观测xk=1出现占总观测的比例。

2.2.1 狄利克雷分布

多项式分布关于参数μk的共轭分布族是狄利克雷分布

Dir(μ|α)=Γ(α0)Γ(α1)...Γ(αK)Kk=1μαk1k

α0=Kk=1αk

用似然2.34乘先验2.38就得到后验,形式为:

p(μ|D,α)p(D|μ)p(μ|α)Kk=1μαk+mk1k

因为形式与先验相同,对比写出归一化系数:

p(μ|D,α)=Dir(μ|α+m) =Γ(α0+N)Γ(α1+m1)...+Γ(αK+mK)Kk=1μαk+mk1k

其中m=(m1,...,mK)T。与二项分布的beta先验一样,可以把狄利克雷分布参数αk当成观测到xk=1的数量。而二元变量就是多项式变量的一个特例。

2.3 高斯分布

单变量x,

N(x|μ,σ2)=1(2πσ2)12exp{12σ2(xμ)2}

D维向量x,

N(x|μ,Σ)=1(2π)D/21|Σ|1/2exp{12(xμ)TΣ1(xμ)}

・一组随机变量之和,概率分布随着项(随机变量个数)增加趋近于高斯分布。(拉普拉斯中心极限定理)比如均匀分布还有之前的二项分布

中心极限定理python

点击查看代码

# 均匀分布
uniform = Uniform(low=0, high=1)
plt.figure(figsize=(10, 5))

plt.subplot(1, 3, 1)
plt.xlim(0, 1)
plt.ylim(0, 5)
plt.annotate("N=1", (0.1, 4.5))
plt.hist(uniform.draw(100000), bins=20, density=True)

plt.subplot(1, 3, 2)
plt.xlim(0, 1)
plt.ylim(0, 5)
plt.annotate("N=2", (0.1, 4.5))
plt.hist(0.5 * (uniform.draw(100000) + uniform.draw(100000)), bins=20, density=True)

plt.subplot(1, 3, 3)
plt.xlim(0, 1)
plt.ylim(0, 5)
sample = 0
for _ in range(10):
    sample = sample + uniform.draw(100000)
plt.annotate("N=10", (0.1, 4.5))
plt.hist(sample * 0.1, bins=20, density=True)

plt.show()

点击查看代码

import numpy as np
from prml.rv.rv import RandomVariable


class Uniform(RandomVariable):
    """
    Uniform distribution
    p(x|a, b)
    = 1 / ((b_0 - a_0) * (b_1 - a_1)) if a <= x <= b else 0
    """

    def __init__(self, low, high):
        """
        construct uniform distribution

        Parameters
        ----------
        low : int, float, or np.ndarray
            lower boundary
        high : int, float, or np.ndarray
            higher boundary
        """
        super().__init__()
        low = np.asarray(low)
        high = np.asarray(high)
        assert low.shape == high.shape
        assert (low <= high).all()
        self.low = low
        self.high = high
        self.value = 1 / np.prod(high - low)

    @property
    def low(self):
        return self.parameter["low"]

    @low.setter
    def low(self, low):
        self.parameter["low"] = low

    @property
    def high(self):
        return self.parameter["high"]

    @high.setter
    def high(self, high):
        self.parameter["high"] = high

    @property
    def ndim(self):
        return self.low.ndim

    @property
    def size(self):
        return self.low.size

    @property
    def shape(self):
        return self.low.shape

    @property
    def mean(self):
        return 0.5 * (self.low + self.high)

    def _pdf(self, X):
        higher = np.logical_and.reduce(X >= self.low, 1)
        lower = np.logical_and.reduce(X <= self.high, 1)
        return self.value * np.logical_and(higher, lower)

    def _draw(self, sample_size=1):
        u01 = np.random.uniform(size=(sample_size,) + self.shape)
        return u01 * (self.high - self.low) + self.low



马氏距离,这里Δ就是xμ之间的马氏距离。

Δ2=(Xμ)TΣ1(Xμ)


这一段推导了高斯分布在图像上的展示,结论如下
高斯分布的图像是一个椭圆
1.Δ

2.μ1,μ2

3.λ1,λ2

|Σ|1/2=Dj=1λ1/2j


高斯分布的期望

E[X]=μ


高斯分布的协方差矩阵

var[X]=Σ


高斯分布的局限性
1.参数较多
一个协方差矩阵ΔD(D+1)2个独立参数
2.高斯分布本质上是单峰的(即只有一个最大值)

需要引入隐变量

条件高斯分布

多元高斯性质:两个变量的联合高斯分布,一个变量为条件的高斯分布也是高斯分布。边缘高斯分布也是高斯分布

p(xa,xb)=p(xa|xb)p(xb)

精度矩阵(precision matrix)

ΛΣ1


x划分为两个不相交的子集xa,xb。令xax的前M个分量,令xb为剩余的DM个分量,得到:

x=(xaxb)

μ=(μaμb)

Σ=(ΣaaΣabΣbaΣbb)

Λ=(ΛaaΛabΛbaΛbb)

条件分布p(xa|xb)由联合分布 p(x)=p(xa,xb)通过固定xb为观测到的值,然后标准化所得到的表达式就可以得到xa上的有效概率。

12(xμ)TΣ1(xμ)= 12(xaμa)TΛaa(xaμa)12(xaμa)TΛab(xbμb) 12(xbμb)TΛba(xaμa)12(xbμb)TΛbb(xbμb)

因为一个通用的的高斯分布N(x|μ,Σ)的指数项可以写成:

12(xμ)TΣ1(xμ)=12xTΣ1x+xTΣ1μ+const

p(xa|xb)的协方差(精度矩阵的逆):

Σa|b=Λ1aa

所以,p(xa|xb)的均值

μab=Σab{ΛaaμaΛab(xbμb)}=μaΛ1aaΛab(xbμb)

其中

Λaa=(ΣaaΣabΣ1bbΣba)1Λab=(ΣaaΣabΣ1bbΣba)1ΣabΣ1bb

故均值和协方差矩阵还可以表示为

μab=μa+ΣabΣ1bb(xbμb)Σab=ΣaaΣabΣ1bbΣba


总结
条件分布:

p(xa|xb)=N(x|μa|b,Λ1aa)

μa|b=μaΛ1aaΛab(xaμb)

2.32 边缘高斯分布

p(xa,xb)=p(xa|xb)p(xb)

p(xa)=p(xa,xb)dxb


E[xa]=μacov[xa]=Σaa

总结
边缘分布:

p(xa)=N(xa|μa,Σaa)

2.33 高斯变量的贝叶斯定理

本节主要求后面这5个分布的均值和方差,p(y|x)p(x)=p(x|y)p(y)=p(x,y)

假设给定高斯边缘分布p(x)和均值是关于x的线性函数且方差与$ x p(y|x)$

p(x)=N(xμ,Λ1)p(yx)=N(yAx+b,L1)

其中$ \mu, A, b \Lambda , L x,y M,D A D \times M $矩阵。


求x,y的联合分布。

z=(xy)


p(x,y)

cov[z]=R1=(Λ1Λ1ATAΛ1L1+AΛ1AT)


p(x,y)

E[z]=R1(ΛμATLbLb)

E[z]=(μAμ+b)


边缘分布y

E[y]=Aμ+bcov[y]=L1+AΛ1AT


条件分布p(xy)

E[xy]=(Λ+ATLA)1{ATL(yb)+Λμ}cov[xy]=(Λ+ATLA)1


总结:

对于$ x y x $的条件高斯分布:

p(x)=N(x|μ,Λ1)

p(y|x)=N(y|Ax+b,L1)

那么$ y x y $的条件高斯分布为:

p(y)=N(y|Aμ+b,L1+AΛ1AT)

p(x|y)=N(x|Σ{ATL(yb)+Λμ},Σ)

其中

Σ=(Λ+ATLA)1


2.3.4 高斯分布的最大似然估计

给定一个数据集X=(x1,...,xN)T,i.i.d.最大似然来估计参数,对数似然为:

lnp(X|μ,Σ)=ND2ln(2π)N2ln|Σ|12Nn=1(xnμ)TΣ1(xnμ)


μ求导(C.19):

μlnp(X|μ,Σ)=Nn=1Σ1(xnμ)


最大似然估计:

μML=1NNn=1xn


ΣML=1NNn=1(xnμML)(xnμML)T


如果我们估计真实概率分布,可以得到有偏的结果。协方差期望小于真实值。

E[μML]=μE[ΣML]=N1NΣ


所以需要补正。

˜Σ=1N1Nn=1(xnμML)(xnμML)T


2.35 顺序估计

处理一个数据,整合进模型,处理完就丢掉。
最大似然估计:

μML=1NNn=1xn

把第$ N \mu_{ML}^{(N)}$,就可以写成:

μ(N)ML=1NNn=1xn =1NxN+1NN1n=1xn =1NxN+N1Nμ(N1)ML =μ(N1)ML+1N(xNμ(N1)ML)

只要增加一个新观测数据的修正量,就可以得到结果了。随着N增加,修正量的影响也在变小


推广到通用层面:Robbins-Monro算法。

考虑一对有联合分布p(z,θ)控制的随机变量θ,z$\thetazf(\theta)$给出:

f(θ)E[z|θ]=zp(z|θ)dz

称之为回归函数。我们假定的目标是找到f(θ)=0的根θ


Robbins-Monro算法的顺序估计就是:

假设:
当$\theta > \theta^∗ f(\theta) > 0,\theta < \theta^∗ f(\theta) < 0$

θ(N)=θ(N1)+aN1z(θ(N1))

其中$ z(\theta^{(N)}) \theta \theta^{(N)} z {a_N} $表示一个满足下列条件的正数序列:

limNaN=0N=1aN=N=1a2N<

2.129以概率1收敛于根。2.130保证修正越来越小,2.131保证不会收敛到不根的值(阻止太快收束),2.132保证累计噪声是有限的(抑制noise发散),会收敛。


顺序估计应用于最大似然问题

最大似然解θML是对数似然函数的驻点。因此满足:

θ{1NNn=1lnp(xnθ)}|θML=0

交换求导与求和顺序,且令极限N得到:

limN1NNn=1θlnp(xn|θ)=Ex[θlnp(xn|θ)]

Ex[θlnp(xn|θ)]就是公式2.127的形式,故Ex=0即是最大似然解
最大似然的解就是回归函数的根

用Robbins-Monro算法:

θ(N)=θ(N1)+aN1θ(N1)lnp(xN|θ(N1))

参数aN=σ2/N


最大似然问题可以转化为顺序更新问题:N个数据点的均值:N-1个数据点的均值和一个数据点xN的贡献。

这里我们看后验分布可以写成:

p(μ|D)[p(μ)N1n=1p(xn|μ)]p(xN|μ)

2.144 是一个比较通用的顺序估计公式
方括号里的是观测N-1个数据点后的后验分布(忽略归一化系数) 可以被看作一个 先验分布


2.3.6 高斯分布的贝叶斯推断

一元高斯分布似然函数

p(xμ)=Nn=1p(xnμ)=1(2πσ2)N2exp{12σ2Nn=1(xnμ)2}


针对参数μ,高斯分布的共轭分布是高斯分布,令先验为

p(μ)=N(μ|μ0,σ20)

且后验分布由:

p(μ|X)p(X|μ)p(μ)

给出。 通过简单的配出指数中二次项的操作,可以得到的后验分布为:

p(μ|X)=N(μ|μN,σ2N)

其中

μN=σ2Nσ20+σ2μ0+Nσ20Nσ20+σ2μML1σ2N=1σ20+Nσ2

其中μMLμ的最大似然解,由样本均值给出:

μML=1NNn=1xn

μ


观察有几个结论:

  • 2.141:后验的均值是先验μ0和似然μML的折中。(2.20推导的时候也有这一点)
    如果N=0,就变成先验;N趋近无穷大,变成似然。
  • 2.142:精度是可以叠加的。每个观测精度之和+先验精度就是后验精度。
    N趋近于无穷大,后验精度趋近于零,在最大似然附近变成尖峰。
  • 2.143:当数据点无穷大,最大似然可以精确地由贝叶斯公式恢复(μN=μML)(通过后验可以计算出来?)。
  • 对于有限的N,如果σ0,先验的方差无穷大,2.141后验均值就是最大似然,2.142方差变为σ2N=σ2/N

假设均值是已知的,推断方差。 同样选择先验是共轭的。定义精度λ1/σ2进行计算是最方便的。关于λ的似然函数为:

p(xλ)=Nn=1N(xnμ,λ1)λN2exp{λ2Nn=1(xnμ)2}

因此,对应的共轭先验正比于λ的幂次数和λ的线性函数的指数。这就是Gamma分布,定义为:

Gam(λ|a,b)=1Γ(a)baλa1exp(bλ)

其中
Γ(x)0ux1eudu是归一化系数。

Gamma分布的均值和方差为:

E[λ]=abvar[λ]=ab2

设先验

Gam(λ|a0,b0)

然后假设先验为Gam(λ|a0,b0)如果乘以似然函数(2.145),那么就得到后验分布:

p(λx)λa01λN2exp{b0λλ2Nn=1(xnμ)2}

整理一下,看成Gam(λ|aN,bN)的gamma分布,其中:

aN=a0+N2bN=b0+12Nn=1(xnμ)2=b0+N2σ2ML

其中σ2ML是对方差的最大似然估计。

观察有几个结论:

  • 2.150:N个数据点的效果是使a增加了N/2。因此我们可以把先验分布中的参数a0看成2a0个“有效”先验观测。
  • 2.151:N个数据点为参数b贡献了Nσ2ML/2其中σ2ML是方差,所以把先验中的参数b0解释为:2a0个方差为2b0/(2a0)=b0/a0“有效”的先验观测的效果。

回忆一下,我们在Dirichlet先验中做过类似的有效观测数的解释。这些分布是指数族的例子,我们将会看到,把共轭先验解释为有效的虚拟数据点是指数族分布的一种通用方法。

σ2Gamma


现在,假设均值和精度都是未知的。为了找到共轭先验,考虑似然函数对μ,λ的依赖:

p(xμ,λ)=Nn=1(λ2π)12exp{λ2(xnμ)2}[λ12exp(λμ22)]Nexp{λμNn=1xnλ2Nn=1x2n}

现在,我们在想找到一个对于μ,λ的依赖与似然函数有着相同的函数形式的先验分布p(μ,λ)因此,假设形式:

p(μ,λ)[λ12exp(λμ22)]βexp{cλμdλ}=exp{βλ2(μcβ)2}λβ2exp{(dc22β)λ}

其中c,d,β是常量。由于总有p(μ,λ)=p(μ|λ)p(λ),我们可以通过观察找到p(μ|λ),p(λ)。特别的,p(μ|λ)是一个精度为关于λ的线性函数的高斯分布,p(λ)是一个gamma分布时,得到的标准化的先验形式为:

p(μ,λ)=N(μ|μ0,(βλ)1)Gam(λ|a,b)

其中,我们的新常数为
μ0=c/β,a=1+β/2,b=dc2/2β(对比就可以得到)。式(2.154)的分布被称为正态-gamma(normal-gamma)或高斯-gamma(Gaussian-gamma)分布,并在图2.14中展示。

μ,σ2gamma


即使我们选择一个μ和\lambda\mu\lambda$值也会相互耦合


对于$ D x \mathcal{N}(x|\mu, \Lambda^{−1}) \mu \Lambda $,共轭先验是Wishart分布:

W(Λ|W,v)=B|Λ|(vD1)/2exp(12Tr(W1Λ))

其中$ v W D \times D TR(.) B $为:

B(W,v)=|W|v/2(2vD/2πD(D1)/4Di=1Γ(v+1i2))1

如果均值和精度同时未知,那么,和一元变量类似的推理得到共轭先验:

p(μ,Λ|μ0,β,W,v)=N(μ|μ0,(βΛ)1)W(Λ|W,v)

这被称为正态-Wishart分布或高斯-Wishart分布。

μWishart


2.3.7 学生t分布

T分布有一个重要的性质,鲁棒性笔记好,在有一些异常值的情况下,正态分布会产生较大的偏离,但是t分布不会

点击查看代码

X = np.random.normal(size=20)
X = np.concatenate([X, np.random.normal(loc=20., size=3)])
plt.hist(X.ravel(), bins=50, density=1., label="samples")

students_t = StudentsT()
gaussian = Gaussian()

gaussian.fit(X)
students_t.fit(X)

print(gaussian)
print(students_t)

x = np.linspace(-5, 25, 1000)
plt.plot(x, students_t.pdf(x), label="student's t", linewidth=2)
plt.plot(x, gaussian.pdf(x), label="gaussian", linewidth=2)
plt.legend()
plt.show()

2.3.9 高斯混合模型

使用kmeans聚类的高斯混合模型

点击查看代码

import numpy as np
import matplotlib.pyplot as plt

from prml.rv import (
    Bernoulli,
    Beta,
    Categorical,
    Dirichlet,
    Gamma,
    Gaussian,
    MultivariateGaussian,
    MultivariateGaussianMixture,
    StudentsT,
    Uniform
)

# 例子
x1 = np.random.normal(size=(100, 2))
x1 += np.array([-5, -5])
x2 = np.random.normal(size=(100, 2))
x2 += np.array([5, -5])
x3 = np.random.normal(size=(100, 2))
x3 += np.array([0, 5])
X = np.vstack((x1, x2, x3))

model = MultivariateGaussianMixture(n_components=3)
model.fit(X)
print(model)

x_test, y_test = np.meshgrid(np.linspace(-10, 10, 100), np.linspace(-10, 10, 100))
X_test = np.array([x_test, y_test]).reshape(2, -1).transpose()
probs = model.pdf(X_test)
Probs = probs.reshape(100, 100)
plt.scatter(X[:, 0], X[:, 1])
plt.contour(x_test, y_test, Probs)
plt.xlim(-10, 10)
plt.ylim(-10, 10)
plt.show()


点击查看代码

import numpy as np
from prml.clustering import KMeans
from prml.rv.rv import RandomVariable


class MultivariateGaussianMixture(RandomVariable):
    """
    p(x|mu, L, pi(coef))
    = sum_k pi_k N(x|mu_k, L_k^-1)
    """

    def __init__(self,
                 n_components,
                 mu=None,
                 cov=None,
                 tau=None,
                 coef=None):
        """
        construct mixture of Gaussians

        Parameters
        ----------
        n_components : int
            number of gaussian component
        mu : (n_components, ndim) np.ndarray
            mean parameter of each gaussian component
        cov : (n_components, ndim, ndim) np.ndarray
            variance parameter of each gaussian component
        tau : (n_components, ndim, ndim) np.ndarray
            precision parameter of each gaussian component
        coef : (n_components,) np.ndarray
            mixing coefficients
        """
        super().__init__()
        assert isinstance(n_components, int)
        self.n_components = n_components
        self.mu = mu
        if cov is not None and tau is not None:
            raise ValueError("Cannot assign both cov and tau at a time")
        elif cov is not None:
            self.cov = cov
        elif tau is not None:
            self.tau = tau
        else:
            self.cov = None
            self.tau = None
        self.coef = coef

    @property
    def mu(self):
        return self.parameter["mu"]

    @mu.setter
    def mu(self, mu):
        if isinstance(mu, np.ndarray):
            assert mu.ndim == 2
            assert np.size(mu, 0) == self.n_components
            self.ndim = np.size(mu, 1)
            self.parameter["mu"] = mu
        elif mu is None:
            self.parameter["mu"] = None
        else:
            raise TypeError("mu must be either np.ndarray or None")

    @property
    def cov(self):
        return self.parameter["cov"]

    @cov.setter
    def cov(self, cov):
        if isinstance(cov, np.ndarray):
            assert cov.shape == (self.n_components, self.ndim, self.ndim)
            self._tau = np.linalg.inv(cov)
            self.parameter["cov"] = cov
        elif cov is None:
            self.parameter["cov"] = None
            self._tau = None
        else:
            raise TypeError("cov must be either np.ndarray or None")

    @property
    def tau(self):
        return self._tau

    @tau.setter
    def tau(self, tau):
        if isinstance(tau, np.ndarray):
            assert tau.shape == (self.n_components, self.ndim, self.ndim)
            self.parameter["cov"] = np.linalg.inv(tau)
            self._tau = tau
        elif tau is None:
            self.parameter["cov"] = None
            self._tau = None
        else:
            raise TypeError("tau must be either np.ndarray or None")

    @property
    def coef(self):
        return self.parameter["coef"]

    @coef.setter
    def coef(self, coef):
        if isinstance(coef, np.ndarray):
            assert coef.ndim == 1
            if np.isnan(coef).any():
                self.parameter["coef"] = np.ones(self.n_components) / self.n_components
            elif not np.allclose(coef.sum(), 1):
                raise ValueError(f"sum of coef must be equal to 1 {coef}")
            self.parameter["coef"] = coef
        elif coef is None:
            self.parameter["coef"] = None
        else:
            raise TypeError("coef must be either np.ndarray or None")

    @property
    def shape(self):
        if hasattr(self.mu, "shape"):
            return self.mu.shape[1:]
        else:
            return None

    def _gauss(self, X):
        d = X[:, None, :] - self.mu
        D_sq = np.sum(np.einsum('nki,kij->nkj', d, self.cov) * d, -1)
        return (
                np.exp(-0.5 * D_sq)
                / np.sqrt(
            np.linalg.det(self.cov) * (2 * np.pi) ** self.ndim
        )
        )

    def _fit(self, X):
        cov = np.cov(X.T)  # 计算协方差矩阵
        kmeans = KMeans(self.n_components)  # 直接用kmeans聚类
        kmeans.fit(X)
        self.mu = kmeans.centers
        self.cov = np.array([cov for _ in range(self.n_components)])
        self.coef = np.ones(self.n_components) / self.n_components
        params = np.hstack(
            (self.mu.ravel(),
             self.cov.ravel(),
             self.coef.ravel())
        )
        while True:
            stats = self._expectation(X)
            self._maximization(X, stats)
            new_params = np.hstack(
                (self.mu.ravel(),
                 self.cov.ravel(),
                 self.coef.ravel())
            )
            if np.allclose(params, new_params):
                break
            else:
                params = new_params

    def _expectation(self, X):
        resps = self.coef * self._gauss(X)
        resps /= resps.sum(axis=-1, keepdims=True)
        return resps

    def _maximization(self, X, resps):
        Nk = np.sum(resps, axis=0)
        self.coef = Nk / len(X)
        self.mu = (X.T @ resps / Nk).T
        d = X[:, None, :] - self.mu
        self.cov = np.einsum(
            'nki,nkj->kij', d, d * resps[:, :, None]) / Nk[:, None, None]

    def joint_proba(self, X):
        """
        calculate joint probability p(X, Z)

        Parameters
        ----------
        X : (sample_size, n_features) ndarray
            input data

        Returns
        -------
        joint_prob : (sample_size, n_components) ndarray
            joint probability of input and component
        """
        return self.coef * self._gauss(X)

    def _pdf(self, X):
        joint_prob = self.coef * self._gauss(X)
        return np.sum(joint_prob, axis=-1)

    def classify(self, X):
        """
        classify input
        max_z p(z|x, theta)

        Parameters
        ----------
        X : (sample_size, ndim) ndarray
            input

        Returns
        -------
        output : (sample_size,) ndarray
            corresponding cluster index
        """
        return np.argmax(self.classify_proba(X), axis=1)

    def classify_proba(self, X):
        """
        posterior probability of cluster
        p(z|x,theta)

        Parameters
        ----------
        X : (sample_size, ndim) ndarray
            input

        Returns
        -------
        output : (sample_size, n_components) ndarray
            posterior probability of cluster
        """
        return self._expectation(X)


点击查看代码

import numpy as np
from scipy.spatial.distance import cdist


class KMeans(object):

    def __init__(self, n_clusters):
        self.n_clusters = n_clusters

    def fit(self, X, iter_max=100):
        """
        perform k-means algorithm

        Parameters
        ----------
        X : (sample_size, n_features) ndarray
            input data
        iter_max : int
            maximum number of iterations

        Returns
        -------
        centers : (n_clusters, n_features) ndarray
            center of each cluster
        """
        I = np.eye(self.n_clusters)
        centers = X[np.random.choice(len(X), self.n_clusters, replace=False)]
        for _ in range(iter_max):
            prev_centers = np.copy(centers)
            D = cdist(X, centers)
            cluster_index = np.argmin(D, axis=1)
            cluster_index = I[cluster_index]
            centers = np.sum(X[:, None, :] * cluster_index[:, :, None], axis=0) / np.sum(cluster_index, axis=0)[:, None]
            if np.allclose(prev_centers, centers):
                break
        self.centers = centers

    def predict(self, X):
        """
        calculate closest cluster center index

        Parameters
        ----------
        X : (sample_size, n_features) ndarray
            input data

        Returns
        -------
        index : (sample_size,) ndarray
            indicates which cluster they belong
        """
        D = cdist(X, self.centers)
        return np.argmin(D, axis=1)


2.5 非参数化估计

  • 直方图估计
  • 核密度估计
  • K近邻估计

直方图方法

在实际应用中,直方图方法对于快速地将一维或者二维的数据可视化很有用,但是并不适用于大多数概率密度估计的应用
一个明显的问题是概率密度故连续

2.5.1 核密度估计

假设观测是D维空间未知概率分布p(x),希望估计。看包含x的小区域R,概率密度是:

P=Rp(x)dx

假设收集了R内部的K个数据点,服从二项分布:x落在区域R中被观测到,数量为K个的概率。

Bin(K|N,P)=N!K!(NK)!PK(1P)NK

使用2.11,2.12(在给定数据集规模N的条件下,x=1的观测出现数量m的概率分布的期望和方差)

E[m]Nm=0mBin(mN,μ)=Nμ

E[m]Nm=0(mE[m])2Bin(mN,μ)=Nμ(1μ)

得到落在区域内部的数据点的平均比例(mean fraction)为$ \mathbb{E}[K/N] = P 2.12 var[K/N] = P(1-P)/N $。

对于大的N值,这个分布将会在均值附近产生尖峰(方差变小),且

KNP

但是,如果同时假定区域R足够小,使得在这个区域内的概率密度p(x)大致为常数,那么就有

Pp(x)V

其中$V R $的体积。结合式(2.244)和(2.245)得到密度估计的形式:

p(x)=KNV

注意,式(2.246)的成立依赖于两个相互矛盾的假设,即区域$ R 使使 K $足够让二项分布达到尖峰。(太少就没点在区域里)

因为p(x)概率函数和N是数据点。我们有两种思路:

  • 固定K,通过数据确定V,就是k邻近算法。
  • 固定V,通过数据确定K的值,就是核方法。

可以证明在极限N下,V随N的增加而收缩,K随N的增加而增大。最终两种方法得到的概率密度都会收敛于真实的概率密度。(Duda and Hart, 1973)。


首先是核方法,我们把区域$ R x K $,定义函数

k(u)={1,|ui|12,i=1,,D0, 其他情况 

这表示一个以原点为中心的单位立方体。函数k(u)就是核函数的一个例子。从式(2.247),如果数据点$ x_n x h k((x-x_n)/h) $等于1,否则它的值为0。位于这个立方体内的数据点的总数为:

K=Nn=1k(xxnh)

把这个表达式代入式(2.246),可以得到点x处的概率密度估计

p(x)=1NNn=11hDk(xxnh)

使用函数$ k(u) N x_n N x$为中心的一个立方体。

但问题跟直方图一样,非连续性。(不是1就是0).这个是由密度估计中立方体的边界带来的。如果我们选择一个平滑的核函数,那么就可以得到一个更加光滑的模型。一个常用的选择是高斯核函数,它给出

p(x)=1NNn=11(2πh2)D2exp{


高斯核密度估计的python实现

点击查看代码

import numpy as np
import matplotlib.pyplot as plt
import math


def gauss(x):
    return (1 / math.sqrt(2 * math.pi)) * math.exp(-0.5 * (x ** 2))


def get_kde(x, data_array, bandwidth=0.1):
    N = len(data_array)
    res = 0
    if len(data_array) == 0:
        return 0
    for i in range(len(data_array)):
        res += gauss((x - data_array[i]) / bandwidth)  # 对应书上的公式 2.250 D=1
    res /= (N * bandwidth)
    return res


input_array = np.random.randn(20000).tolist()
print(input_array)
bandwidth = 1.05 * np.std(input_array) * (len(input_array) ** (-1 / 5))  # 带宽也就是书里面的h 是个超参数
x_array = np.linspace(min(input_array), max(input_array), 50)
y_array = [get_kde(x_array[i], input_array, bandwidth) for i in range(x_array.shape[0])]
# 这里可以看到,每次计算kde 都要代入整个观测数据集,运算很慢

plt.subplot(1, 2, 1)
plt.hist(input_array, bins=40, density=True)
plt.plot(x_array.tolist(), y_array, color='red', linestyle='-')

bandwidth = 0.005
x_array = np.linspace(min(input_array), max(input_array), 50)
y_array = [get_kde(x_array[i], input_array, bandwidth) for i in range(x_array.shape[0])]

plt.subplot(1, 2, 2)
plt.hist(input_array, bins=40, density=True)
plt.plot(x_array.tolist(), y_array, color='red', linestyle='-')
plt.show()


2.5.2

与之前固定$ V 然后从数据中确定 K 的值不同,我们考虑固定 K 的值然后使用数据来确定合适的

其实这部分很简单,回到公式2.246

p(x)=\frac{K}{NV}

K,N都是常数,计算V就是计算概率密度了,但是
注意,由K近邻方法得到的模型不是真实的概率密度模型,因为它在整个空间的积分是发散的

posted @   筷点雪糕侠  阅读(501)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 清华大学推出第四讲使用 DeepSeek + DeepResearch 让科研像聊天一样简单!
· 推荐几款开源且免费的 .NET MAUI 组件库
· 易语言 —— 开山篇
· 实操Deepseek接入个人知识库
· Trae初体验
点击右上角即可分享
微信分享提示