我有一个涉及连续概率分布函数集合的问题，其中大部分是根据经验确定的(例如出发时间、运输时间)。我需要的是获取其中两个 PDF 并对它们进行算术运算的某种方法。例如。如果我有两个取自 PDF X 的值 x 和取自 PDF Y 的 y，我需要获取 (x+y) 或任何其他操作 f(x,y) 的 PDF。
分析解决方案是不可能的，所以我正在寻找的是允许此类事情的 PDF 的某种表示形式。一个明显(但计算量大)的解决方案是蒙特卡洛法:生成大量 x 和 y 值，然后仅测量 f(x, y)。但这会占用太多 CPU 时间。
我确实考虑过将 PDF 表示为范围列表，其中每个范围的概率大致相等，从而有效地将 PDF 表示为均匀分布列表的并集。但我看不出如何组合它们。
有没有人有解决这个问题的好方法？
编辑:目标是创建一种用于操作 PDF 的微型语言(也称为领域特定语言)。但首先我需要梳理底层的表示和算法。
编辑 2: dmckee 建议使用直方图实现。这就是我的统一分布列表的目的。但我不知道如何将它们组合起来创建新的发行版。最终，我需要在 P(x < y) 可能非常小的情况下找到类似 P(x < y) 的东西。
编辑 3: 我有一堆直方图。它们不是均匀分布的，因为我是从出现数据生成它们的，所以基本上如果我有 100 个样本并且我想在直方图中有 10 个点，那么我会为每个条分配 10 个样本，并使条的宽度可变但面积不变。
我发现要添加 PDF 需要对它们进行卷积，并且我已经为此加深了数学知识。当你对两个均匀分布进行卷积时，你会得到一个包含三个部分的新分布:较宽的均匀分布仍然存在，但在较窄的均匀分布的每一侧都有一个三角形。因此，如果我对 X 和 Y 的每个元素进行卷积，我将得到一堆这些元素，它们全部重叠。现在我想弄清楚如何将它们全部相加，然后得到一个与其最接近的直方图。
我开始怀疑蒙特卡洛到底是不是个坏主意。
编辑 4: This paper详细讨论了均匀分布的卷积。一般来说，你会得到一个“梯形”分布。由于直方图中的每个“列”都是均匀分布的，我曾希望通过对这些列进行卷积并对结果求和来解决问题。
然而，结果比输入复杂得多，并且还包括三角形。 编辑 5: [删除了错误的内容]。但是，如果将这些梯形近似为具有相同面积的矩形，那么您就会得到正确答案，并且减少结果中的矩形数量看起来也非常简单。这可能是我一直在努力寻找的解决方案。
编辑 6: 已解决!这是此问题的最终 Haskell 代码:
-- | 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
其他算子留给读者练习。
最佳答案
不需要直方图或符号计算:如果采取正确的观点，一切都可以在语言级别以封闭形式完成。
[我将交替使用术语“度量”和“分发”。另外，我的 Haskell 已经生锈了，我请你原谅我在这方面的不精确。]
概率分布实际上是余数据。
令 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)
因此，给定两个度量，您可以针对它们的卷积进行集成。
此外，给定法则 mu 的随机变量 X，a * X 的法则由下式给出:
rescale :: Double -> Measure -> Measure
rescale a mu f = mu $ \x -> f(a * x)
此外，phi(X) 的分布由 image measure 给出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
它也是一个 monad(练习——提示:这很像 continuation monad。什么是 return？什么是 call/cc 的类比？ ).
此外，结合微分几何框架，这可能会变成自动计算贝叶斯后验分布的东西。
在一天结束的时候，你可以写这样的东西
m = mean $ apply cos ((from_pdf gauss) ** (empirical data))
计算 cos(X + Y) 的平均值，其中 X 具有 pdf gauss 并且 Y 已通过 MC 方法采样，其结果在 data 中。
关于algorithm - 表示连续概率分布，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/395981/