algorithm - uniforme - Representando distribuciones de probabilidad continuas

tipos de distribuciones de probabilidad

Question

Tengo un problema que involucra una colección de funciones de distribución de probabilidad continua, la mayoría de las cuales se determinan empíricamente (por ejemplo, tiempos de salida, tiempos de tránsito). Lo que necesito es alguna forma de tomar dos de estos archivos PDF y hacer aritmética sobre ellos. Por ejemplo, si tengo dos valores x tomados de PDF X ey tomados de PDF Y, necesito obtener el PDF para (x + y) o cualquier otra operación f (x, y).

Una solución analítica no es posible, entonces lo que estoy buscando es una representación de archivos PDF que permita tales cosas. Una solución obvia (pero computacionalmente costosa) es montecarlo: genere muchos valores de xey, y luego simplemente mida f (x, y). Pero eso lleva demasiado tiempo de CPU.

Pensé en representar el PDF como una lista de rangos donde cada rango tiene una probabilidad aproximadamente igual, representando efectivamente el PDF como la unión de una lista de distribuciones uniformes. Pero no puedo ver cómo combinarlos.

¿Alguien tiene alguna buena solución para este problema?

Editar: El objetivo es crear un mini-lenguaje (también conocido como Dominio Específico del Idioma) para manipular archivos PDF. Pero primero necesito ordenar la representación y los algoritmos subyacentes.

Editar 2: dmckee sugiere una implementación de histograma. A eso es a lo que me refería con mi lista de distribuciones uniformes. Pero no veo cómo combinarlos para crear nuevas distribuciones. En última instancia, necesito encontrar cosas como P (x <y) en los casos en que esto puede ser bastante pequeño.

Edición 3: Tengo un montón de histogramas. No están distribuidos uniformemente porque los estoy generando a partir de datos de ocurrencia, así que básicamente si tengo 100 muestras y quiero diez puntos en el histograma, entonces asigno 10 muestras a cada barra y hago que las barras tengan ancho variable pero área constante.

Me di cuenta de que para agregar archivos PDF los convivía, y me he recuperado de las matemáticas. Cuando convinas dos distribuciones uniformes, obtienes una nueva distribución con tres secciones: la distribución uniforme más amplia todavía está allí, pero con un triángulo pegado a cada lado del ancho del más estrecho. Entonces, si convierto cada elemento de X e Y obtendré un montón de estos, todos superpuestos. Ahora estoy tratando de encontrar la forma de sumarlos todos y luego obtener un histograma que sea la mejor aproximación.

Me estoy empezando a preguntar si Montecarlo no era tan mala idea después de todo.

Edición 4: Este documento discute convoluciones de distribuciones uniformes con cierto detalle. En general, obtienes una distribución "trapezoidal". Dado que cada "columna" en los histogramas es una distribución uniforme, esperaba que el problema pudiera resolverse mediante la convolución de estas columnas y la suma de los resultados.

Sin embargo, el resultado es considerablemente más complejo que las entradas, y también incluye triángulos. Editar 5: [Se borraron cosas incorrectas]. Pero si estos trapecios se aproximan a rectángulos con la misma área, entonces obtienes la respuesta correcta, y la reducción del número de rectángulos en el resultado parece bastante directa también. Esta podría ser la solución que he estado tratando de encontrar.

Editar 6: ¡Resuelto! Aquí está el código final de Haskell para este problema:

-- | 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

Otros operadores se dejan como un ejercicio para el lector.

Answer 1

¿Hay algo que te impida emplear un mini-lenguaje para esto?

Con eso quiero decir, defina un lenguaje que le permita escribir f = x + y y evalúe f por usted tal como está escrito. Y de manera similar para g = x * z , h = y(x) , etc. ad nauseum . (La semántica sugiero que se solicite al evaluador que seleccione un número aleatorio en cada PDF interno que aparezca en el RHS en el momento de la evaluación, y que no intente comprender la forma compostada de los PDF resultantes. Puede que esto no sea lo suficientemente rápido ... .)

Suponiendo que comprende los límites de precisión que necesita, puede representar un PDF de forma bastante simple con un histograma o spline (el primero es un caso degenerado del último). Si necesita mezclar PDF definidos analíticamente con los determinados experimentalmente, tendrá que agregar un mecanismo de tipo.

