1. 对数样条模型

1.1 从密度函数引入：

假设要求的密度函数是 \(f(\cdot)\)，则可以通过 \(\hat{l}=\hat{s}+c(\hat{s})\) 的形式来估计 \(l=\log(f)\)。其中 \(c(\hat{s})\) 是归一化常数，可以用极大似然法去估计 \(\hat{s}\)。使得最终的密度估计 \(\hat{f}=\exp(\hat{l})\) 为正且积分为1。

对数样条（logspline）密度估计是将 B 样条拟合到密度函数的对数上，从而保证密度函数的非负和平滑估计，通常选用的是三次样条。因此对数样条密度估计实际上是一种非参数方法和参数方法的结合。

因此，对数样条密度估计可以写成：

设 \(f(x)\) 是未知密度函数， \(\hat{f}(x)\) 是其对数样条估计。则对数样条密度估计为：

\[ \hat{f}(x)=\exp\left(\sum_{j=1}^{\bar{j}} C_j B_j (x) \right) \]

其中 \(B_j (x)\) 是 B 样条基函数， \(C_j\) 是系数，求和是在所有 B 样条上进行的。系数 \(C_j\) 是通过最小化惩罚似然函数确定的，该函数平衡模型的拟合程度和估计密度的平滑度。

1.2 选用对数样条的两个原因

1.2.1 与核密度相比

我们首先从直方图谈起，直方图可以在两个维度上清晰地表示数据的分布密度，但其缺点也非常明显。首先，用直方图表示的密度分布并不光滑，并且区间的选择对于直方图的形态影响很大，当区间足够大的时候，甚至会产生相反的分布结果，这是不可以接受的，第三点，利用直方图法可以通过图像清晰地表示数据分布，但是一般我们只考虑了一个维度的数据，最多可以考虑平面（二维）数据的分布状况，对于更高维度的数据，直方图就显得乏力了。针对这种情况，我们可以考虑连续数据集中任意点\(x\)附近的分布情况。设区间为\(h\)，概率密度分布为\(F(\cdot)\), 对于\(x\)附近\(2h\)区域，我们有

\[ f(x)=\lim_{n \to 0} \frac {F(x+h)-F(x-h)}{2h} \]

我们把该区域的密度函数值当作\(x\)点的密度函数值，从数量分布上来讲，我们可以粗略地改写上述公式。我们设样本集总数量为\(n=N_{\text{total}}\), 落在邻域\([x-h,x+h]\)内的样本数量为\(N_{x_i \in [x-h,x+h]}\)，再将这一范围缩小，我们就得到了\(x\)点附近的大致分布情况：

\[ \hat{f}(x)=\frac{1}{2h} \lim_{n \to 0} \frac{N_{x_i}}{N_{\text{total}}} \]

记\(x_1,x_2,\dots,x_n\)为\(F(\cdot)\)上独立同分布的\(n\)个样本点，我们有

\[ \hat{f}(x)=\frac{1}{2nh} \sum_{i=x-h}^{x+h} x_i =\frac{1}{nh} \sum_i \frac{|x-x_i|}{2h}<1, h \to 0 \]

请注意到\(\hat{f}(x)\)并不光滑，但我们可以分析出，当\(h\)比较大的时候，\(x\)附近的分布会比较平缓，与总体密度分布函数的差距较大，当\(h\)比较小的时候，\(x\)附近的分布会尖锐起来，更趋近总体密度分布函数。这里我们记

\[ K(x)=\frac{1}{2} 1_{\{x<1\}} \]

带入概率密度函数我们有

\[ \hat{f}(x)=\frac{1}{nh} \sum_i K\left(\frac{|x-x_i|}{h}\right) \]

由于需要满足概率密度函数的积为1，所以

\[ \int \hat{f}(x) = \frac{1}{nh} \sum_i \int K\left(\frac{|x-x_i|}{h}\right) dx=\frac{1}{nh} \sum_i K(t)dt=\int K(t)dt \]

核密度函数的构建过程中利用极限的思想摒除了直方图来表示概率分布的缺陷，使其成为了使用最广的非参数概率密度估计方法。核密度估计主要取决于带宽的选择，其中普遍使用的方法是交叉验证法和插值法。而插值法则需要基于正态分布。此外，核密度函数的标准核密度估计也有一些局限性，首先，在一维中，当\(N\)为样本量时，没有正的核密度的平均累计平方误差衰减得比\(N^{-4/5}\)快。在更高的维度中，收敛的速度甚至更慢。其次，核密度估计在有界或半无限域中存在边界误差，因为核密度函数没有考虑有关密度函数的支撑的信息。第三，由于构造缺陷，标准核函数缺乏局部自适应，这往往导致虚假凸起的出现，并倾向于使概率密度函数的波峰和波谷变平。针对这些问题，有学者引入了高阶核被用来提高精度，但它们的缺点是不能保证提供适当的非负密度估计。有一些边界核估计和边界校正方案的论文已经解决了边界偏差。局部适应性的不足已经被自适应核估计器解决了。对于这三个问题，最系统和全面的解决方法是扩散估计器，它提高了渐近精度，自然地处理域边界，具有数据自适应平滑，并始终提供真正的概率密度，即非负并积分到单位。

在此基础上，我们应该考虑到核密度函数的边界误差问题的产生是因为并没有考虑到支集的情况，而是直接将其均匀分开然后进行极限化，再进行积分计算，在这种方法之上进行边界矫正的方案已经有很多，但都是针对某一类情况的修补，显得比较麻烦，在没有可靠的修正方案或者已有的修正方案误差较大时，我们应该考虑其他的方法。相比于此，对数样条估计就是一个很好的替代方案。

1.2.2 最大熵估计

