Dev Faq
10
votes

Générer des nombres aléatoires avec une distribution probabiliste

OK, alors voici mon problème. Nous envisageons d'acheter un ensemble de données d'une entreprise pour augmenter notre ensemble de données existant. Aux fins de cette question, disons que ce jeu de données se classe avec un numéro organique (ce qui signifie que le numéro attribué à un endroit n'a aucune incidence sur le numéro attribué à un autre). La plage technique est de 0 à l'infini, mais à partir d'ensembles d'échantillons que j'ai vus, il est de 0 à 70. Sur la base de l'échantillon, il n'est certainement pas une distribution uniforme (sur 10 000, il y a peut-être 5 places avec un score de plus de 40 ans, 50 avec un score de plus de 10 et 1000 avec un score supérieur à 1). Avant de décider d'acheter cet ensemble, nous aimerions le simuler afin que nous puissions voir à quel point cela pourrait être utile.

Donc, pour le simuler, j'ai réfléchi à générer un nombre aléatoire pour chaque lieu (environ 150 000 nombres aléatoires). Mais, je veux aussi garder à l'esprit des données et garder la répartition relativement la même (ou au moins assez proche). Je saisis mon cerveau toute la journée en essayant de penser à un moyen de le faire et de venir vider.

Une pensée que j'avais devait carré le nombre aléatoire (entre 0 et sqrt (70)). Mais cela favoriserait à la fois moins de 1 et plus grand nombre.

Je pense que la répartition réelle devrait être hyperbolique dans le premier quadrant ... Je viens simplement sur la manière de transformer une distribution linéaire, même la distribution de nombres aléatoires en une distribution hyperbolique (si l'hyperbolie est même ce que je veux En premier lieu).

Toutes les pensées?

Donc, pour résumer, voici la distribution que je voudrais (environ):

  • 40 - 70: 0,02% - 0,05%
  • 10 - 40: 0,5% - 1%
  • 1 - 10: 10% - 20%
  • 0 - 1: reste (78,95% - 89,48%)
3 commentaires

J'ai trouvé ce glossaire de statistiques [ STATS.GLA.AC.UK/ étapes / glossaire / ... ]. Cela pourrait aider.

Je ne comprends pas vraiment. Avez-vous des chiffres de points flottants de 10k entre 0 et 70 que vous souhaitez distribuer sur un ensemble de 150K?

@Jonas elfström: Eh bien, l'inverse. Je veux générer des numéros de point variant aléatoires 150K avec la distribution spécifiée ...

4 Réponses :

8
votes

Il y a des années écrites pour PHP4, choisissez simplement votre distribution: 

<?php

define( 'RandomGaussian',           'gaussian' ) ;          //  gaussianWeightedRandom()
define( 'RandomBell',               'bell' ) ;              //  bellWeightedRandom()
define( 'RandomGaussianRising',     'gaussianRising' ) ;    //  gaussianWeightedRisingRandom()
define( 'RandomGaussianFalling',    'gaussianFalling' ) ;   //  gaussianWeightedFallingRandom()
define( 'RandomGamma',              'gamma' ) ;             //  gammaWeightedRandom()
define( 'RandomGammaQaD',           'gammaQaD' ) ;          //  QaDgammaWeightedRandom()
define( 'RandomLogarithmic10',      'log10' ) ;             //  logarithmic10WeightedRandom()
define( 'RandomLogarithmic',        'log' ) ;               //  logarithmicWeightedRandom()
define( 'RandomPoisson',            'poisson' ) ;           //  poissonWeightedRandom()
define( 'RandomDome',               'dome' ) ;              //  domeWeightedRandom()
define( 'RandomSaw',                'saw' ) ;               //  sawWeightedRandom()
define( 'RandomPyramid',            'pyramid' ) ;           //  pyramidWeightedRandom()
define( 'RandomLinear',             'linear' ) ;            //  linearWeightedRandom()
define( 'RandomUnweighted',         'non' ) ;               //  nonWeightedRandom()



function mkseed()
{
    srand(hexdec(substr(md5(microtime()), -8)) & 0x7fffffff) ;
}   //  function mkseed()




/*
function factorial($in) {
    if ($in == 1) {
        return $in ;
    }
    return ($in * factorial($in - 1.0)) ;
}   //  function factorial()


function factorial($in) {
    $out = 1 ;
    for ($i = 2; $i <= $in; $i++) {
        $out *= $i ;
    }

    return $out ;
}   //  function factorial()
*/




function random_0_1()
{
    //  returns random number using mt_rand() with a flat distribution from 0 to 1 inclusive
    //
    return (float) mt_rand() / (float) mt_getrandmax() ;
}   //  random_0_1()


function random_PN()
{
    //  returns random number using mt_rand() with a flat distribution from -1 to 1 inclusive
    //
    return (2.0 * random_0_1()) - 1.0 ;
}   //  function random_PN()




function gauss()
{
    static $useExists = false ;
    static $useValue ;

    if ($useExists) {
        //  Use value from a previous call to this function
        //
        $useExists = false ;
        return $useValue ;
    } else {
        //  Polar form of the Box-Muller transformation
        //
        $w = 2.0 ;
        while (($w >= 1.0) || ($w == 0.0)) {
            $x = random_PN() ;
            $y = random_PN() ;
            $w = ($x * $x) + ($y * $y) ;
        }
        $w = sqrt((-2.0 * log($w)) / $w) ;

        //  Set value for next call to this function
        //
        $useValue = $y * $w ;
        $useExists = true ;

        return $x * $w ;
    }
}   //  function gauss()


function gauss_ms( $mean,
                   $stddev )
{
    //  Adjust our gaussian random to fit the mean and standard deviation
    //  The division by 4 is an arbitrary value to help fit the distribution
    //      within our required range, and gives a best fit for $stddev = 1.0
    //
    return gauss() * ($stddev/4) + $mean;
}   //  function gauss_ms()


function gaussianWeightedRandom( $LowValue,
                                 $maxRand,
                                 $mean=0.0,
                                 $stddev=2.0 )
{
    //  Adjust a gaussian random value to fit within our specified range
    //      by 'trimming' the extreme values as the distribution curve
    //      approaches +/- infinity
    $rand_val = $LowValue + $maxRand ;
    while (($rand_val < $LowValue) || ($rand_val >= ($LowValue + $maxRand))) {
        $rand_val = floor(gauss_ms($mean,$stddev) * $maxRand) + $LowValue ;
        $rand_val = ($rand_val + $maxRand) / 2 ;
    }

    return $rand_val ;
}   //  function gaussianWeightedRandom()


function bellWeightedRandom( $LowValue,
                             $maxRand )
{
    return gaussianWeightedRandom( $LowValue, $maxRand, 0.0, 1.0 ) ;
}   //  function bellWeightedRandom()


function gaussianWeightedRisingRandom( $LowValue,
                                       $maxRand )
{
    //  Adjust a gaussian random value to fit within our specified range
    //      by 'trimming' the extreme values as the distribution curve
    //      approaches +/- infinity
    //  The division by 4 is an arbitrary value to help fit the distribution
    //      within our required range
    $rand_val = $LowValue + $maxRand ;
    while (($rand_val < $LowValue) || ($rand_val >= ($LowValue + $maxRand))) {
        $rand_val = $maxRand - round((abs(gauss()) / 4) * $maxRand) + $LowValue ;
    }

    return $rand_val ;
}   //  function gaussianWeightedRisingRandom()


function gaussianWeightedFallingRandom( $LowValue,
                                        $maxRand )
{
    //  Adjust a gaussian random value to fit within our specified range
    //      by 'trimming' the extreme values as the distribution curve
    //      approaches +/- infinity
    //  The division by 4 is an arbitrary value to help fit the distribution
    //      within our required range
    $rand_val = $LowValue + $maxRand ;
    while (($rand_val < $LowValue) || ($rand_val >= ($LowValue + $maxRand))) {
        $rand_val = floor((abs(gauss()) / 4) * $maxRand) + $LowValue ;
    }

    return $rand_val ;
}   //  function gaussianWeightedFallingRandom()


function logarithmic($mean=1.0, $lambda=5.0)
{
    return ($mean * -log(random_0_1())) / $lambda ;
}   //  function logarithmic()


function logarithmicWeightedRandom( $LowValue,
                                    $maxRand )
{
    do {
        $rand_val = logarithmic() ;
    } while ($rand_val > 1) ;

    return floor($rand_val * $maxRand) + $LowValue ;
}   //  function logarithmicWeightedRandom()


function logarithmic10( $lambda=0.5 )
{
    return abs(-log10(random_0_1()) / $lambda) ;
}   //  function logarithmic10()


function logarithmic10WeightedRandom( $LowValue,
                                      $maxRand )
{
    do {
        $rand_val = logarithmic10() ;
    } while ($rand_val > 1) ;

    return floor($rand_val * $maxRand) + $LowValue ;
}   //  function logarithmic10WeightedRandom()


function gamma( $lambda=3.0 )
{
    $wLambda = $lambda + 1.0 ;
    if ($lambda <= 8.0) {
        //  Use direct method, adding waiting times
        $x = 1.0 ;
        for ($j = 1; $j <= $wLambda; $j++) {
            $x *= random_0_1() ;
        }
        $x = -log($x) ;
    } else {
        //  Use rejection method
        do {
            do {
                //  Generate the tangent of a random angle, the equivalent of
                //      $y = tan(pi * random_0_1())
                do {
                    $v1 = random_0_1() ;
                    $v2 = random_PN() ;
                } while (($v1 * $v1 + $v2 * $v2) > 1.0) ;
                $y = $v2 / $v1 ;
                $s = sqrt(2.0 * $lambda + 1.0) ;
                $x = $s * $y + $lambda ;
            //  Reject in the region of zero probability
            } while ($x <= 0.0) ;
            //  Ratio of probability function to comparison function
            $e = (1.0 + $y * $y) * exp($lambda * log($x / $lambda) - $s * $y) ;
        //  Reject on the basis of a second uniform deviate
        } while (random_0_1() > $e) ;
    }

    return $x ;
}   //  function gamma()


function gammaWeightedRandom( $LowValue,
                              $maxRand )
{
    do {
        $rand_val = gamma() / 12 ;
    } while ($rand_val > 1) ;

    return floor($rand_val * $maxRand) + $LowValue ;
}   //  function gammaWeightedRandom()


function QaDgammaWeightedRandom( $LowValue,
                                 $maxRand )
{
    return round((asin(random_0_1()) + (asin(random_0_1()))) * $maxRand / pi()) + $LowValue ;
}   //  function QaDgammaWeightedRandom()


function gammaln($in)
{
    $tmp = $in + 4.5 ;
    $tmp -= ($in - 0.5) * log($tmp) ;

    $ser = 1.000000000190015
            + (76.18009172947146 / $in)
            - (86.50532032941677 / ($in + 1.0))
            + (24.01409824083091 / ($in + 2.0))
            - (1.231739572450155 / ($in + 3.0))
            + (0.1208650973866179e-2 / ($in + 4.0))
            - (0.5395239384953e-5 / ($in + 5.0)) ;

    return (log(2.5066282746310005 * $ser) - $tmp) ;
}   //  function gammaln()


function poisson( $lambda=1.0 )
{
    static $oldLambda ;
    static $g, $sq, $alxm ;

    if ($lambda <= 12.0) {
        //  Use direct method
        if ($lambda <> $oldLambda) {
            $oldLambda = $lambda ;
            $g = exp(-$lambda) ;
        }
        $x = -1 ;
        $t = 1.0 ;
        do {
            ++$x ;
            $t *= random_0_1() ;
        } while ($t > $g) ;
    } else {
        //  Use rejection method
        if ($lambda <> $oldLambda) {
            $oldLambda = $lambda ;
            $sq = sqrt(2.0 * $lambda) ;
            $alxm = log($lambda) ;
            $g = $lambda * $alxm - gammaln($lambda + 1.0) ;
        }
        do {
            do {
                //  $y is a deviate from a Lorentzian comparison function
                $y = tan(pi() * random_0_1()) ;
                $x = $sq * $y + $lambda ;
            //  Reject if close to zero probability
            } while ($x < 0.0) ;
            $x = floor($x) ;
            //  Ratio of the desired distribution to the comparison function
            //  We accept or reject by comparing it to another uniform deviate
            //  The factor 0.9 is used so that $t never exceeds 1
            $t = 0.9 * (1.0 + $y * $y) * exp($x * $alxm - gammaln($x + 1.0) - $g) ;
        } while (random_0_1() > $t) ;
    }

    return $x ;
}   //  function poisson()


function poissonWeightedRandom( $LowValue,
                                $maxRand )
{
    do {
        $rand_val = poisson() / $maxRand ;
    } while ($rand_val > 1) ;

    return floor($x * $maxRand) + $LowValue ;
}   //  function poissonWeightedRandom()


function binomial( $lambda=6.0 )
{
}


function domeWeightedRandom( $LowValue,
                             $maxRand )
{
    return floor(sin(random_0_1() * (pi() / 2)) * $maxRand) + $LowValue ;
}   //  function bellWeightedRandom()


function sawWeightedRandom( $LowValue,
                            $maxRand )
{
    return floor((atan(random_0_1()) + atan(random_0_1())) * $maxRand / (pi()/2)) + $LowValue ;
}   //  function sawWeightedRandom()


function pyramidWeightedRandom( $LowValue,
                               $maxRand )
{
    return floor((random_0_1() + random_0_1()) / 2 * $maxRand) + $LowValue ;
}   //  function pyramidWeightedRandom()


function linearWeightedRandom( $LowValue,
                               $maxRand )
{
    return floor(random_0_1() * ($maxRand)) + $LowValue ;
}   //  function linearWeightedRandom()


function nonWeightedRandom( $LowValue,
                            $maxRand )
{
    return rand($LowValue,$maxRand+$LowValue-1) ;
}   //  function nonWeightedRandom()




function weightedRandom( $Method,
                         $LowValue,
                         $maxRand )
{
    switch($Method) {
        case RandomGaussian         :
            $rVal = gaussianWeightedRandom( $LowValue, $maxRand ) ;
            break ;
        case RandomBell             :
            $rVal = bellWeightedRandom( $LowValue, $maxRand ) ;
            break ;
        case RandomGaussianRising   :
            $rVal = gaussianWeightedRisingRandom( $LowValue, $maxRand ) ;
            break ;
        case RandomGaussianFalling  :
            $rVal = gaussianWeightedFallingRandom( $LowValue, $maxRand ) ;
            break ;
        case RandomGamma            :
            $rVal = gammaWeightedRandom( $LowValue, $maxRand ) ;
            break ;
        case RandomGammaQaD         :
            $rVal = QaDgammaWeightedRandom( $LowValue, $maxRand ) ;
            break ;
        case RandomLogarithmic10    :
            $rVal = logarithmic10WeightedRandom( $LowValue, $maxRand ) ;
            break ;
        case RandomLogarithmic      :
            $rVal = logarithmicWeightedRandom( $LowValue, $maxRand ) ;
            break ;
        case RandomPoisson          :
            $rVal = poissonWeightedRandom( $LowValue, $maxRand ) ;
            break ;
        case RandomDome             :
            $rVal = domeWeightedRandom( $LowValue, $maxRand ) ;
            break ;
        case RandomSaw              :
            $rVal = sawWeightedRandom( $LowValue, $maxRand ) ;
            break ;
        case RandomPyramid          :
            $rVal = pyramidWeightedRandom( $LowValue, $maxRand ) ;
            break ;
        case RandomLinear           :
            $rVal = linearWeightedRandom( $LowValue, $maxRand ) ;
            break ;
        default                     :
            $rVal = nonWeightedRandom( $LowValue, $maxRand ) ;
            break ;
    }

    return $rVal;
}

?>

6 commentaires

Merci pour le code. Cependant, j'ai essayé de regarder toutes les méthodes que vous avez fournies et que je n'ai pas vu cela semblait correspondre à mon modèle. Les statistiques n'ont jamais été mon point fort. Si vous pouviez signaler un modèle que vous pensez peut-être, je serais toutes les oreilles ... merci ...

Une option serait d'essayer de générer une série de valeurs et de les traquer sur un graphique à l'aide de chacune des différentes distributions prédéfinies pour voir à quoi ressemblent les courbes. Wikipedia dispose également de nombreuses entrées sur plusieurs de ces distributions ..... Cependant, pour ce que vous décrivez (si je l'ai interprété correctement), essayez GaussianPetherisingRandom si vous souhaitez plus de valeurs supérieures de la gamme ou de la chute de poids gaussiian si vous souhaitez plus de valeurs de plage plus basse. .. bien que Poisson est souvent une méthode utile pour de nombreuses situations du monde réel

Ok, j'ai essayé chacun. Le GaussianEdPoldedFallingRomingRandom est le plus proche, mais il ne tombe toujours pas assez vite (200 au lieu de 5 plus de 40, 5700 au lieu de 50 plus de 10 et 9500 au lieu de 1000 sur 1. J'ai essayé CSCH et il semble beaucoup plus proche ( Comme il correspond aux hautes plages), mais tombe trop vite au milieu. Pensées?

Dans ce cas, cherchez-vous à utiliser GaussianEdRandom ($ lowvalue, $ maxrand, $ moyenne, $ STDDEV), mais définissant vos propres valeurs pour une écart moyen et standard, ou modifiez l'appel à Gauss () dans GaussianEdPartedFallingRandom () pour appeler Gauss_Ms ($ moyenne , $ STDDEV) avec vos propres valeurs pour une écart-type moyen et type. Cela peut prendre des expériences ... mais regardez la page Wikipedia pour voir comment les modifications apportées à ces paramètres affectent la forme de la courbe

@MarkBaker excellente ressource !! Je vois que la poste est déjà assez ancienne, mais puis-je demander si une fonction de distribution discrète basée sur les données empiriques pourrait ne pas être meilleure qu'une fonction théorique ici?

@ hek2mgl - Si je l'écrivais maintenant, je le ferais très différemment

10
votes

Regardez les distributions utilisées dans l'analyse de fiabilité - ils ont tendance à avoir ces longues queues. Une possibilité relativement simplement simplement est la distribution Weibull avec p (x> x) = exp [- (x / b) ^ a].

ajustement de vos valeurs comme p (x> 1) = 0,1 et p (x> 10) = 0,005, je reçois un = 0,36 et b = 0,1. Cela impliquerait que p (x> 40) * 10000 = 1,6, qui est un peu trop bas, mais p (x> 70) * 10000 = 0,2 qui est raisonnable.

Modifier Oh, et pour générer une variable aléatoire distribuée de Weibull à partir d'une valeur uniforme (0,1) u, ​​calculez simplement B * [- journal (1-u)] ^ (1 / a). Ceci est la fonction inverse de 1-p (x> x) au cas où je suis mal calculé quelque chose.

2 commentaires

WOW, qui semble presque identique au résultat que je suis après (4> 40, 60> 10, 1030> 1). Excellent! Merci!

Vous pouvez mettre à jour votre réponse avec le code PHP de travail?

0
votes

Cette façon naïve de le faire va probablement incliner la distribution d'une certaine manière que je ne peux pas voir en ce moment. L'idée est simplement de parcourir votre premier jeu de données, trié et comme paires. Puis randomise 15 nouveaux chiffres entre chaque paire pour obtenir le nouveau tableau.

EXEMPLE RUBY, car je ne parle pas beaucoup de PHP. Espérons qu'une telle idée simple devrait être facile à traduire en php. xxx

0 commentaires
3
votes

Le moyen le plus facile (mais pas très efficace) de générer des nombres aléatoires qui suivent une distribution donnée est une technique appelée VON NEUMANN Rejection .

La simple exploration de la technique est celle-ci. Créez une case qui renferme complètement votre distribution. (Appelons votre distribution f ) puis choisissez un point aléatoire (x, y) dans la case. Si y , puis utilisez x comme numéro aléatoire. Si y> f (x) , supprimez les deux x et y et choisissez un autre point. Continuez jusqu'à ce que vous ayez une quantité suffisante de valeurs à utiliser. Les valeurs de x que vous ne rejetez pas seront distribuées selon f .

1 commentaires

Sauf si je ne me trompe pas, n'est-ce pas que d'obtenir des points aléatoires sous une courbe définie par f (x) ? Considérant que ma courbe a l'air hyperbolique, la plus grande densité de points serait autour de l'origine, de sorte que les chiffres générés ne seraient pas biaisés vers le milieu de la boîte bornée créée entre l'origine et le sommet (et ne favorisent donc pas les nombres inférieurs. J'en ai besoin de)?

Articles qui pourrait vous intéresser :

Pourquoi est-ce que j'obtiens une erreur laravel lors de la création d'un nouveau projet?
Select2 ne s'affiche pas correctement dans le mode bootstrap
Laravel Echo n'écoute pas
La validation de la demande de formulaire ne fonctionne pas. Il charge juste la page d'accueil dans Postman