Un histograma es solo una matriz, cuyo contenido representa la incidencia en una región particular del rango de entrada. No has dicho si tienes una preferencia de idioma, así que asumiré algo parecido a un c. Necesita conocer la estructura bin (los tamaños de las funciones son fáciles pero no siempre los mejores), incluidos los límites alto y bajo y posiblemente la normalización:

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[];
};

Este tipo de cosas es muy común en el software de análisis científico, y es posible que desee utilizar una implementación existente.

Answer 2

La robótica móvil autónoma trata un problema similar en localización y navegación, en particular la localización de Markov y el filtro de Kalman (fusión del sensor). Ver Una comparación experimental de métodos de localización continuó, por ejemplo.

Otro enfoque que podría tomar prestado de los robots móviles es la planificación de rutas usando campos potenciales.

Answer 3

Si quieres un poco de diversión, intenta representarlos simbólicamente como lo harían Maple o Mathemetica. Maple usa gráficos acíclicos dirigidos, mientras que Matematica usa una lista / ceceo como appoach (creo, pero ha sido mucho tiempo, ya que incluso pensé en esto).

Haz todas tus manipulaciones simbólicamente, luego al final empuja a través de valores numéricos. (O simplemente encuentre una forma de iniciarse en un shell y hacer los cálculos).

Pablo.

Answer 4

Un par de respuestas:

1) Si ha determinado PDF de manera empírica, o tiene histogramas o tiene una aproximación a un PDF paramétrico. Un PDF es una función continua y no tienes datos infinitos ...

2) Supongamos que las variables son independientes. Luego, si haces que el PDF sea discreto, entonces P (f (x, y)) = f (x, y) p (x, y) = f (x, y) p (x) p (y) sumado a todas las combinaciones de x y y tal que f (x, y) se encuentra con tu objetivo.

Si va a ajustar los archivos PDF empíricos a archivos PDF estándar, por ejemplo, la distribución normal, puede usar las funciones ya determinadas para calcular la suma, etc.

Si las variables no son independientes, entonces tiene más problemas en sus manos y creo que debe usar cópulas .

Creo que definir tu propio mini lenguaje, etc., es excesivo. puedes hacer esto con arreglos ...

Answer 5

Algunos pensamientos iniciales:

Primero, Mathematica tiene una buena facilidad para hacer esto con distribuciones exactas.

En segundo lugar, la representación como histogramas (es decir, archivos PDF empíricos) es problemática ya que debe tomar decisiones sobre el tamaño del contenedor. Eso se puede evitar almacenando una distribución acumulativa, es decir, una CDF empírica. (De hecho, conserva la capacidad de recrear el conjunto de datos completo de las muestras en las que se basa la distribución empírica).

Aquí hay un feo código de Mathematica para tomar una lista de muestras y devolver una CDF empírica, es decir, una lista de pares de valor-probabilidad. Ejecute el resultado de esto a través de ListPlot para ver un diagrama del CDF empírico.

