Comment générer des points correspondant à un histogramme?

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

Question

Je travaille sur un système de simulation. J'aurai bientôt des données expérimentales (histogrammes) pour la distribution des valeurs dans le monde réel pour plusieurs entrées de simulation.

Lors de l'exécution de la simulation, j'aimerais pouvoir générer des valeurs aléatoires correspondant à la distribution mesurée. Je préférerais le faire sans stocker les histogrammes d'origine. Quels sont les bons moyens de

  1. Mapper un histogramme sur un ensemble de paramètres représentant la distribution?
  2. Générer des valeurs basées sur ces paramètres lors de l'exécution?

EDIT: Les données en entrée sont des durées d’événement pour plusieurs types d’événements. Je pense que différents types auront différentes fonctions de distribution.

Était-ce utile?

La solution

Au moins deux options:

  1. Intégrez l'histogramme et inversez numériquement.
  2. Rejet

Intégration numérique

Extrait du calcul en physique moderne de William R. Gibbs:

  

On peut toujours intégrer numériquement [la fonction] et inverser le [ cdf ]   mais cela n’est souvent pas très satisfaisant, surtout si le pdf est en train de changer   rapidement.

Vous construisez littéralement une table qui traduit la plage [0-1) en plages appropriées dans la distribution cible. Ensuite, jetez votre PRNG habituel (haute qualité) et traduisez-le avec la table. C'est lourd, mais il est clair, pratique et tout à fait général.

Rejet:

Normalisez l'histogramme cible, puis

  1. Lancez les dés pour choisir une position ( x ) le long de la plage de manière aléatoire.
  2. Lancez à nouveau et sélectionnez ce point si le nouveau nombre aléatoire est inférieur à l'histogramme normalisé de cette case. Sinon, allez à (1).

Encore une fois, simple d'esprit mais clair et efficace. Il peut être lent pour une distribution avec une probabilité très faible (pics avec de longues queues).

Avec ces deux méthodes, vous pouvez approximer les données avec des ajustements polynomiaux par morceaux ou des splines pour générer une courbe lisse si un histogramme de fonction d'étape n'est pas souhaité --- mais laissez-le pour plus tard ce peut être une optimisation prématurée.

De meilleures méthodes peuvent exister pour des cas particuliers.

Tout cela est assez standard et devrait apparaître dans n'importe quel manuel d'analyse numérique si j'ai besoin de plus de détails.

Autres conseils

Plus d'informations sur le problème seraient utiles. Par exemple, sur quelle sorte de valeurs les histogrammes sont-ils traités? Sont-ils catégoriques (couleurs, lettres, par exemple) ou continus (hauteur, temps, par exemple)?

Si les histogrammes sont trop de données catégorielles, je pense qu'il peut être difficile de paramétrer les distributions sauf s'il existe de nombreuses corrélations entre les catégories.

Si les histogrammes contiennent des données continues, essayez d’ajuster la distribution à l’aide de mélanges de gaussiennes. En d’autres termes, essayez d’ajuster l’histogramme en utilisant un $ \ sum_ {i = 1} ^ n w_i N (m_i, v_i) $ où m_i et v_i sont la moyenne et la variance. Ensuite, lorsque vous souhaitez générer des données, vous devez d'abord échantillonner un i de 1..n avec une probabilité proportionnelle aux poids w_i, puis échantillonner un x ~ n (m_i, v_i) comme vous le feriez avec n'importe quel gaussien.

Dans tous les cas, vous voudrez peut-être en savoir plus sur les modèles de mélange .

Il semble donc que ce que je veux pour générer une distribution de probabilité donnée soit une Fonction Quantile , qui est l'inverse de la fonction de distribution cumulative , comme l'indique @dmckee.

La question est la suivante: quel est le meilleur moyen de générer et de stocker une fonction quantile décrivant un histogramme continu donné? J'ai le sentiment que la réponse dépendra beaucoup de la forme de l'entrée - si elle suit n'importe quel modèle, il devrait y avoir des simplifications par rapport au cas le plus général. Je mettrai à jour ici que je vais.

Modifier:

J'ai eu une conversation cette semaine qui m'a rappelé ce problème. Si je ne décris pas l'histogramme comme une équation et que je ne stocke que le tableau, puis-je effectuer des sélections dans le temps O (1)? Il s'avère que vous pouvez, sans perte de précision, au prix de O (N lgN) le temps de construction.

Créez un tableau de N éléments. Une sélection aléatoire uniforme dans le tableau trouvera un élément avec une probabilité de 1 / N. Pour chaque élément, stockez la fraction de correspondances pour laquelle cet élément doit réellement être sélectionné, ainsi que l'index d'un autre élément qui sera sélectionné si celui-ci ne l'est pas.

Échantillonnage aléatoire pondéré, mise en œuvre du code 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;
}

Modifier 2:   Après quelques recherches, j'ai trouvé qu'il était même possible de faire la construction en temps O (N). Avec un suivi minutieux, vous n'avez pas besoin de trier le tableau pour trouver les grands et les petits bacs. mise en œuvre mise à jour ici

Si vous devez extraire un grand nombre d’échantillons avec une distribution pondérée de points discrets, consultez une réponse à une question similaire .

Toutefois, si vous devez approximer une fonction aléatoire continue à l'aide d'un histogramme, votre meilleure option est probablement la réponse d'intégration numérique de dmckee. Vous pouvez également utiliser le crénelage, stocker le point à gauche et choisir un nombre uniforme entre les deux points.

Pour choisir parmi un histogramme (original ou réduit), La méthode d'alias de Walker est rapide et simple.

Pour une distribution normale, les éléments suivants peuvent aider:

http://en.wikipedia.org/wiki/Normal_distribution#Generating_values_for_normal_for_normal_rand_page

Non affilié à StackOverflow
scroll top