表示连续概率分布

22

我有一个问题,涉及到一组连续概率分布函数,其中大部分是经验确定的(例如:出发时间、运输时间)。我需要的是找到一种方法,能够对两个这样的PDF进行算术计算。例如,如果我从PDF X中取得两个值x和y,那么我需要得到(x+y)或任何其他操作f(x,y)的PDF。

由于不可能存在解析解,因此我正在寻找一些能够处理这些操作的PDF表示方式。一个显而易见(但计算复杂度较高)的解决方案是通过蒙特卡罗方法:生成大量的x和y值,然后测量f(x, y)。但这需要太多的CPU时间。

我曾想过将PDF表示为一系列范围的列表,其中每个范围具有大致相等的概率,以有效地将PDF表示为一系列均匀分布的并集。但我无法看到如何将它们结合起来。

是否有人有解决这个问题的好方法?

编辑:目标是创建一个用于操作PDF的小型语言(也称为特定领域语言)。但首先我需要解决底层表示和算法。

编辑2:dmckee提出了一种直方图实现方法。这就是我通过均匀分布列表所得到的结果。但我看不出如何将它们结合起来以创建新的分布。最终,我需要找到像P(x < y)这样的东西,在这种情况下可能非常小。

编辑3:我有一堆直方图。它们不是均匀分布的,因为我从发生数据中生成它们。因此,如果我有100个样本,并且我想在直方图中有十个点,则我会将每个条分配10个样本,并使条形变量宽度但面积恒定。

我已经弄清楚了如何添加PDF,即对它们进行卷积,我已经学习了相关数学知识。当您对两个均匀分布进行卷积时,您将得到一个新的分布,其中包含三个部分:较宽的均匀分布仍然存在,但在较窄的分布两侧各有一个三角形。因此,如果我将X和Y的每个元素进行卷积,我将得到许多这样的重叠部分。现在我正在努力弄清楚如何对它们求和,然后得到最佳近似的直方图。

我开始怀疑蒙特卡罗并不是一个坏主意。

编辑4: 这篇论文详细讨论了均匀分布的卷积问题。总的来说,你会得到一个“梯形”分布。由于直方图中的每个“柱子”都是均匀分布,因此我希望这个问题可以通过对这些列进行卷积并对结果求和来解决。

然而,结果比输入复杂得多,并且还包括三角形。但是,如果将这些梯形近似为具有相同面积的矩形,则会得到正确的答案,并且减少结果中的矩形数量

-- | Continuous distributions of scalars are represented as a
-- | histogram where each bar has approximately constant area but
-- | variable width and height.  A histogram with N bars is stored as
-- | a list of N+1 values.
data Continuous = C {
      cN :: Int,
      -- ^ Number of bars in the histogram.
      cAreas :: [Double],
      -- ^ Areas of the bars.  @length cAreas == cN@
      cBars :: [Double]
      -- ^ Boundaries of the bars.  @length cBars == cN + 1@
    } deriving (Show, Read)


{- | Add distributions.  If two random variables @vX@ and @vY@ are
taken from distributions @x@ and @y@ respectively then the
distribution of @(vX + vY)@ will be @(x .+. y).

This is implemented as the convolution of distributions x and y.
Each is a histogram, which is to say the sum of a collection of
uniform distributions (the "bars").  Therefore the convolution can be
computed as the sum of the convolutions of the cross product of the
components of x and y.

When you convolve two uniform distributions of unequal size you get a
trapezoidal distribution. Let p = p2-p1, q - q2-q1.  Then we get:


>   |                              |
>   |     ______                   |
>   |     |    |           with    |  _____________
>   |     |    |                   |  |           |
>   +-----+----+-------            +--+-----------+-
>         p1   p2                     q1          q2
> 
>  gives    h|....... _______________
>            |       /:             :\
>            |      / :             : \                1
>            |     /  :             :  \     where h = -
>            |    /   :             :   \              q
>            |   /    :             :    \
>            +--+-----+-------------+-----+-----
>             p1+q1  p2+q1       p1+q2   p2+q2

However we cannot keep the trapezoid in the final result because our
representation is restricted to uniform distributions.  So instead we
store a uniform approximation to the trapezoid with the same area:

>           h|......___________________
>            |     | /               \ |
>            |     |/                 \|
>            |     |                   |
>            |    /|                   |\
>            |   / |                   | \
>            +-----+-------------------+--------
>               p1+q1+p/2          p2+q2-p/2

-}
(.+.) :: Continuous -> Continuous -> Continuous
c .+. d = C {cN     = length bars - 1,
             cBars  = map fst bars, 
             cAreas = zipWith barArea bars (tail bars)}
    where
      -- The convolve function returns a list of two (x, deltaY) pairs.
      -- These can be sorted by x and then sequentially summed to get
      -- the new histogram.  The "b" parameter is the product of the
      -- height of the input bars, which was omitted from the diagrams
      -- above.
      convolve b c1 c2 d1 d2 =
          if (c2-c1) < (d2-d1) then convolve1 b c1 c2 d1 d2 else convolve1 b d1 
