Ik moet willekeurige getallen genereren volgens de normale verdeling binnen het interval $ (a, b) $. (Ik werk in R.)
Ik weet dat de functie rnorm(n,mean,sd)
willekeurige getallen genereert volgens de normale verdeling, maar hoe stel je de intervallimieten daarbinnen in? Zijn daar bepaalde R-functies voor beschikbaar?
Reacties
Antwoord
Het klinkt alsof u wilt simuleren vanuit een afgekapte distributie , en in uw specifieke voorbeeld , een afgekapt normaal .
Er zijn verschillende methoden om dit te doen, sommige eenvoudig, sommige relatief efficiënt.
Ik zal enkele benaderingen aan uw normale voorbeeld illustreren.
-
Hier is een zeer eenvoudige methode om een voor een te genereren (in een soort pseudocode ):
$ \ tt {repeat} $ genereer $ x_i $ van N (mean, sd) $ \ tt {tot} $ lager $ \ leq x_i \ leq $ hoger
Als de meeste distributie binnen de grenzen valt, is dit redelijk, maar het kan behoorlijk traag worden als je genereert bijna altijd buiten de limieten.
In R zou je de een-voor-een-lus kunnen vermijden door het gebied binnen de grenzen te berekenen en voldoende waarden te genereren dat je er vrijwel zeker van kunt zijn dat je na het weggooien de waarden buiten de grenzen had je nog steeds zoveel waarden als nodig.
-
Je zou accept-reject kunnen gebruiken met een geschikte majorizing-functie over het interval (in sommige gevallen goed genoeg zijn). Als de limieten redelijk smal waren ten opzichte van de s.d. maar je zat “niet ver in de staart, een uniforme majorisering zou bijvoorbeeld goed werken met de normale.
-
Als je een redelijk efficiënte cdf en inverse cdf hebt (zoals
pnorm
enqnorm
voor de normale distributie in R) kunt u de inverse-cdf-methode gebruiken die wordt beschreven in de eerste alinea van het simulatiegedeelte van de Wikipedia-pagina op de afgekapte normaal . [In feite dit is hetzelfde als het nemen van een afgekapt uniform (afgekapt op de vereiste kwantielen, wat eigenlijk helemaal geen afwijzingen vereist, aangezien dat gewoon een ander uniform is) en daarop de inverse normale cdf toepassen. Merk op dat dit kan mislukken als je “ver in de staart zit] -
Er zijn andere benaderingen; dezelfde Wikipedia-pagina vermeldt de aanpassing van de ziggurat -methode, die zou moeten werken voor verschillende distributies.
De dezelfde Wikipedia-link vermeldt twee specifieke pakketten (beide op CRAN) met functies voor het genereren van afgekapte normalen:
Het
MSM
-pakket in R heeft een functie,rtnorm
, die de trunciatie berekent op basis van een afgekapte normaal. Hettruncnorm
-pakket in R heeft ook functies om te tekenen uit een afgekapte normaal.
Als je rondkijkt, wordt veel hiervan behandeld in antwoorden op andere vragen (maar niet precies duplicaten aangezien deze vraag algemener is dan alleen de afgekapte normaal) … zie aanvullende discussie in
a. Dit antwoord
b. Xi “an” s antwoord hier , dat een link heeft naar zijn arXiv-paper (samen met enkele andere waardevolle reacties).
Answer
De snelle aanpak is om de 68-95-99.7 regel .
In een normale verdeling valt 99,7% van de waarden binnen 3 standaarddeviaties van het gemiddelde. Dus als u uw gemiddelde instelt op het midden van uw gewenste minimumwaarde en maximumwaarde, en uw standaarddeviatie instelt op 1/3 van uw gemiddelde, krijgt u (meestal) waarden die binnen het gewenste interval vallen. Dan kun je de rest gewoon opschonen.
minVal <- 0 maxVal <- 100 mn <- (maxVal - minVal)/2 # Generate numbers (mostly) from min to max x <- rnorm(count, mean = mn, sd = mn/3) # Do something about the out-of-bounds generated values x <- pmax(minVal, x) x <- pmin(maxVal, x)
Ik heb onlangs hetzelfde probleem ondervonden toen ik probeerde te genereren willekeurige studentcijfers voor testgegevens. In de bovenstaande code “heb ik pmax
en pmin
gebruikt om waarden buiten het bereik te vervangen door de min of max in-bounds waarde.Dit werkt voor mijn doel, omdat ik vrij kleine hoeveelheden gegevens genereer, maar voor grotere hoeveelheden krijg je merkbare hobbels bij de min- en max-waarden. Dus afhankelijk van je doeleinden is het misschien beter om die waarden te negeren, ze te vervangen met NA
s, of “herrol” ze totdat ze “binnenkomen.
Reacties
- Waarom zou je dit doen? Het is zo eenvoudig om normale willekeurige getallen te genereren en de getallen die moeten worden afgekapt weg te laten, dat het niet ‘ nodig is om er ingewikkeld over te zijn, tenzij de gewenste afkapping bijna 100% van de oppervlakte bedraagt van de dichtheid.
- Misschien heb ik ‘ de oorspronkelijke vraag verkeerd geïnterpreteerd. Ik kwam deze vraag tegen toen ik probeerde uit te zoeken hoe ik een niet-direct-stats-gerelateerde programmeertaak in R kon bereiken, en ik ‘ heb nu pas opgemerkt dat deze pagina een statistiek-uitwisseling is , niet een programmeer-stackexchange. 🙂 In mijn geval wilde ik een bepaald aantal willekeurige gehele getallen genereren, met waarden variërend van 0 tot 100, en ik wilde dat de gegenereerde waarden over dat bereik op een mooie belcurve zouden vallen. Sinds ik dit heb geschreven, ‘ realiseerde ik me dat
sample(x=min:max, prob=dnorm(...))
misschien een gemakkelijkere manier is om dat te doen. - @Glen_b Aaron Wells noemt
sample(x=min:max, prob=dnorm(...))
wat een beetje korter lijkt dan je antwoord. - Merk echter op dat de
sample()
-truc alleen nuttig is als je ‘ probeert willekeurige gehele getallen te kiezen, of een andere set discrete, vooraf gedefinieerde waarden.
Antwoord
Geen van de antwoorden hier geeft een efficiënte methode voor het genereren van afgekapte normale variabelen die geen afwijzing van willekeurig grote aantal gegenereerde waarden. Als u waarden wilt genereren uit een afgekapte normale verdeling, met gespecificeerde onder- en bovengrenzen $ a < b $ , dan kan — zonder afwijzing — worden gedaan door uniforme kwantielen te genereren over het kwantielbereik dat is toegestaan door de afkapping, en inverse transformatie-steekproeven te gebruiken om overeenkomstige normale waarden te krijgen .
Laat $ \ Phi $ de CDF van de standaard normale distributie aangeven. We willen $ X_1, …, X_N $ genereren op basis van een afgekapte normale verdeling (met gemiddelde parameter $ \ mu $ en variantieparameter $ \ sigma ^ 2 $ ) $ ^ \ dagger $ met lagere en bovenste truncatiegrenzen $ a < b $ . Dit kan als volgt worden gedaan:
$$ X_i = \ mu + \ sigma \ cdot \ Phi ^ {- 1} (U_i) \ quad \ quad \ quad U_1, …, U_N \ sim \ text {IID U} \ Big [\ Phi \ Big (\ frac {a- \ mu} {\ sigma} \ Big), \ Phi \ Big (\ frac {b- \ mu} {\ sigma} \ Big) \ Big]. $$
Er is geen ingebouwde functie voor gegenereerde waarden uit de afgekapte distributie, maar het is triviaal om deze methode te programmeren met de gewone functies voor het genereren van willekeurige variabelen. Hier is een simpele R
functie rtruncnorm
die deze methode implementeert in een paar regels code.
rtruncnorm <- function(N, mean = 0, sd = 1, a = -Inf, b = Inf) { if (a > b) stop("Error: Truncation range is empty"); U <- runif(N, pnorm(a, mean, sd), pnorm(b, mean, sd)); qnorm(U, mean, sd); }
Dit is een gevectoriseerde functie die N
IID willekeurige variabelen zal genereren uit de afgekapte normale verdeling. Het zou gemakkelijk zijn om via dezelfde methode functies voor andere ingekorte distributies te programmeren. Het zou ook niet zo moeilijk zijn om bijbehorende dichtheids- en kwantielfuncties te programmeren voor de afgekapte distributie.
$ ^ \ dagger $ Merk op dat de afkapping verandert het gemiddelde en de variantie van de verdeling, dus $ \ mu $ en $ \ sigma ^ 2 $ zijn niet het gemiddelde en de variantie van de afgekapte verdeling.
Antwoord
Drie manieren hebben voor mij gewerkt:
-
door sample () te gebruiken met rnorm ():
sample(x=min:max, replace= TRUE, rnorm(n, mean))
-
met behulp van het msm-pakket en de rtnorm-functie:
rtnorm(n, mean, lower=min, upper=max)
-
de rnorm () gebruiken en de onder- en bovengrenzen specificeren, zoals Hugh hierboven heeft gepost:
sample <- rnorm(n, mean=mean); sample <- sample[x > min & x < max]
x <- rnorm(n, mean, sd); x <- x[x > lower.limit & x < upper.limit]