連続確率分布関数のコレクションに関連する問題があります。そのほとんどは経験的に決定されます（出発時間、通過時間など）。私が必要としているのは、これらのPDFのうちの2つを取得し、それらを計算する方法です。たとえば、PDFXから取得した2つの値xとPDFYから取得した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個のポイントが必要な場合は、各バーに10個のサンプルを割り当て、バーを可変幅で一定の面積にします。
PDFを追加するには、それらを畳み込むことがわかったので、そのための計算に骨を折った。2つの一様分布を畳み込むと、3つのセクションを持つ新しい分布が得られます。広い一様分布はまだ存在しますが、狭い方の幅の両側に三角形が貼り付けられています。したがって、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
他の演算子は、読者の練習問題として残されています。