d2 c1 c2
      convolve1 b p1 p2 q1 q2 = 
          [(p1+q1+halfP, h), (p2+q2-halfP, (-h))]
               where 
                 halfP = (p2-p1)/2
                 h = b / (q2-q1)
      outline = map sumGroup $ groupBy ((==) `on` fst) $ sortBy (comparing fst) 
$ concat
                [convolve (areaC*areaD) c1 c2 d1 d2 |
                 (c1, c2, areaC) <- zip3 (cBars c) (tail $ cBars c) (cAreas c),
                 (d1, d2, areaD) <- zip3 (cBars d) (tail $ cBars d) (cAreas d)
                ]
      sumGroup pairs = (fst $ head pairs, sum $ map snd pairs)

      bars = tail $ scanl (\(_,y) (x2,dy) -> (x2, y+dy)) (0, 0) outline
      barArea (x1, h) (x2, _) = (x2 - x1) * h

其他运算符留给读者作为练习。

10个回答

16
不需要直方图或符号计算:如果采用正确的观点,一切都可以在语言层面上以封闭形式完成。
“测度”和“分布”这两个术语可以互换使用。请注意,测度只能与测试函数进行积分(这是“测度”的一个可能的数学定义)。最终你将会做的事情也就是这个:例如,对恒等式积分就是求平均值。
概率分布实际上是协程数据。
设mu为概率测度。你可以对测度执行的唯一操作就是将其与测试函数进行积分。
mean :: Measure -> Double
mean mu = mu id

另一个例子:

variance :: Measure -> Double
variance mu = (mu $ \x -> x ^ 2) - (mean mu) ^ 2

另一个例子,计算P(mu < x)的方法如下:
cdf :: Measure -> Double -> Double
cdf mu x = mu $ \z -> if z < x then 1 else 0

这表明了一种通过二元性的方法。
因此,类型“Measure”应表示类型“(Double -> Double) -> Double”。 这使您可以对MC模拟结果、针对PDF的数值/符号积分等进行建模。例如,函数
empirical :: [Double] -> Measure
empirical h:t f = (f h) + empirical t f

返回对 f 进行积分的结果,该积分针对通过 MC 抽样获得的 经验测度。此外

from_pdf :: (Double -> Double) -> Measure
from_pdf rho f = my_favorite_quadrature_method rho f

从(正常)密度构建测量。

现在,好消息是,如果mu和nu是两个测量,则卷积 mu ** nu 如下所示:

(mu ** nu) f = nu $ \y -> (mu $ \x -> f $ x + y)

因此,给定两个度量,您可以对它们的卷积进行积分。

另外,给定一个随机变量X,其律为mu,a * X的律为:

rescale :: Double -> Measure -> Measure
rescale a mu f = mu $ \x -> f(a * x)

此外,在我们的框架中,phi(X)的分布由图像测度phi_* X给出:(推进测量)
apply :: (Double -> Double) -> Measure -> Measure
apply phi mu f = mu $ f . phi

现在,您可以轻松地为度量制定嵌入式语言。这里还有很多事情要做,特别是与实线以外的样本空间、随机变量之间的依赖关系、条件等方面有关,但我希望您能理解重点。

特别地,推进是函子性的:

newtype Measure a = (a -> Double) -> Double
instance Functor Measure a where
    fmap f mu = apply f mu

它也是一个单子(练习提示:它非常类似于 continuation monad。什么是 return?什么是 call/cc 的类比?)。

此外,与微分几何框架结合起来,这可能可以变成一些自动计算贝叶斯后验分布的东西。

最终,你可以编写像这样的内容:

m = mean $ apply cos ((from_pdf gauss) ** (empirical data))

计算cos(X + Y)的平均值,其中X的概率密度函数为高斯分布,Y通过MC方法进行采样,其结果存储在data中。


这非常有趣。 - jtobin
我基于这个答案编写了一个小的概念验证库,供您参考。您可以在此处找到一个示例:http://github.com/jtobin/measurable/blob/master/tests/Test.hs。 - jtobin
@jtobin:太好了,我会看一下的! - Alexandre C.
2
一个特别酷的事情值得注意的是,(>>=) 贝叶斯推断。即priorMeasure >>= likelihoodMeasure产生了后验预测度量。 - jtobin
2
以防有人读了我的上面的评论并感到困惑:我也过于草率地得出了那个结论。Bind实际上是一个整合运算符,因为priorMeasure >>= likelihoodMeasure产生了先验预测分布。当您传入后验时,后验预测自然会跟随,但>>=不会自动处理后验。 - jtobin
显示剩余2条评论

