¿Cómo genero puntos que coinciden con un histograma?

StackOverflow https://stackoverflow.com/questions/423006

Pregunta

Estoy trabajando en un sistema de simulación. Pronto tendré datos experimentales (histogramas) para la distribución de valores en el mundo real para varias entradas de simulación.

Cuando se ejecuta la simulación, me gustaría poder producir valores aleatorios que coincidan con la distribución medida. Preferiría hacer esto sin almacenar los histogramas originales. ¿Cuáles son algunas buenas maneras de

  1. ¿Asignar un histograma a un conjunto de parámetros que representan la distribución?
  2. ¿Generando valores basados ??en esos parámetros en tiempo de ejecución?

EDITAR: Los datos de entrada son duraciones de eventos para varios tipos diferentes de eventos. Espero que diferentes tipos tengan diferentes funciones de distribución.

¿Fue útil?

Solución

Al menos dos opciones:

  1. Integre el histograma e inviértalo numéricamente.
  2. Rechazo

Integración numérica

De Computación en física moderna por William R. Gibbs:

  

Uno siempre puede integrar numéricamente [la función] e invertir el [ cdf ]   pero esto a menudo no es muy satisfactorio, especialmente si el pdf está cambiando   rápidamente.

Literalmente construye una tabla que traduce el rango [0-1) a los rangos apropiados en la distribución de destino. Luego, lance su PRNG habitual (de alta calidad) y traduzca con la tabla. Es engorroso, pero claro, realizable y completamente general.

Rechazo:

Normalice el histograma de destino, luego

  1. Lance los dados para elegir una posición ( x ) a lo largo del rango al azar.
  2. Lance nuevamente y seleccione este punto si el nuevo número aleatorio es menor que el histograma normalizado en esta bandeja. De lo contrario goto (1).

De nuevo, de mente simple pero clara y trabajadora. Puede ser lento para la distribución con mucha probabilidad muy baja (picos con colas largas).

Con estos dos métodos, puede aproximar los datos con ajustes polinomiales por partes o splines para generar una curva suave si no se desea un histograma de función escalonada --- pero deje eso para más adelante Puede ser una optimización prematura.

Pueden existir mejores métodos para casos especiales.

Todo esto es bastante estándar y debería aparecer en cualquier libro de texto de Análisis Numérico si se necesitan más detalles.

Otros consejos

Más información sobre el problema sería útil. Por ejemplo, ¿en qué tipo de valores se encuentran los histogramas? ¿Son categóricos (por ejemplo, colores, letras) o continuos (por ejemplo, alturas, tiempo)?

Si los histogramas están sobre datos categóricos, creo que puede ser difícil parametrizar las distribuciones a menos que haya muchas correlaciones entre categorías.

Si los histogramas están sobre datos continuos, podría intentar ajustar la distribución utilizando mezclas de gaussianos. Es decir, intente ajustar el histograma utilizando $ \ sum_ {i = 1} ^ n w_i N (m_i, v_i) $ donde m_i y v_i son la media y la varianza. Luego, cuando desee generar datos, primero muestrea una i de 1..n con probabilidad proporcional a los pesos w_i y luego muestrea una x ~ n (m_i, v_i) como lo haría con cualquier gaussiano.

De cualquier manera, es posible que desee leer más sobre modelos de mezcla .

Entonces, parece que lo que quiero para generar una distribución de probabilidad dada es una Función de Quantile , que es el inverso de la función de distribución acumulada , como dice @dmckee.

La pregunta es: ¿Cuál es la mejor manera de generar y almacenar una función de cuantiles que describa un histograma continuo dado? Tengo la sensación de que la respuesta dependerá en gran medida de la forma de la entrada: si sigue algún tipo de patrón, debería haber simplificaciones en el caso más general. Voy a actualizar aquí a medida que voy.

Editar:

Esta semana tuve una conversación que me recordó este problema. Si renuncio a describir el histograma como una ecuación y solo almaceno la tabla, ¿puedo hacer selecciones en O (1) tiempo? Resulta que puede, sin ninguna pérdida de precisión, a costa del tiempo de construcción de O (N lgN).

Crea una matriz de N elementos. Una selección aleatoria uniforme en la matriz encontrará un elemento con probabilidad 1 / N. Para cada elemento, almacene la fracción de aciertos para los que este elemento debería seleccionarse realmente, y el índice de otro elemento que se seleccionará si este no lo está.

Muestreo aleatorio ponderado, implementación de C: 

//data structure
typedef struct wrs_data {
  double share; 
  int pair;
  int idx;
} wrs_t;


//sort helper
int wrs_sharecmp(const void* a, const void* b) {
  double delta = ((wrs_t*)a)->share - ((wrs_t*)b)->share;
  return (delta<0) ? -1 : (delta>0);
}


//Initialize the data structure
wrs_t* wrs_create(int* weights, size_t N) {
  wrs_t* data = malloc(sizeof(wrs_t));
  double sum = 0;
  int i;
  for (i=0;i<N;i++) { sum+=weights[i]; }
  for (i=0;i<N;i++) {
    //what percent of the ideal distribution is in this bucket?
    data[i].share = weights[i]/(sum/N); 
    data[i].pair = N;
    data[i].idx = i;
  }
  //sort ascending by size
  qsort(data,N, sizeof(wrs_t),wrs_sharecmp);

  int j=N-1; //the biggest bucket
  for (i=0;i<j;i++) {
    int check = i;
    double excess = 1.0 - data[check].share;
    while (excess>0 && i<j) {
      //If this bucket has less samples than a flat distribution,
      //it will be hit more frequently than it should be.  
      //So send excess hits to a bucket which has too many samples.
      data[check].pair=j; 
      // Account for the fact that the paired bucket will be hit more often,
      data[j].share -= excess;  
      excess = 1.0 - data[j].share;
      // If paired bucket now has excess hits, send to new largest bucket at j-1
      if (excess >= 0) { check=j--;} 
    }
  }
  return data;
}


int wrs_pick(wrs_t* collection, size_t N)
//O(1) weighted random sampling (after preparing the collection).
//Randomly select a bucket, and a percentage.
//If the percentage is greater than that bucket's share of hits, 
// use it's paired bucket.
{
  int idx = rand_in_range(0,N);
  double pct = rand_percent();
  if (pct > collection[idx].share) { idx = collection[idx].pair; }
  return collection[idx].idx;
}

Edición 2:   Después de una pequeña investigación, encontré que incluso es posible hacer la construcción en tiempo O (N). Con un seguimiento cuidadoso, no es necesario ordenar la matriz para encontrar los contenedores grandes y pequeños. Implementación actualizada aquí

Si necesita obtener una gran cantidad de muestras con una distribución ponderada de puntos discretos, mire una respuesta a una pregunta similar .

Sin embargo, si necesita aproximarse a alguna función aleatoria continua utilizando un histograma, entonces su mejor apuesta es probablemente la respuesta de integración numérica de dmckee. Alternativamente, puede usar el alias, y almacenar el punto a la izquierda, y elegir un número uniforme entre los dos puntos.

Para elegir un histograma (original o reducido), método de alias de Walker Es rápido y simple.

No afiliado a StackOverflow
scroll top