首先引入一个熵的定义。1948年， Claude E. Shannon引入信息（熵），将其定义为离散随机事件的出现概率。一个系统越是有序，信息熵就越低；反之，一个系统越是混乱，信息熵就越高。所以说，信息熵可以被认为是系统有序化程度的一个度量。

事件的概率越小则事件的不确定性就高，即信息量越高。因此信息量函数\(f\) 关于概率\(p\)是减函数；两个独立事件所产生的信息量应等于各自的信息量之和，即:

\[ f(P1,P2)=f(P1)+f(P2) \]

同时满足这两个条件的函数\(f\)是对数函数，即：

\[ f(P)=\log(1/P)=-\log(P) \]

信息熵则为信息量的期望：

\[ H(P)=E[-\log(P)]=H(P)=-\sum_{X} P(X) \log P(X) \]

可以证明: \(0 \leq H(P) \leq \log|X|\) ，其中 \(|X|\) 为 \(X\) 的取值的个数。

最大熵Max Entropy原理:学习概率模型时，在所有可能的概率模型（即概率分布）中，熵最大的模型是最好的模型。通常还有其他已知条件来确定概率模型的集合，因此最大熵原理为：在满足已知条件的情况下，选取熵最大的模型。在满足已知条件前提下，如果没有更多的信息，则那些不确定部分都是“等可能的”。而等可能性通过熵最大化来刻画。

最大熵原理选取熵最大的模型，而决策树的划分目标选取熵最小的划分。原因在于：最大熵原理认为在满足已知条件之后，选择不确定性最大（即：不确定的部分是等可能的）的模型。也就是不应该再施加任何额外的约束。因此这是一个求最大不确定性的过程，所以选择熵最大的模型。决策树的划分目标是为了通过不断的划分从而不断的降低实例所属的类的不确定性，最终给实例一个合适的分类。因此这是一个不确定性不断减小的过程，所以选取熵最小的划分。

一种常见的约束为期望的约束：

\[ E[f(X)]=\sum_{X} P(X)f(X) = \tau \]

其中\(f\)代表随机变量\(X\)的某个函数（其结果是另一个随机变量）。

假设有\(k\)个约束条件：

\[ \begin{align*} E[f_1 (X)]=\sum_{X} P(X)f_1 (X) =\tau_1 \\ \vdots \\ E[f_k (X)]=\sum_{X} P(X)f_k (X) =\tau_k \end{align*} \]

则需要求解约束最优化问题：

\[ \begin{align*} \max_{P(X)} -\sum_{X} P(X)\log P(X) \\ s.t. \sum_{X} P(X) =1, 0 \leq P \leq 1 \\ E[f_1 (X)]=\sum_{X} P(X)f_1 (X) =\tau_1 \\ \vdots \\ E[f_k (X)]=\sum_{X} P(X)f_k (X) =\tau_k \end{align*} \]

利用拉格朗日乘子法求解：

\[ L(P)=-\sum_{X} P(X)\log P(X) +\lambda_0 (\sum_{X} P(X) -1) +\lambda_1 (\sum_{X} P(X)f_1 (X) -\tau_1) +\cdots +\lambda_k (\sum_{X} P(X)f_k (X) -\tau_k) \]

解得：

\[ P(X)=\frac{1}{Z} \exp\left(-\sum_{i=1}^{k} \lambda_i f_i (X)\right) \]

\[ Z=\sum_{X} \exp\left(-\sum_{i=1}^{k} \lambda_i f_i (X)\right) \]

例如：当有两个约束条件时：\(f_1 (X)=X\)，\(f_2 (X)=X^2\)时。表示约束了变量的期望和方差。即：\(\sum_{X} P(X)X =\tau_1, \sum_{X} P(X)X^2 =\tau_2\)。

此时有：

\[ P(X)=\frac{\exp(-\lambda_1 X-\lambda_2 X)}{\sum_{X} \exp(-\lambda_1 X-\lambda_2 X)} \]

代入约束解得：

\[ P(X)=\sqrt{\frac{1}{2\pi(\tau_2-\tau_1^2)}} \exp\left(-\frac{(X-\tau_1)^2}{2(\tau_2-\tau_1^2)}\right) \]

它是均值为\(\tau_1\)，方差为\(\tau_2-\tau_1^2\)的正态分布。即：约束了随机变量期望、方差的分布为正态分布。事实上从以最大熵作为先验的分布都是指数族分布。

1.3 对数样条模型

考虑区间 \((L, U)\) 上的单变量密度估计，这里 \(L\) 与 \(U\) 可以为无穷。假设有 \(M \geq 4\) 个节点 \(L < t_1 < t_2 < \cdots < t_k < U\)。设函数 \(\{1, B_1(x), \cdots, B_P(x)\}\) 为 \(P+1\) 维空间上的一组基函数。

现在考虑用参数化模型 \(f_{(x|\theta)}(x)\) 对密度函数 \(f(x)\) 进行建模，即令：

\[ \log f_{(x|\theta)}(x | \theta) = \theta_1 B_1(x) + \cdots + \theta_P B_P(x) - C(\theta) \]

其中：

\[ \exp\{C(\theta)\} = \int_{L}^{U} \exp\{\theta_1 B_1(x) + \cdots + \theta_P B_P(x)\} dx \]

其中，\(\theta\) 表示所有列向量的集合，\(\theta = (\theta_1, \ldots, \theta_P)^T\)。由于参数化模型 \(f_{(x|\theta)}(x | \theta)\) 的形式是由密度函数的性质所决定的。所以这里要求 \(C(\theta)\) 是有限的。

在这一参数化的结构下，参数 \(\theta\) 的对数似然函数为：

\[ l(\theta | x_1, \ldots, x_n) = \sum_{i=1}^{n} \log f_{(x|\theta)}(x_i | \theta) \]