6
概率分布形成一个单子范畴;参见Claire Jones的著作以及LICS 1989论文,但这些思想可以追溯到Giry于1982年发表的文章(DOI 10.1007/BFb0092872)和Lawvere于1962年发表的一篇笔记,我无法追踪其来源(链接:http://permalink.gmane.org/gmane.science.mathematics.categories/6541)。
但我没有看到余单子:没有办法从“(a->Double)->Double”中取出“a”。也许如果将其变成多态的 - 对于所有的r,(a->r)->r?(那就是延续单子。)

2
我曾经在我的论文中研究过类似的问题。
计算近似卷积的一种方法是对密度函数(在这种情况下为直方图)进行傅里叶变换,将它们相乘,然后进行逆傅里叶变换以得到卷积。
请参阅我的论文附录C中各种特殊情况概率分布操作的公式。你可以在以下网址找到论文:http://riso.sourceforge.net 我编写了Java代码来执行这些操作。你可以在以下网址找到代码:https://sourceforge.net/projects/riso

2

有什么阻止你使用小语言来完成这个任务吗?

我的意思是,定义一种语言,让你可以写出 f = x + y 并像写出的那样计算 f。同样的,对于 g = x * z, h = y(x) 等等,都可以这样做。 (我建议的语义是,在评估时,计算器在每个最内层 PDF 上选择一个随机数,并不尝试理解所得到的 PDF 的组合形式。这可能不够快...)


假设您了解所需的精度限制,您可以用直方图或样条函数(前者是后者的退化形式)相当简单地表示 PDF。如果您需要将分析定义的 PDF 与实验确定的 PDF 混合在一起,则必须添加类型机制。


直方图只是一个数组,其内容表示输入范围中特定区域的发生率。您没有说您是否有语言偏好,因此我们假设您使用类 C 语言。您需要知道 bin 结构(统一大小很容易,但并不总是最佳),包括高低限制和可能的标准化:

struct histogram_struct {
  int bins; /* Assumed to be uniform */
  double low;
  double high;
  /* double normalization; */    
  /* double *errors; */ /* if using, intialize with enough space, 
                         * and store _squared_ errors
                         */
  double contents[];
};

这种事情在科学分析软件中非常常见,您可能希望使用现有的实现。


我确实想编写一种迷你语言。但是我希望底层语义比蒙特卡罗更高效。 - Paul Johnson
抱歉,我表达不够清晰。您真的需要知道复合PDF是什么,还是只需要能够高效地从中提取数字?后一种情况实际上只需要在代码内部使用一个良好的PDF组合解释器即可。 - dmckee --- ex-moderator kitten
我想知道什么是复合PDF。最终,我需要能够确定像P(x<y)这样的事情。 - Paul Johnson
你可以用直方图相对简单地表示PDF。是的。你知道有没有相关算法吗? - Paul Johnson
四处寻找保存这个想法的方法... 呃,嗯,嗯。怎么样,建立一个集群。是的,你需要一个集群... - dmckee --- ex-moderator kitten

1

1

几个回复:

1)如果您已经确定了经验PDF,则要么您有直方图,要么您有参数PDF的近似值。 PDF是连续函数,您没有无限的数据...

2)假设变量是独立的。然后,如果将PDF离散化,则P(f(x,y))= f(x,y)p(x,y)= f(x,y)p(x)p(y),总和为所有组合的x和y,使得f(x,y)满足您的目标。

如果您要将经验PDF拟合为标准PDF,例如正态分布,则可以使用已确定的函数来计算总和等。

如果变量不独立,则您会遇到更多麻烦,我认为您必须使用copulas

我认为定义自己的小语言等是过度的。 您可以使用数组完成此操作...


1
一些初步想法:
首先,Mathematica具有使用精确分布进行计算的良好工具。
其次,用直方图(即经验PDF)表示存在问题,因为您必须对箱子大小进行选择。这可以通过存储累积分布来避免,即经验CDF。(事实上,您还保留了重新创建基于经验分布的完整样本数据集的能力。)
以下是一些丑陋的Mathematica代码,用于将样本列表转换为经验CDF,即值-概率对的列表。将此输出通过ListPlot运行以查看经验CDF的绘图。
最后,这里有一些关于组合离散概率分布的信息:

http://www.dartmouth.edu/~chance/teaching_aids/books_articles/probability_book/Chapter7.pdf


我同意你关于箱子大小的观点,但历史数据的量很大,所以我会坚持使用箱子。不过,了解经验累积分布函数的内容还是很有用的。我已经阅读了你指出的那一章。谢谢。 - Paul Johnson

1

我认为直方图或1/N区域列表是个好主意。为了论证,我假设所有分布都有一个固定的N。

使用您链接的第四版编辑的论文来生成新的分布。然后,用一个新的N元素分布来近似它。

如果您不想让N固定,那就更容易了。将新生成的分布中的每个凸多边形(梯形或三角形)用均匀分布来近似。


1
另一个建议是使用核密度。特别是如果您使用高斯核,那么它们可以相对容易地处理...除非在不小心的情况下,分布会迅速增大。根据应用程序,还有其他近似技术,如重要性抽样可供使用。

0

如果你想玩一下,可以尝试像Maple或Mathematica那样用符号表示它们。Maple使用有向无环图,而Matematica使用类似于列表/ lisp的方法(我相信,但这已经是很久以前的事情了,自从我甚至想过这个问题)。

在符号上进行所有操作,然后最终推导出数值。 (或者只需找到一种在shell中启动并执行计算的方法)。

保罗。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接