我有一个涉及连续概率分布函数集合的问题,其中大部分是根据经验确定的(例如出发时间、过境时间)。我需要的是某种方法来获取其中两个 PDF 并对它们进行算术运算。例如。如果我有两个值 x 取自 PDF X,y 取自 PDF Y,我需要获取 (x+y) 的 PDF,或任何其他操作 f(x,y)。
分析解决方案是不可能的,因此我正在寻找允许此类操作的 PDF 的某种表示形式。一个明显(但计算成本昂贵)的解决方案是蒙特卡罗:生成大量 x 和 y 值,然后仅测量 f(x, y)。但这需要太多的CPU时间。
我确实考虑过将 PDF 表示为范围列表,其中每个范围具有大致相等的概率,从而有效地将 PDF 表示为均匀分布列表的并集。但我不知道如何将它们结合起来。
有人对这个问题有什么好的解决办法吗?
Edit:目标是创建一种用于操作 PDF 的迷你语言(也称为领域特定语言)。但首先我需要理清底层的表示和算法。
Edit 2:dmckee 建议使用直方图实现。这就是我通过统一分布列表得到的结果。但我不知道如何将它们组合起来创建新的发行版。最终,在 P(x
Edit 3:我有一堆直方图。它们分布不均匀,因为我是根据出现数据生成它们的,所以基本上,如果我有 100 个样本并且我想要直方图中有 10 个点,那么我会为每个条形分配 10 个样本,并使条形宽度可变但面积恒定。
我发现要添加 PDF,你需要对它们进行卷积,并且我已经为此加强了数学知识。当您对两个均匀分布进行卷积时,您会得到一个具有三个部分的新分布:较宽的均匀分布仍然存在,但每侧都有一个三角形,其宽度为较窄的分布。因此,如果我对 X 和 Y 的每个元素进行卷积,我将得到一堆这些元素,全部重叠。现在我试图弄清楚如何将它们全部相加,然后得到一个最接近它的直方图。
我开始怀疑蒙特卡洛到底是不是一个坏主意。
Edit 4: 这张纸 http://www.heldermann-verlag.de/eqc/eqc01_16/eqc16002.pdf详细讨论了均匀分布的卷积。一般来说,您会得到“梯形”分布。由于直方图中的每个“列”都是均匀分布的,因此我希望可以通过对这些列进行卷积并对结果求和来解决该问题。
然而,结果比输入复杂得多,并且还包括三角形。Edit 5:[错误内容已删除]。但是,如果这些梯形近似于具有相同面积的矩形,那么您就会得到正确的答案,并且减少结果中的矩形数量看起来也非常简单。这可能是我一直在寻找的解决方案。
Edit 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
其他运算符留给读者作为练习。