其中 \((x_1, \ldots, x_n)\) 为样本观测值。只要确定节点的位置，且每两个节点所构成的区间内包含用于估计的充分多的观测点，那么在 \(C(\theta)\) 有限的约束条件下最大化似然函数即可得到参数 \(\theta\) 的最大似然估计值。且因为 \(l(\theta | x_1, \ldots, x_n)\) 是凹函数，所以估计值是唯一的。在估计得到模型参数后，取：

\[ \hat{f}(x) = f_{(x|\theta)}(x | \hat{\theta}) \]

则 \(\hat{f}(x)\) 即为密度函数 \(f(x)\) 的最大似然对数样条密度估计。

2. 对数样条模型的似然估计

2.1 前情补充

牛顿迭代

设 \(x^*\) 是 \(f(x)=0\) 的根，选取 \(x_0\) 作为 \(x^*\) 的初始近似值。

过点 \((x_0, f(x_0))\) 处作曲线 \(y=f(x)\) 的切线 \(L_0\)：

\[ y = f(x_0) + f'(x_0) \cdot (x-x_0) \]

则 \(L_0\) 与 \(X\) 轴的交点横坐标为 \(x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}\)，则称 \(x_1\) 为 \(x^*\) 的一次近似值。

过点 \((x_1, f(x_1))\) 处作曲线 \(y=f(x)\) 的切线 \(L_1\)：

\[ y = f(x_1) + f'(x_1) \cdot (x-x_1) \]

则 \(L_1\) 与 \(X\) 轴的交点横坐标为 \(x_2 = x_1 - \frac{f(x_1)}{f'(x_1)}\)，则称 \(x_2\) 为 \(x^*\) 的二次近似值。

重复以上过程得到 \(x^*\) 的近似值序列。其中 \(x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}\) 称为 \(x^*\) 的 \(n+1\) 次近似值。

牛顿迭代公式：\(x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}\)

牛顿迭代示意图
牛顿迭代示意图

Hessian 矩阵

首先考虑一个一元函数求极值的问题。例如函数：\(f(x)=x^2\)。极值点处的一阶导数一定等于 0。但对于 \(f(x)=x^3\) 来说，只检查一阶导数就不能确定是否为极值了。因此则需要再求一次导，如果二阶导 \(f''(x)<0\)，那么说明函数在该点取得了局部极大值；如果二阶导数 \(f''(x)>0\)，则说明函数在该点取得了局部极小值；如果 \(f''(x)=0\)，则结果仍然是不确定的，就不得不通过其它方式来确定函数的极值性了。

而考虑在多元函数中求极值点，方法与此类似。以 \(f=f(x, y, z)\) 为例。首先，对于函数中的每个变量分别求偏导数，这是可知函数的极值点可能出现在哪里，即：

\[ \frac{\partial f}{\partial x}=0, \quad \frac{\partial f}{\partial y}=0, \quad \frac{\partial f}{\partial z}=0 \]

接下来，继续求二阶导数，此时包含混合偏导数的情况一共有 9 个，用矩阵形式来表示则有：

\[ H=\left[\begin{array}{lll} \frac{\partial^2 f}{\partial x \partial x} & \frac{\partial^2 f}{\partial x \partial y} & \frac{\partial^2 f}{\partial x \partial z} \\ \frac{\partial^2 f}{\partial y \partial x} & \frac{\partial^2 f}{\partial y \partial y} & \frac{\partial^2 f}{\partial y \partial z} \\ \frac{\partial^2 f}{\partial z \partial x} & \frac{\partial^2 f}{\partial z \partial y} & \frac{\partial^2 f}{\partial z \partial z} \end{array}\right] \]

也就是说, 设多元实函数 \(f\left(x_1, x_2, \cdots, x_n\right)\) 在点 \(M_0\left(a_1, a_2, \cdots, a_n\right)\) 的邻域内有二阶连续偏导, 若 有:

\[ \left.\frac{\partial f}{\partial x_j}\right|_{\left(a_1, a_2, \cdots, a_n\right)}=0, j=1,2, \cdots, n \]

并且:

\[ A=\left[\begin{array}{cccc} \frac{\partial^2 f}{\partial x_2^2} & \frac{\partial^2 f}{\partial x_1 \partial x_2} & \cdots & \frac{\partial^2 f}{\partial x_1 \partial x_n} \\ \frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x_2^2} & \cdots & \frac{\partial^2 f}{\partial x_2 \partial x_n} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^2 f}{\partial x_n \partial x_1} & \frac{\partial^2 f}{\partial x_n \partial x_2} & \cdots & \frac{\partial^2 f}{\partial x_n^2} \end{array}\right] \]

则可以得出: 1. 当 \(\mathrm{A}\) 为正定矩降时, \(f\left(x_1, x_2, \cdots, x_n\right)\) 在 \(M_0\left(a_1, a_2, \cdots, a_n\right)\) 处是极小值 2. 当 \(\mathrm{A}\) 为负定矩降时, \(f\left(x_1, x_2, \cdots, x_n\right)\) 在 \(M_0\left(a_1, a_2, \cdots, a_n\right)\) 处是极大值 3. 当 \(\mathrm{A}\) 为不定矩阵时, \(M_0\left(a_1, a_2, \cdots, a_n\right)\) 不是极值点 4. 当 \(\mathrm{A}\) 为半正定矩阵或半负定矩降时, \(M_0\left(a_1, a_2, \cdots, a_n\right)\) 是 “可疑 “的极值点, 仍需要其他 方法来判定。

2.2 对数样条模型的似然估计

设 \(H(\theta)\) 为 \(C(\theta)\) 的 Hessian 矩阵。其为 \(P \times P\) 矩阵，第 \((j,k)\) 个元素为：

\[ \frac{\partial^2 C(\theta)}{\partial \theta_j \partial \theta_k} = - \int_R B_j (x) B_k (x) f_{x| \theta} (x) dx + \int_R B_j (x) f_{x| \theta} (x) dx \int_R B_k (x) f_{x| \theta} (x) dx. \]

