Pau*_*son 22 algorithm math statistics haskell probability
我有一个涉及连续概率分布函数集合的问题,其中大多数是根据经验确定的(例如出发时间,运输时间).我需要的是一些方法来获取这些PDF中的两个并对它们进行算术运算.例如,如果我有两个值x取自PDF X,而y取自PDF Y,我需要得到(x + y)的PDF或任何其他操作f(x,y).
分析解决方案是不可能的,所以我正在寻找的是允许这些事情的PDF的一些表示.一个明显的(但计算上很昂贵的)解决方案是monte-carlo:生成大量的x和y值,然后测量f(x,y).但这需要太多的CPU时间.
我确实考虑将PDF表示为范围列表,其中每个范围具有大致相等的概率,有效地将PDF表示为统一分布列表的并集.但我看不出如何将它们结合起来.
有没有人对这个问题有任何好的解决方案？
编辑:目标是创建一种用于处理PDF的迷你语言(又称域特定语言).但首先我需要理清基础表示和算法.
编辑2: dmckee建议直方图实现.这就是我对统一分布列表的看法.但我不知道如何将它们组合起来创建新的发行版.最终我需要找到像P(x <y)这样的东西,这可能是非常小的.
编辑3:我有一堆直方图.它们不是均匀分布的,因为我是根据出现数据生成的,所以基本上如果我有100个样本并且我想在直方图中有10个点,那么我会为每个条形分配10个样本,并使条形变宽但是恒定区域.
我已经想出要添加PDF,你就可以对它们进行卷积,而且我已经为这个数学做了数学鉴定.当您对两个均匀分布进行卷积时,您会得到一个包含三个部分的新分布:更宽的均匀分布仍然存在,但是每边都有一个三角形,较窄的宽度.因此,如果我对X和Y的每个元素进行卷积,我会得到一堆这些,都是重叠的.现在我想弄清楚如何将它们全部相加,然后得到一个最佳近似直方图.
我开始怀疑蒙特卡洛毕竟不是一个坏主意.
编辑4: 本文详细讨论了均匀分布的卷积.一般来说,你得到一个"梯形"分布.由于直方图中的每个"列"是均匀分布,我曾希望通过卷积这些列并对结果求和来解决问题.
然而,结果比输入复杂得多,并且还包括三角形. 编辑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
其他操作员留给读者练习.
Ale*_* C. 16
不需要直方图或符号计算:如果采用正确的观点,一切都可以在封闭形式的语言层面完成.
[我将互换地使用术语"措施"和"分配".另外,我的Haskell生锈了,我请你原谅我在这方面不精确.]
概率分布实际上是codata.
设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
它也是一个单子(运动 - 暗示:这看起来像是延续单子.这是return什么？模拟的是call/cc什么？).
此外,结合微分几何框架,这可能会变成自动计算贝叶斯后验分布的东西.
在一天结束时,你可以写一些像
m = mean $ apply cos ((from_pdf gauss) ** (empirical data))
计算cos(X + Y)的平均值,其中X具有pdf gauss,Y已经由结果所在的MC方法采样data.
小智 6
概率分布形成一个单子; 例如,看看Claire Jones的作品以及LICS 1989年的论文,但这些想法可以追溯到Giry(DOI 10.1007/BFb0092872)的1982年论文和Lawvere的1962年的一篇文章,我无法追踪(http://永久链接. gmane.org/gmane.science.mathematics.categories/6541).
但是我没有看到comonad:没有办法从"(a-> Double) - > Double"获得"a".也许如果你把它变成多态 - (a-> r) - > r对所有r？(这是继续monad.)
有什么因素阻止您为此使用迷你语言吗？
我的意思是，定义一种语言，让您可以像编写时一样编写f = x + y和评估。f类似地g = x * z，h = y(x)等等。令人作呕。（我建议的语义要求评估者在评估时在 RHS 上出现的每个最里面的 PDF 上选择一个随机数，而不是尝试理解生成的 PDF 的合成形式。这可能不够快。 .)
假设您了解所需的精度限制，您可以使用直方图或样条曲线相当简单地表示 PDF（前者是后者的简并情况）。如果您需要将分析定义的 PDF 与实验确定的 PDF 混合在一起，则必须添加类型机制。
直方图只是一个数组，其内容表示输入范围的特定区域中的发生率。您没有说您是否有语言偏好，所以我会假设类似 c 的语言。您需要了解 bin 结构（uniorm 大小很容易，但并不总是最好），包括上限和下限以及可能的标准化：
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[];
};
这种事情在科学分析软件中很常见，您可能想使用现有的实现。