张春成
2023/03/04阅读：17主题：默认主题
精度消失
精度消失
本文尝试说明一个案例，该案例通过适当的数学约束，解决在概率密度函数计算过程中遇到的浮点数计算精度消失的问题。
-
精度消失[1] -
问题的定义[2] -
问题的解决[3] -
代码分享[4]
-
问题的定义
在前文中我们遇到了这样一个积分式，该用于计算二阶原点矩，如下式所示
其中用到的概率密度函数如下
其中， 代表标准正态分布的累积分布的概率密度函数， 代表采样次数， 代表全部 个样本的前一半元素的数量。如果您熟悉高等数学基本概念的话就会发现，这是一种典型的不定式问题，它类似于下式
我们不仅知道上式有界，还能知道它的值是确定的。但用计算机对它求解时会得到三个结果，它们分别是
-
分子 -
分母 ，该式用 float64 表示时为零； -
分母 ，该式用 float64 表示时为零。
也就是说，如果用计算机硬算的话，我们除了一堆 warning 和 0 之外，什么都得不到。也就是说，计算机的强大算力和精度无力在这个层面解决问题。
问题的解决
下面我将使用全概率公式的约束来尝试解决精度消失的问题，它非常简单
为了能够对它进行计算，首先将其进行离散化。对于计算机来说，这个水平的精度是可以保证的，因为离散值的粒度可以规定，且它们通常满足日常统计计算的精度要求。举个例子来说，这个精度通常代表集体中成员的数量或测量的毫米数等，因此 的规模可以在几百、几千或几万不等，对于计算机的算力来说绰绰有余。
由上式可知，我们不再需要清楚地知道每个 的值，而只需要知道它们之间的比例关系，就可以将它们求出来。而为了确定比例关系，我们可以绕开绝对值的计算，转而求它们的对数，即
因此，它们之间的比例关系可以用对数运算实现
在这个尺度上，计算机的浮点数完全可以胜任它的工作。
代码分享
以上的计算思路可以通过如下代码实现。由于概率密度函数的计算需要在不同尺度上进行，因此我们需要找到一个合适的转换函数将不同尺度上的概率值进行转换，使其可以在统一的尺度上进行比较和运算。本文给出的 convert 函数即实现了这一功能，它通过取对数然后进行归一化，实现了不同尺度上的概率值之间的转换。
由于开发过程很仓促，convert 函数写得像是 convert.min.js 的样子，没有什么可读性可言，希望大家能够理解。
def pesudo_probability(k, n):
'''
Compute the pesudo probability.
The probability is for the distribution of iid. n-sampled
from standard normalized distribution,
and the median value is k.
'''
m = (n - 1) / 2
f = stats.norm.cdf(k)
f1 = (1 - f)
return m * np.log(f * f1)
def convert_to_probability(pesudo_prob):
'''
Convert pesudo probability into probability
'''
a = np.array(pesudo_prob)
b = np.argmax(a)
c = a[b]
d = np.exp(a - c)
e = np.sum(d)
f = d / e
a, b, c, d, e, f
return f
expectation_median_variance = dict()
for n in tqdm(SAMPLES, 'Simulation'):
ks = np.linspace(-1, 1, 100)
pesudo_prob = [pesudo_probability(k, n) for k in ks]
prob = convert_to_probability(pesudo_prob)
var = np.sum((ks ** 2) * prob)
expectation_median_variance[n] = var
expectation_median_variance
参考资料
精度消失: #精度消失
[2]问题的定义: #问题的定义
[3]问题的解决: #问题的解决
[4]代码分享: #代码分享
作者介绍