设 \(Y_1, \ldots, Y_n\) 是函数 \(f\) 中大小为 \(n\) 的样本。\(S(\theta)\) 为得分函数，是一个 \(p\) 维的向量，其第 \(j\) 个元素为：

\[ \frac{\partial l}{\partial \theta_j} (\theta) = \sum_{i=1}^n B_j (x_i ) - n \frac{\partial C}{\partial \theta_j} (\theta). \]

则充分统计量 \(b_j\) 为：

\[ b_j = \sum B_j (Y_i). \]

\(\hat{\theta}\) 的极大似然估计为：

\[ S(\hat{\theta}) = 0. \]

令 \(I(\theta) = -nH(\theta)\) 表示对应于随机样本的信息矩阵。并用牛顿迭代方法进行最大似然估计。即从一个起始的猜测值 \(\hat{\theta}^{(0)}\) 出发，通过迭代公式：

\[ \hat{\theta}^{(m+1)} = \hat{\theta}^{(m)} + I^{-1} (\hat{\theta}^{(m)}) S(\hat{\theta}^{(m)}), \]

进行迭代计第, 直到序列 \(\{\widehat{\theta}^{(m)}\}\) 的值稳定为止。

即当 \(l(\hat{\theta}^{(m+1)}) \leq l(\hat{\theta}^{(m)})\), 则应将 \(\hat{\theta}^{(m+1)}\) 替换为 \(\hat{\theta}^{(m+1)}-\hat{\theta}^{(m)}+\alpha I^{-1}(\hat{\theta}^{(m)}) S(\hat{\theta}^{(m)})\) 。其中常 数 \(\alpha \in(0,1)\), 通常选择为 \(\alpha=2^{-k}\), \(k\) 是使 \(l(\hat{\theta}^{(m)}+2^{-k} I^{-1}(\hat{\theta}^{(m)})-S(\hat{\theta}^{(m)}))>l(\hat{\theta}^{(m)})\) 成立 的最小正整数。

当 \(l(\hat{\theta}^{(m)})\) 越接近 \(\hat{\theta}\), 则选伐的收敛速度越快。通常当 \(\epsilon<10^{-6}\), 则认为是一个可行的收敛程度, 即:

\[ \epsilon=\Sigma \frac{\left(I^{-1}(\hat{\theta}^{(m)}) S(\hat{\theta}^{(m)})\right)_i^2}{\max \left((\hat{\theta}_i(m))^2, 0.00001\right)} \]

2.3 一阶二阶条件

因为

\[ f(y; \theta) = \exp \left( \theta_1 B_1 (y) + \ldots + \theta_p B_p (y) - C(\theta) \right), \quad L < y < U. \]

对于任意的 \(\theta \in \Theta\), \(f(; \theta)\) 是 \((L, U)\) 上的密度函数。

\[ F(y; \theta) = \int_L^y f(z; \theta) \, dz, \quad L < y < U, \]

\[ Q(p; \theta) = F^{-1} (p; \theta), \quad 0 < p < 1, \]

\(F(y; \theta)\) 和 \(Q(p; \theta)\) 是对应的分布函数与分位数函数。

对于给定的 \(y \in (L, U)\)，设

\[ \varphi(y; \theta) = \log(f(y; \theta)) = \sum_j \theta_j B_j (y) - C(\theta), \quad \theta \in \Theta. \]

给定一个长度为正的 \((L, U)\) 的子区间 \(A\)。设

\[ \varphi(A; \theta) = \log \left( \int_A f(y; \theta) \, dy \right) = \log \left( \int_A e^{\varphi(y; \theta)} \, dy \right), \quad \theta \in \Theta. \]

给定 \(\Theta\) 上的函数 \(g\)，设

\[ g_j (\theta) = \frac{\partial g(\theta)}{\partial \theta_j}, \quad g_{jk} (\theta) = \frac{\partial^2 g(\theta)}{\partial \theta_j \partial \theta_k}, \quad \theta \in \Theta, \]

其中 \(1 \leq j, k \leq p\)。那么就有

\[ \varphi_j (y; \theta) = B_j (y) - C_j (\theta) \quad \text{and} \quad \varphi_{jk} (y; \theta) = -C_{jk} (\theta), \quad \theta \in \Theta. \]

此时，当 \(A\) 为正时，有：

\[ \varphi_j (A; \theta) = \frac{\int_A \varphi_j (y; \theta) f(y; \theta) \, dy}{\int_A f(y; \theta) \, dy}, \quad \theta \in \Theta. \]

另外有：

\[ \begin{align*} \varphi_{jk} (A; \theta) &= \frac{\int_A \varphi_{jk} (y; \theta) f(y; \theta) \, dy}{\int_A f(y; \theta) \, dy} + \frac{\int_A \varphi_j (y; \theta) \varphi_k (y; \theta) f(y; \theta) \, dy}{\int_A f(y; \theta) \, dy} \\ &- \frac{\left( \int_A \varphi_j (y; \theta) f(y; \theta) \, dy \right) \left( \int_A \varphi_k (y; \theta) f(y; \theta) \, dy \right)}{\left( \int_A f(y; \theta) \, dy \right)^2}, \quad \theta \in \Theta. \end{align*} \]

3. 节点数量

选择节点数量的逻辑与选择带宽相似，节点过多可能会导致过度拟合并且波动较大，节点过少会导致过度平滑。这里介绍一个 Kooperberg & Stone节点选择方式。Kooperberg & Stone（1991）提出，节点的选择应该遵循以下原则：

  1. 节点的选择应该有一个简单且“自动”的规则
  2. 节点应该放置在统计量的位置，或统计量附近
  3. 节点应该在 (n + 1) / 2 左右近似对称分布
  4. 最低统计量和最最高统计量上应该有节点
  5. 极端的节点不应该与样本量有关
  6. 位于中间的节点应该基本等间隔