empiricalCDF[t_] := Flatten[{{#[[2,1]],#[[1,2]]},#[[2]]}&/@Partition[Prepend[Transpose[{#[[1]], Rest[FoldList[Plus,0,#[[2]]]]/Length[t]}&[Transpose[{First[#],Length[#]}&/@ Split[Sort[t]]]]],{Null,0}],2,1],1]

Finalmente, aquí hay información sobre cómo combinar distribuciones de probabilidad discretas:

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

Answer 6

Creo que los histogramas o la lista de regiones de área 1 / N es una buena idea. Por el bien del argumento, supondré que tendrá una N fija para todas las distribuciones.

Utilice el papel que ha vinculado edit 4 para generar la nueva distribución. Luego, aproximémoslo con una nueva distribución de elementos N.

Si no quiere que se arregle N, es aún más fácil. Tome cada polígono convexo (trapecio o triángulo) en la nueva distribución generada y aproxiínelo con una distribución uniforme.

Answer 7

Otra sugerencia es usar densidades de kernel . Especialmente si usa núcleos gaussianos, entonces puede ser relativamente fácil trabajar con ellos ... excepto que las distribuciones explotan rápidamente en tamaño sin cuidado. Dependiendo de la aplicación, existen técnicas de aproximación adicionales como el muestreo de importancia que se puede usar.

Answer 8

Trabajé en problemas similares para mi disertación.

Una forma de calcular circunvoluciones aproximadas es tomar la transformada de Fourier de las funciones de densidad (histogramas en este caso), multiplicarlas, luego tomar la transformada de Fourier inversa para obtener la convolución.

Mire el Apéndice C de mi disertación para las fórmulas de varios casos especiales de operaciones sobre distribuciones de probabilidad. Puede encontrar la disertación en: http://riso.sourceforge.net

Escribí código Java para llevar a cabo esas operaciones. Puede encontrar el código en: https://sourceforge.net/projects/riso

Answer 9

No hay necesidad de histogramas o cálculos simbólicos: todo se puede hacer en el nivel del idioma en forma cerrada, si se toma el punto de vista correcto.

[Utilizaré el término "medida" y "distribución" indistintamente. Además, mi Haskell está oxidado y te pido que me perdones por ser impreciso en esta área.]

Las distribuciones de probabilidad son realmente codata .

Deje que mu sea una medida de probabilidad. Lo único que puede hacer con una medida es integrarla contra una función de prueba (esta es una posible definición matemática de "medida"). Tenga en cuenta que esto es lo que eventualmente hará: por ejemplo, integrarse en contra de la identidad es tomar la media:

mean :: Measure -> Double
mean mu = mu id

otro ejemplo:

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

otro ejemplo, que calcula P (mu <x):

cdf :: Measure -> Double -> Double
cdf mu x = mu $ /z -> if z < x then 1 else 0

Esto sugiere un acercamiento por dualidad.

El tipo Measure indicará, por lo tanto, el tipo (Double -> Double) -> Double . Esto le permite modelar resultados de simulación MC, cuadratura numérica / simbólica contra un PDF, etc. Por ejemplo, la función

empirical :: [Double] -> Measure
empirical h:t f = (f h) + empirical t f

devuelve la integral de f contra una medida empírica obtenida por ej. Muestreo MC también

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

construir medidas a partir de densidades (regulares).

Ahora, las buenas noticias. Si mu y nu son dos medidas, la convolución mu ** nu viene dada por:

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

Entonces, dadas dos medidas, puedes integrarte contra su convolución.

Además, dada una variable aleatoria X de ley mu , la ley de a * X viene dada por:

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

Además, la distribución de phi (X) viene dada por la medida de imagen phi_ * X, en nuestro marco:

apply :: (Double -> Double) -> Measure -> Measure
apply phi mu f = mu $ f . phi

Entonces ahora puede calcular fácilmente un lenguaje incrustado para las medidas. Aquí hay muchas cosas más que hacer, particularmente con respecto a espacios de muestra que no sean la línea real, dependencias entre variables aleatorias, condicionamiento, pero espero que entienda el punto.

En particular, el pushforward es funcionario:

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

También es una mónada (ejercicio - pista: esto se parece mucho a la mónada de continuación. ¿Qué es el return ? ¿Cuál es el análogo de call/cc ?).

Además, combinado con un marco de geometría diferencial, esto probablemente pueda convertirse en algo que calcule distribuciones posteriores Bayesianas automáticamente.

Al final del día, puedes escribir cosas como

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

para calcular la media de cos (X + Y) donde X tiene pdf gauss e Y ha sido muestreado por un método de MC cuyos resultados están en data .

Answer 10

Las distribuciones de probabilidad forman una mónada; ver, por ejemplo, el trabajo de Claire Jones y también el de LICS 1989, pero las ideas se remontan a un documento de Giry (DOI 10.1007 / BFb0092872) de 1982 y a una nota de 1962 de Lawvere que no puedo rastrear ( http: // permalink. gmane.org/gmane.science.mathematics.categories/6541 ).

Pero no veo el comonad: no hay forma de obtener una "a" de un "(a-> Double) -> Double". Tal vez si lo haces polimórfico - (a-> r) -> r para todo r? (Esa es la mónada de continuación).