要满足上述六个条件，我们有以下规则：当 \((L, U) = (-∞, ∞)\) 时，取 \(t_1\)，\(t_K\) 分别为样本的最小和最大次序统计量的观测值，即 \(t_1 = x_{(1)}\)，\(t_K = x_{(n)}\)。当 \(1 < i < K/2\) 时，取 \(t_i = q(r_i)\)，即为样本观测值的 \(r_i\) 分位点。其中 \(0 < r_2 < r_3 < ⋯ < r_K < 1\)，其中 \(r_i\) 满足：

\[ n(r_{i+1} - r_i) = 4 \max⁡\{4 - \epsilon, 1\} \max⁡\{4 - 2\epsilon, 1\} \cdots \max⁡\{4 - (i - 1)\epsilon, 1\} \]

由此可以确定中位节点。接着，为了使剩余节点取得与之前节点对应的分位点对称，我们有：

\[ r_{K + 1 - i} - r_{K - i} = r_{i + 1} - r_i, \quad K / 2 ≤ i ≤ K - 1 \]

k 1 2 3 4 5 6 7
\(r_k\) 1 5 20 75 131 146 150
k 1 2 3 4 5 6 7 8 9 10
\(r_k\) 1 5 19 60 158 343 441 482 496 500

4. 平滑参数

4.1 加惩罚项

对于带惩罚项的 Log-spline 模型，正则化参数 \(\lambda\) 控制了拟合优度和估计密度函数平滑性之间的权衡。通过使用交叉验证或信息准则等方法来优化 \(\lambda\)，模型可以通过将不重要节点的系数收缩至零来有效地选择适当的节点数。一般地，在计算 Log-spline 的拟合误差的时候，我们使用：

\[ RSS(f, \lambda) = \sum_{i=1}^{N} (y_i - f(x_i))^2 \]

引入 \(\lambda\) 光滑因子式均方误差，我们有：

\[ RSS(f, \lambda) = \sum_{i=1}^{N} (y_i - f(x_i))^2 + \lambda \int (f''(t))^2 dt \]

当我们将 \(f(x)\) 带入到 \(RSS(f, \lambda)\) 中，\((f''(t))^2\) 分解，进行矩阵的乘积变换，我们可以得到：

\[ RSS(\lambda, \theta) = (y - N\theta)^T (y - N\theta) + \lambda \theta^T \Omega_N \theta \]

然后根据最小二乘思想，对 \(\theta\) 求导使得上式等于0，我们得到：

\[ \hat{\theta} = (N^T N + \lambda \Omega_N)^{-1} N^T y \]

将得到的 \(\hat{\theta}\) 带入原来的拟合函数可得：

\[ \hat{f}(x) = \sum_{j=1}^{N} N_j(x) \hat{\theta}_j \]

4.2 信息准则

BIC 通常会产生更简单的模型，而 AIC 则可能会产生更复杂的模型。

  • AIC：赤池信息量准则 (AIC) 是评估统计模型的复杂度和衡量统计模型“拟合”资料之优良性的一种标准，是由日本统计学家赤池弘次创立和发展的。赤池信息量准则建立在信息熵的概念基础上。在一般的情况下，AIC 可以表示为：

    \[ AIC = 2k - 2\ln(L) \]

    其中, \(k\) 是参数的数量, \(L\) 是似然函数。假设条件是模型的误差服从独立正态分布。让 \(n\) 为观察数，\(RSS\) 为残差平方和，那么 AIC 变为：

    \[ AIC = 2k + n \ln(RSS/n) \]

    增加自由参数的数目提高了拟合的优良性，AIC 鼓励数据拟合的优良性但尽量避免出现过度拟合的情况。所以优先考虑的模型应是 AIC 值最小的那一个。

  • BIC：贝叶斯信息准则 (BIC) 是主观贝叶斯派归纳理论的重要组成部分。公式为：

    \[ BIC = \ln(n)k - 2\ln(L) \]

    其中，\(k\) 为模型参数个数, \(n\) 为样本数量, \(L\) 为似然函数。在 Logspline 方法中，假设节点的选择方案共有 \(S\) 个，第 \(s\) 中方案的节点数为 \(m_s\)，\(s = 1, \cdots, S\)。则使：

    \[ BIC(s) = -2l(\|\hat{\theta}_s\| | x_1, \cdots, x_n) + (m_s - 1) \log n \]

    达到最小的节点选择方案为最终选择方案。在带惩罚项的对数样条密度估计中，

    \[ BIC_m = -2l(\hat{\theta}_m) + \text{dim}_m \log N/N, \quad m = 1, \cdots, M \]

AIC 和 BIC 的原理是不同的，AIC 是从预测角度，选择一个好的模型用来预测，BIC 是从拟合角度，选择一个对现有数据拟合最好的模型。AIC 和 BIC 的共性在于构造这些统计量所遵循的统计思想是一致的，就是在考虑拟合残差的同时，依自变量个数施加“惩罚”。AIC 和 BIC 的不同点在于，BIC 的惩罚项比 AIC 大，考虑了样本个数，样本数量多，可以防止模型精度过高造成的模型复杂度过高。当 \(n \geq 8\) 时，\(k \ln(n) \geq 2k\)，所以，BIC 相比 AIC 在大数据量时对模型参数惩罚得更多，导致 BIC 更倾向于选择参数少的简单模型。

5. 条件对数样条模型

回顾一下对数样条模型： \[ \begin{gathered} \log f_{x \mid \theta}(x \mid \theta)=\theta_1 B_1(x)+\ldots+\theta_P B_p(x)-C(\theta) \\ \exp \{C(\theta)\}=\int_L^U \exp \left\{\theta_1 B_1(x)+\ldots+\theta_P B_p(x)\right\} \mathrm{d} x \end{gathered} \]

假设引入一个协变承 \(\mathrm{X}\) 使新的密度滔数为 \(f(y \mid x) 。 f(y \mid x)\) 就是一个关于 \(\mathrm{Y}\) 在给定 \(X=x\) 条件 下的条件密度遥数。则它的条件对数样条模型可以写成:

\[ \begin{gathered} \log f(y \mid x ; \beta)=\sum_{k=0}^P \theta_k(x) B_k(y) \\ \theta_0(x)=-\log \int_{-\infty}^{\infty} \exp \left(\sum_{R=1}^P \theta_k(x) B_k(y)\right) \mathrm{d} y \end{gathered} \]

这里使用的仍然是 \(\mathrm{B}\) 样条, 归一化第数为 \(\theta_0(x)\) 。其中:

\[ \theta_k(x)=\sum_{j=0}^q \beta_{j k} H_j(x), k=1, \cdots, p \]

由于这里用的 \(\mathrm{B}\) 样条侈然是三次样条, 但是仅在节点 \(t_{k-1}, t_k, t_{k+1}, t_{k+2},(2 \leq k \leq p-2)\) 之间是三次样条。因此, \(B_k\) 在区间 \(\left(t_{k-1}, t_{k+3}\right)\) 为正, 其余均为 0 . \(\mathrm{B}\) 样条的表达式可写成:

\[ B_{p-1}(y)=\left(y-t_{K-3}\right)_{+}^3+a\left(y-t_{K-2}\right)_{+}^3+b\left(y-t_{K-2}\right)_{+}^3 c\left(y-t_K\right)_{+}^3 \]

其中 \(\mathrm{a}, \mathrm{b}, \mathrm{c}\), 使 \(B_{p-1}\) 的一二三阶导在 \(\left[t_K, \infty\right)\) 为 0 。

则:

\[ B_p(y)=\left(y-t_{K-2}\right)_{+}^3+d\left(y-t_{K-1}\right)_{+}^3+e\left(y-t_K\right)_{+}^3 \]

其中 \(\mathrm{d}\), 它, 使 \(B_p\) 的二三阶导在 \(\left[t_{K-2}, \infty\right)\) 为 \(0_0\) 由此, \(B_p\) 管要在 \(\left[t_{K-2}, \infty\right)\) 为正且递增, 并且在 \(\left[t_K, \infty\right)\) 为线性且科率为正。由此推导出来: \(B_1\) 在 \(\left(-\infty, t_3\right)\) 为正且递減, 并且在 \(\left(-\infty, t_1\right]\) 为线 性且科率为负。定义 \(B_0(y) \equiv 1\) 接者, 构站一个关于 \(\theta\) 的基遥数 \(H_0, \cdots, H_q, q=J-1\) 。则向桑 \(\beta\) 可以表达为:

\[ \int_{-\infty}^{\infty} \exp \left(\sum_{j, k=1}^p \beta_{j k} H_j(x) B_k(y)\right) \mathrm{d} y<\infty \text { 对一所有 } \mathrm{x} \]

且对于所有 \(x\) :

\[ \begin{aligned} & \theta_1(x)=\sum_j \beta_{j 1} H_j(x)<0 \\ & \theta_p(x)=\sum_j \beta_{j p} H_j(x)<0 \end{aligned} \]

因为 \(B_1, \ldots, B_P\) 具有 “compact supports”。 因此要满足上式，则需要添加条件：

\[ \beta_{01}<0,\cdots,\beta_{q1}<0,\beta_{0p}<0,\cdots,\beta_{qp}<0 \]

由此定义：

\[ f(y|x;\beta)=\exp\left(\sum_{k=1}^{P}\sum_{j=0}^{q}\beta_{jk}H_j(x)B_k(y)-C(x,\beta)\right) \]

其中：

\[ C(x,\beta)=log\int_{-\infty}^{\infty}exp\left(\sum_{k=1}^{F}\sum_{j=0}^{q}\beta_{jk}H_j(x)B_k(y)\right)\mathrm{d}y \]

则这里 \(f(y|x; \beta)\) 是一个条件对数样条模型。

补充： C. de Boor (1978). A Practical Guide to Splines. Springer-Verl

6. 总结和应用

6.1 总结

对数样条密度估计相对于其他密度估计技术（如核密度估计和基于直方图的方法）具有几个优点，包括：

  1. 灵活性：对数样条密度估计可以捕捉到潜在密度函数中的各种形状和模式，适用于各种应用。
  2. 平滑性：由于使用了B样条，对数样条估计本质上是平滑的，这可以提高估计的可解释性和可视化效果。
  3. 自动模型选择：确定最佳节点数和相应的B样条模型的过程内置于估计过程中，简化了模型选择。

6.2 经济学中的应用

对数样条模型在统计学中主要是用于非参数密度估计和分布函数估计。例如有人曾将对数样条密度估计用于肺癌生存期的统计。

Chen, Y. (2004). Logspline density estimation with an application to the study of survival data of lung cancer patients. East Tennessee State University.

在经济学中，对数样条模型可以用于估计价格的分布函数或者收入的密度函数，从而可以更准确地研究价格或者收入的分布特征。此外，对数样条模型还可以用于探索经济数据中的非线性关系，如非线性回归问题。例如，可以使用对数样条模型来建模收入对消费支出的非线性关系，从而更好地了解收入和消费支出之间的联系。

例子：

Matthee, M., & Naudé, W. (2007). Export diversity and regional growth: Empirical evidence from South Africa (No. 2007/11). WIDER Research Paper.

使用三次样条密度函数来估计次国家出口多样性与港口距离之间的关系，从而确定距离(可以代表运输成本)是否与出口多样性有显着关联。

6.3 收入分布估计的最大熵估计（例1）

  1. Ryu, H., 2008. Maximum entropy estimation of income distributions from bonferroni indices. In: Chotikapanich, Duangkamon (Ed.), Modeling Income Distributions and Lorenz Curves. Springer, pp. 193–211.
图6.1 Income Distribution & Gini Goefficient
图6.1 Income Distribution & Gini Goefficient

文章定义了Bonferroni Indices（BI），也就是在gini系数的基础上给更贫穷的人赋予更大的权重。在更加不平均的国家（例如巴西），BI和gini反应的结果可能差别很大，BI比gini更能反应真实情况。

  1. Ryu, H., 2013. A bottom poor sensitive gini coefficient and maximum entropy estimation of income distributions. Econom. Lett. 118, 370–374.
  1. Pgini（底层贫困敏感基尼系数）定义

过往的基尼系数直接使用收入作为度量对最贫困的人的考虑较少，本文使用生活中受到的不便作为度量。定义收入的倒数为受到不便的程度（收入越高，受到的不便越少），由此增加对最贫困的人所占的权重。

  • 定义pgini为：\(pg\text{ini}=\frac{1}{2n^2\mu}\sum_{i=1}^n\sum_{j=1}^n\left|\frac{1}{y_i}-\frac{1}{y_j}\right|\)

  • 社会存在的不便总量是\(\sum_{i=1}^n \frac{1}{y_i}\)

  • 平均不便是\(\mu=\frac{1}{n}\sum_{i=1}^n \frac{1}{y_i}\)

  • 设每个个体拥有的不便占\(r_i\)份额,\(r_i=\frac{1}{n\mu}\cdot\frac{1}{y_i}\)

则pgini可写成\(pg\text{ini}\equiv\frac{1}{2n}\sum_{i=1}^n\sum_{j=1}^n\left|r_i-r_j\right|\), 也就是下图中椭圆部分的面积。

图6.2 Inconvenience Lorenz Curve
图6.2 Inconvenience Lorenz Curve
  • Bonferroni index (收入不平等指数) \(B I=1-\int_0^1 \frac{L(z)}{z} \mathrm{~d} z\), BI 对穷人和富人有不同的权重, 穷人权重更大。其中 \(L(z)\) 是洛伦兹曲线。

对于个体 (离散) 的 \(\mathrm{BI}\) 可以写成: \(B I=\frac{1}{n-1} \sum_{i=1}^n \frac{P_i-Q_i}{P_i}\), 其中 \(P_i=\frac{i}{n}, Q_i=\sum_{j=1}^i \frac{x_j}{n \gamma}\) \(\gamma\) 是平均收入, \(x_j\) 是第 \(\mathrm{j}\) 个个体的收入。

  • 这里有一个标量极化指数, 用于比较贫困群体对收入份额变化的般感性: \(P=\frac{4 \gamma}{m} \mid 0.5-L(0,5)-0.5\) Gini \(\mid\), 其中 \(\mathrm{m}\) 是中位数收入,
  1. 用最大战估计求解
  • 洛伦兹曲线可被定义为: \(P L \equiv \int_0^z r\left(z^{\prime}\right) \mathrm{d} z^{\prime}\), 其中 \(r\left(z^{\prime}\right)\) 是关于不便份额的方程, \(z\) 是其中的人, \(z=0\) 表示最富有的人受到的不便最少, \(z=1\) 表示最贫穷的人有最多的不便。

则 \(P L \equiv \int_0^z r\left(z^{\prime}\right) \mathrm{d} z^{\prime}\) 的右边积分部分可以写成:

\[ \int_0^1 z \mathrm{~d} P L=z P L(z)_0^1-\int_0^1 P L(z) \mathrm{d} z=1-g \]

其中设 \(g=\int_0^1 P L(z) \mathrm{d} z=\frac{1-p g \text { ini }}{2}\)

因为: \(d P L(z)=r(z) d z\)

由此, 不便的均值方程可写成: \(\mu_1=\int_0^1 z r(z) \mathrm{d} z=1-g=\frac{1+p g i n i}{2}\)

  • 由此, 可以用最大化熵估计来求解不便份额 \(r(z)\) 的分布函数, 即: \(M A X_r W \equiv-\int r(z) \log r(z) d z\), 并满足 \(\int z r(z) \mathrm{d} z=\mu_1 \mid\)

  • 以上可用拉格期日方式求解可得:

\[ r(z) =\exp [a+b z]=\left[\frac{b}{e^b-1}\right] \exp [b z] \]

\[ \mu_1 =\left[\frac{b}{e^b-1}\right] \int_0^1 z \exp [b z] d z=\frac{1+\text { pgini }}{2} \]

\[ \mu_1\equiv-\frac{1}{b}+\frac{e^b}{e^b-1}=\frac{1+\text{pgini}}{2} \]

总而言之，在给定条件下最大化熵，从而推导出概率密度函数。由此得到关于不便分布的（也就是pgini）的函数。

  1. 实证

本文使用了2000-2009年美国家庭做样本，使用了五等分位的收入统计，由此可以算出五等分位的pgini系数。通过使用二次多项式对其进行光滑，产生“真实的”洛伦兹曲线和基尼系数。另外计算pgini系数曲线。发现对较贫困的（5%，10%）人群收入份额的变化高度相关，但对中产不那么相关。

Figure6.3 Gini&Pgini对比
Figure6.3 Gini&Pgini对比

通过最大化熵的方法，从pgini估算出来的收入分布对非常贫穷和非常富有的人的拟合效果更好。

Figure6.4 拟合结果
Figure6.4 拟合结果

6.4 收入分布的最大熵估计（例2）

  • Huang, Jianhua Z, et al. 2016, Estimation of a Probability Density Function Using Interval Aggregated Data. Journal of Statistical Computation and Simulation, vol. 86, no. 15

这篇文章中，作者用核密度估计建立起来了一种用来估计样本总体概率密度函数的方法，在这个方法建立的过程中，他们克服了核密度估计在此种问题下的准确度低调的缺点，因为他们关注到了数据聚合的过程。把样条方法和到核密度估计的应用问题结合起来，这个过程中产生了一些新的问题，比如说非标准的惩罚函数的计算问题，通过增加了数值计算的方法得以解决，以下是基于总体数据的核密度估计：

\[ \widehat{f}_{kae}(x)=\frac{1}{Ih}\sum\limits_{i=1}^I K(\frac{\widehat{m}_i-x}{h}) \]

在这篇文章中，在进行过数学模型和计算机模拟以后，作者尝试用这种新的模型解决了两个现实问题，分别是美国家庭收入数据分布和中国城镇家庭收入分布的估计。对于前者，惩罚样条估计的效果非常好，几乎可以做到密度曲线的重合，可以由此得出有用的信息。同时，在中国的数据中，作者发现了很有趣的点，作者通过作图发现26年之内的收入增长非常迅速，针对这一点，作者又使用了GIC来进行全体样本上的分布估计：

\[ GIC(u)=\left\{\frac{Q_{t+s}(u)}{Q_s(u)}\right\}^{\frac{1}{s}}-1,u\in[0,1] \]

作者从中发现了收入的不平等性，然后利用洛伦兹曲线，非常遗憾地发现，洛伦兹曲线在对角线下方，距离对角线越来越远甚至没有重叠，这代表收入不平等现象越发严重，作者们通过对不同问题的分析都得出了非常有价值的结论，特别是后一种，对于感兴趣的点进行深入讨论，最后可以得出意想不到的结果。

\[ L(u)=\frac{\int_0^u Q(u)du}{\int_0^1Q(u)du},u\in[0,1] \]

这种通过现实问题发现既有方法有所欠缺的情况有很多，这篇文章非常典范地展示了方法扩充之后灵活运用核密度估计分析金融问题的实例，这种对于方法扩充，然后追问，过程中使用不同工具，最后得出结论。

7. R中实现

安装logspline包

install.packages("logspline")
library(logspline)

7.1 dlogspline

从标准正态分布中生成一个100个值的随机样本。 使用logspline函数估计x中数值的密度。

x <- rnorm(100)
fit <- logspline(x)

使用qlogspline函数计算估计密度的量值。 使用plot函数创建一个估计密度与标准正态分布的图。

qq <- qlogspline((1:99)/100, fit)
plot(qnorm((1:99)/100, qq))
图7.1 R数据图
图7.1 R数据图

使用plogspline函数计算估计密度的累积分布函数。 使用plot函数创建一个估计的累积分布函数图。 在前面的的图中添加一条线，代表标准正态分布的累积分布函数。该图评估了估计的累积分布函数与标准正态累积分布函数的拟合程度。

pp <- plogspline((-250:250)/100, fit) 
plot((-250:250)/100, pp, type = "l") 
lines((-250:250)/100, pnorm((-250:250) /100))
图7.2 估计的分布函数与标准正态分布函数的拟合情况
图7.2 估计的分布函数与标准正态分布函数的拟合情况

使用 dlogspline 函数计算估计分布的密度 使用 plot 函数创建一个估计密度的图。 在前面的图中添加一条线，代表标准正态分布的密度。

dd <- dlogspline((-250:250)/100, fit) 
plot((-250:250)/100, dd, type = "l") 
lines((-250:250)/100, dnorm((-250:250)/100))
图7.3 估计的密度函数与标准正态密度的拟合情况
图7.3 估计的密度函数与标准正态密度的拟合情况

7.2 doldlogspline

qq <- qoldlogspline((1:99)/100, fit)
plot(qnorm((1:99)/100, qq))
图7.4 R数据图
图7.4 R数据图
pp <- poldlogspline((-250:250)/100, fit) 
plot((-250:250)/100, pp, type = "l") 
lines((-250:250)/100, pnorm((-250:250) /100))
图7.5 估计的分布函数与标准正态分布函数的拟合情况
图7.5 估计的分布函数与标准正态分布函数的拟合情况
dd <- doldlogspline((-250:250)/100, fit) 
plot((-250:250)/100, dd, type = "l") 
lines((-250:250)/100, dnorm((-250:250)/100))
图7.6 估计的密度函数与标准正态密度的拟合情况
图7.6 估计的密度函数与标准正态密度的拟合情况

7.3 logspline

y <- rnorm(100)
fit <- logspline(y)
plot(fit)
图7.7 logspline
图7.7 logspline

7.4 oldlogspline

y <- rnorm(100)
fit <- oldlogspline(y)
plot(fit)
图7.8 oldlogspline
图7.8 oldlogspline

7.5 其他计算方式

  • oldlogspline.to.logspline
  • plot.logspline
  • plot.oldlogspline
  • summary.logspline
  • summary.oldlogspline

参考文献

Ataspinar. “Regression, Logistic Regression and Maximum Entropy.” ML Fundamentals, 28 Mar. 2016, ataspinar.com/2016/03/28/regression-logistic-regression-and-maximum-entropy/. Accessed 13 June 2023.

C. Kooperberg & C. J. Stone (1991). A study of log-spline density estimation. Computational Statistics & Data Analysis, 12, 327

C. Kooperberg (1991). Smoothing images, curves, and densities. PhD Dissertation, Dept. of Statistics, University of California at Berkeley.

K. Bak, J. Jhong, J Lee, JShin, J Koo (2021). Penalized logspline density estimation using total variation penalty, Computational Statistics & Data Analysis, 153

K. Gu (1995). Smoothing spline density estimation: conditional distribution. Statistica S

R. Benoît. Mâsse, Young K. Truong (1999) Conditional Logspline Density Estimation, The Canadian Journal of Statistics , 27(4)