Random numbers
You are encouraged to solve this task according to the task description, using any language you may know.
- Task
Generate a collection filled with 1000 normally distributed random (or pseudo-random) numbers with a mean of 1.0 and a standard deviation of 0.5
Many libraries only generate uniformly distributed random numbers. If so, you may use one of these algorithms.
- Related task
Ada
<lang ada>with Ada.Numerics; use Ada.Numerics; with Ada.Numerics.Float_Random; use Ada.Numerics.Float_Random; with Ada.Numerics.Elementary_Functions; use Ada.Numerics.Elementary_Functions;
procedure Normal_Random is
function Normal_Distribution ( Seed : Generator; Mu : Float := 1.0; Sigma : Float := 0.5 ) return Float is begin return Mu + (Sigma * Sqrt (-2.0 * Log (Random (Seed), 10.0)) * Cos (2.0 * Pi * Random (Seed))); end Normal_Distribution; Seed : Generator; Distribution : array (1..1_000) of Float;
begin
Reset (Seed); for I in Distribution'Range loop Distribution (I) := Normal_Distribution (Seed); end loop;
end Normal_Random;</lang>
ALGOL 68
<lang algol68>PROC random normal = REAL: # normal distribution, centered on 0, std dev 1 # (
sqrt(-2*log(random)) * cos(2*pi*random)
);
test:(
[1000]REAL rands; FOR i TO UPB rands DO rands[i] := 1 + random normal/2 OD; INT limit=10; printf(($"("n(limit-1)(-d.6d",")-d.5d" ... )"$, rands[:limit]))
)</lang>
- Output:
( 0.693461, 0.948424, 0.482261, 1.045939, 0.890818, 1.467935, 0.604153, 0.804811, 0.690227, 0.83462 ... )
Arturo
<lang rebol>rnd: function []-> (random 0 10000)//10000
rands: map 1..1000 'x [
1 + (sqrt neg 2 * ln rnd) * (cos 2 * pi * rnd)
]
print rands</lang>
- Output:
0.6219537961087694 1.279396486161406 1.619019280815647 2.119538294228789 -0.1598383851981044 2.67797211803156 0.9304703037226587 1.629254364659528 -0.4171704717398712 0.9082342931486092 -0.5929704390625219 2.117000897984871 -0.1981633787460266 0.01132471973856153 2.102359263212924 0.2408823232884222 2.046195035792376 0.6374831627030295 0.000839808324124558 1.117061838266626 0.7413355299469649 0.4485598815755762 2.999434800016997 2.560580541932842 1.703197984879731 2.889159248353575 -1.800157205708138 1.756810020187321 0.7136708180852145 0.5929151678321705 0.332519993787973 2.660212054362758 0.5835660585480075 0.8527946892567934 1.640573993747053 0.09471843345263908 1.051402997891346 1.116149156137905 -0.7400139019343499 1.782572831979232 2.531779039786426 0.5240268064639871 0.07099232630526586 -0.854892656700071 1.54381929430469 -0.4416899008614745 0.4274356035015117 0.7350027625573482 2.153583935076981 1.461215281535983 -1.041723064151266 2.338060763553139 -0.1492967916030414 0.3799517724040202 0.4577924541353815 0.673317567666373 -2.27731583876462 1.28480355806061 -0.6925811023772748 -0.2642224122781984 0.6590513830891744 2.55537133425143 1.67933335469247 0.8659013395355968 -1.211026941441126 0.9524579534222226 -0.1931750631835656 -0.5119479869237693 0.1814749003063878 3.03139579963414 0.8041947090501279 -0.01865777860269424 -1.345117694187226 1.172047503199659 -0.005768598737427855 0.4430753828101858 -0.04896736045307204 0.4895169217506786 1.048229080758918 2.149038787105586 2.027034431271556 2.432282801435372 2.53693570310459 -0.5968142206719835 1.445180766446883 1.993843935917086 0.3246477684427571 0.2214120156342799 2.123978740745019 0.7660225685637785 -0.4679548030089944 1.892651338960027 -1.082672493660939 1.388650546130418 0.8586892232800707 0.1688680513907505 0.2031720169757011 1.0622153495636 2.237489272627846 0.489619531018731 2.067023562273573 3.232122397451699 -0.03837283616929521 -0.8782247037247171 2.155998134902104 0.5475989225547847 1.468927649628817 1.27719993843459 1.237163843541324 1.73794547561703 0.8804315237620751 0.02508158115251735 1.65857224365674 1.078989935949705 0.7164357441628163 1.71375577930128 0.6629629974248089 1.584782749537475 0.5168549101430115 2.435706381815515 2.458518058181382 0.4983726815666553 1.662610001403874 0.5100242881482009 0.424252537753539 2.874829969968073 2.739428722137513 1.568599755919953 1.17450447863347 0.2147990014354625 0.94960001217684 1.530842927153068 1.033241575011026 0.7136369721788152 -1.752151570685991 -0.1825011364213669 0.7632351752546898 -0.705217801387372 1.146249503264144 0.5024360662432521 -0.05514204232437669 2.751909795371251 1.404999848001164 0.682474910779916 2.333698852277054 -0.3385164481257155 0.7477589125734123 0.7499986390595292 1.357508529051189 0.5882691367699802 0.6512207531472127 0.7682407661401032 1.256155994922453 0.5223880933125757 4.397210245005706 0.1876049987895129 -0.4165321555980093 0.1177916886895712 1.418593319932079 0.6245130605510126 0.03422663901416945 -0.126508344317811 0.6975972877438945 1.322572198342049 1.386458671539238 2.18031821792089 0.718913947303726 -0.4766346702601931 3.536427477250057 1.346365414270038 1.846533142607893 1.569841840872839 1.802939152281028 0.8756670850485474 0.951865353598341 0.958791363972665 1.157743287178253 1.172276335245585 0.7118514688347454 0.04596669908479134 2.160216781880242 -0.4558470112795174 -0.1408698236013977 2.55863514213935 -0.9762323021069998 0.3023059074783826 0.1006358104014331 1.798240854639 1.796303359312389 0.2854804374023457 1.695153128919529 -0.7461161123527316 -0.5717889809615559 -1.22969188461553 -0.3124307656341108 -0.4572231441070826 2.196150179180091 0.9920452455196852 0.7298941327650366 -0.009138175032543883 1.196326721003037 3.029519138874344 2.053193497411489 0.3020035894567694 1.885130148538693 0.3571024527100635 -0.006202951752455332 2.422152014091419 1.953976188093006 -0.3035872780193656 1.038990844201456 1.15759942528698 0.8840631785390038 1.006992878003682 0.7844353853636307 0.8614644612780162 0.9999999999999998 1.37819707688161 1.451946740030631 -0.2701312811593382 2.470927310260342 0.188450176887656 1.552249295530555 -0.4011035144371413 0.4812612270534884 3.533385275058706 0.01878389160958915 2.03329739376812 1.497360243363507 0.01291962460839935 2.441040065770387 2.45750172056598 1.759298261864295 1.828459550040343 0.1185298953877549 0.7656071379341856 1.303441760675321 0.3049733552767356 0.3402213592635629 2.653472185103094 0.5024574596798626 1.455823625970243 0.2952856940037123 1.716812084833785 1.01761544215552 -0.06454420761699553 -0.1136659519769139 0.7720535027210765 0.2077001368790136 2.906882616021471 0.5692803618983897 1.402678021067997 1.281703298994574 1.312617838849445 2.571754608952619 1.30457840948769 0.302675318800333 1.14007511822789 2.200676495613032 0.02502240784087517 1.402893483424553 2.061369194121731 2.430616761925077 1.019909066034401 1.392911511673118 1.880834072596762 1.65887993618586 1.854477465005829 1.581964280941516 0.4524338473246255 1.055389551643304 1.698333284152199 0.09466761685818237 1.682866360474479 0.9203714703780533 0.7883175376961334 0.7735400964254449 3.412029603373202 0.9389264766934377 -0.2286519426867513 1.144899531664923 1.959607319634346 0.6820149079212618 0.3659211651228045 1.24457081853047 1.240565634184464 -0.808996907369538 -0.4095289441050487 -0.03984319060471475 1.407376888599285 2.50506971416251 0.94866459792695 1.377655556676116 0.8621283462953613 1.45572695894971 1.337481186615599 -0.03581875865121709 0.6192721038839956 3.077072925413614 1.175569494850847 0.5891259026809106 0.8610329904326544 2.435659738681764 -0.6747189852569966 2.031022199948546 0.5566425448896293 1.447685798202472 -0.02473867507343686 2.007858891101517 1.263674437306774 1.713629327117102 1.254938386649121 1.892944756136862 0.4684330377337796 -0.4027056816821775 1.11704801106262 -1.414532689933904 0.9485616889624722 2.358307923865745 0.04632719508016703 -0.6085524733548022 -0.3156332989897632 -0.3994652220730563 -2.026887037281629 -0.7817493723936182 0.0836673467640805 -0.3829577698394881 1.29819866462604 1.084558245499732 0.9864835029701788 1.119071535438783 -0.1405294307137892 2.214607285013491 -0.006327316715517561 2.212590467370691 1.26588548344322 -0.5429082556043867 0.09490211216100741 2.549335338093492 0.8121250547881319 2.174788321516417 0.6011032880870217 2.834863443688122 1.019969423481439 1.189008892235732 3.880782990535681 -0.4613835224709339 1.30080684608793 0.3341630366858428 2.558988066567454 -0.3541816158487625 1.992892726585691 1.57301742825892 0.7891172056259123 1.335270370293899 1.099317126413363 1.94850228775292 0.7180743145007634 0.2732012480611574 1.88701768090489 0.919039111854304 0.5164906832658123 -1.035409137727921 1.531827498628754 0.8721386145599138 2.057111508572006 0.5829095702425009 1.219918106522843 1.910275635648726 -0.6316897441802187 -0.5740899163490216 -0.4301836663620986 1.433373775047487 0.1971735647599688 1.101488459356943 1.066928958382982 0.6455351033321124 -0.546989207700902 2.366065132501835 3.587260717546521 1.157362400953183 2.326308940608198 3.695299542189807 1.718391676775385 -0.5572607835669456 1.774638640170457 0.6427803541090709 0.8179656256527835 -1.427179248433601 2.347050988513118 0.6939815297602763 1.225605599505495 -0.9044219759169925 -1.142731033609628 1.418193974248724 1.178306680323816 -1.201815224697584 0.04591899739153427 -0.02887756889841087 2.140995206014014 1.095581224518306 -0.2437858265225779 -0.5795306559652984 -0.5941933636597605 1.547264303474973 -1.642462107162598 0.6960281312476969 -0.008382857825509848 -0.725508676535545 1.51600770180333 1.539926920247545 1.388250699055266 2.224010740525287 1.649318957661719 -0.6147625607082601 0.5515733586034677 0.558258455676736 0.7189388536036069 0.5158843706883898 2.153889332982802 0.3828076118746291 0.3677019184304079 1.659632545441951 -0.8225887804579344 1.748981685320509 1.426549705601257 1.406016845030578 0.8899866163193882 0.6809714511623861 2.427223674550403 3.466490892897378 1.342136503754375 0.06841851946696453 1.373415909361131 0.7652623797518443 0.8396093766790385 1.938274360789231 1.904246521995326 -0.03650121825884889 0.1027546072033909 1.566785952884902 3.06898052727411 1.214927572139952 1.050208409427337 2.976858768261151 0.8105352839625762 0.805432928329994 0.9354516434670886 -1.487850555829779 0.5383045120858965 1.236632524549359 2.061736620589021 1.234899946130458 1.795316488842313 0.6678628022755861 -0.4519810377511588 -0.02987276372091996 1.375123105925975 1.315580464196835 1.022030750123769 0.09034629544573369 1.469521469802323 1.748079437898673 0.4963024850608647 0.7332648207644714 1.335852841324077 0.6439092622804476 0.6781507496429851 0.6543782683697399 1.668742462695305 -0.04875535039821921 -1.214398420486452 2.354469182113136 0.03891042224179864 0.6134475222352157 0.4336358043279851 1.326013407781256 -0.117168836983129 1.759926434163588 0.4791213767042255 0.4365430457578929 0.5690695258146682 -0.0531786578714506 2.367319974594577 0.4491257894231035 1.385175892949472 0.7138088403476408 0.2064240103282239 0.2015159906384059 0.5934245273683694 -0.6713871540982557 1.027224205460117 2.552271735323127 1.65150506372262 1.904433424822259 -0.5620851860805303 0.06669369871967024 -0.819153914633457 0.7998890249813373 1.20972782458538 -0.6143752935638263 1.224218276521915 1.52530528404342 1.693008208210148 -0.7709484090466059 1.598971350634212 2.244904936766772 0.3058929035697816 1.227014894672552 1.12644227640822 1.667369525656748 -0.7886982200436146 0.5671573707196975 0.2960941368496872 0.9718998931972267 1.1431859413311 -0.6172018382535209 2.131221940354825 0.8101330855984981 0.2044766135092729 2.920543355522538 1.728484134270207 -0.3428597341096586 1.867400758332775 2.142384246429031 1.39864944981526 1.782769914485578 1.48733335784685 1.069399026465054 0.8236770291736946 1.895725315388256 1.183307434295155 1.687756226992983 1.816280365106959 1.541983281114427 0.674184937519053 1.155707182973492 0.8896722420006936 1.348739935468164 0.07945443115763218 0.4015850998285639 0.2940033118937194 0.8144776439829758 1.010277227746017 1.685498799452889 -0.838982275528326 -1.18018020407638 -0.868412444924779 0.7400610340786427 1.707052905550757 1.919288250721198 1.251273405114258 0.4644987613273633 1.334692331775329 2.538605205036323 2.485722427011433 0.856324109447091 -0.7136542582953915 2.825465101313566 0.08694272044256701 -0.4285086279494847 -0.2502875438659184 2.437339001762728 3.023974046395975 0.8106101645444662 0.6521151035084938 0.09644813760464621 0.06353816989227168 2.161572265357913 -0.9281891337178096 0.8322866265625576 -0.6624871844665274 2.055961578390782 0.01587797975935534 0.7997631265641372 0.4317059172430303 1.531143804662626 2.562595209034343 0.1901457709002661 0.5061057646508257 0.5780694614211119 -0.5401165798621181 -0.8177535916998482 1.866859409187512 0.6109486305688941 1.11376363312526 0.03183485263470687 0.8165023556504525 0.8693515787128633 1.734754815439754 0.1004730547138618 0.6178010884581566 0.4824611720210606 1.182049877268903 -1.018185671997708 1.585419468001058 1.357963678117527 1.397876718473706 -0.2619430180440445 1.005324640919879 0.7964759685407805 2.099693653056788 -0.3755007476680732 0.8312975399984183 1.396105881808023 -0.1896598582949165 2.002961271703618 3.583366494426321 1.330851465240229 1.130103694295259 0.7688013635923769 -0.1286103296916599 -0.1434832589799135 1.643611062741993 1.121930289267041 0.7369263085398459 -0.399645103166467 0.574005355245518 0.5202189972705118 2.985653630789157 2.072293741185865 1.07616339901463 -0.1765205754710266 -0.7930859474858163 1.157926218550625 0.7843350672317007 0.4196496832593138 -0.07300266330033933 1.438911595136333 2.064411160627518 2.726359221689226 2.54805504911209 1.100018848512457 1.043593019354887 1.653632122087634 1.080919973623599 3.385779270214321 1.683797773663598 0.9532190916410197 2.40329493259908 1.283772772485365 1.67695018692607 1.787999220336419 0.7531209642719869 1.65533954344822 0.3041272538354668 2.703140157448808 2.865577528054054 1.083660166482454 3.183503862026781 0.1176591644524645 0.2033963524096726 0.7949491847411757 0.923829221952576 1.324175554457839 -0.8785271277470366 -0.2022031237603816 1.731298686205138 0.9075127909346693 1.47241382695859 2.765055927301861 1.662551233574523 2.229702341089688 2.420339239959411 0.01856912204413241 2.379803526212374 1.194682883418982 0.7314613660992719 2.619104084085683 1.194266390038433 1.755229084507919 -0.7717968320741864 2.775922639288721 3.426301255877113 2.261015872096904 0.651964538515456 1.345939947315443 0.2105043028449838 0.3328698838274554 3.642706438211579 -0.3904994645123505 0.9844749844701601 0.2924099967396088 0.06560500662064384 -0.06615896167074142 -0.7604495935580171 0.4200986436670405 1.442699027263202 0.4610606186888362 1.921740459008894 2.176654783825957 0.2766904269817657 0.7407258065608351 0.1817407648628061 0.9644631817009489 0.6065483491804464 0.2324125184457826 1.274911601852321 -0.8986056273266818 1.404120595961857 0.04156799293478086 2.287418670386505 2.05840669246076 1.799679819125657 2.004151444242522 0.6644434702942055 -0.1966604233793228 0.8685309527545859 0.8165934975554585 1.953237274081445 2.042743234387201 0.2318316739303269 1.288867345651097 0.9215959041223284 -0.07665186592476791 -0.8030549321239862 0.3653138225572382 1.381729284002746 2.190014342076591 0.688056952439958 0.5071921430901954 -0.6041002388472794 1.459369139552627 -0.8416502124388883 1.001052745331846 1.596638837915149 0.5348346491908362 0.2901786382189061 1.342720144385798 1.589092844018214 -0.2283721036511688 0.8305215735059919 2.165485152240601 -0.03252788961942565 0.02328071227927819 1.746912156602784 2.051175246634993 1.914325258851608 -0.2794919436484127 1.426645034300429 0.6929458022750473 1.394410764453883 1.284846736041881 1.090896414615004 -0.7625203833382994 3.081651907582955 0.7665151300770845 0.6821018419088138 0.5603338674426759 0.6572856151612331 1.392476589756923 -0.04647067345709188 1.08147339871604 0.7078394909697348 0.01086012000244907 0.1384793254931365 2.043919155687828 0.4134384610983507 -0.7080618667851521 0.2352748837640597 1.812832124001387 -0.1268862569762619 1.685355392529269 -0.144857368553704 0.8889031624158457 1.949061940584198 1.768733336991409 0.7623401939543458 0.767937305445493 0.6550831207704331 0.3362298871685868 0.7764073163080257 2.359765159287158 1.429089373463981 1.309205351542116 -0.3268723432746896 0.8104497326422911 2.956965849238324 0.358524762718711 1.589412679232686 2.615427889946867 1.610693624201783 1.55191690831084 2.809198581573606 1.330755065953391 -0.1283170853071223 -0.474750683055517 0.9803923739926476 1.346235159202204 1.392911694970242 0.6346043347252104 2.370365713880602 2.698217365974108 0.001521859597514053 0.7471014323465719 2.117297049631915 2.928504152930254 1.908277373323209 1.500462929674708 1.133317338887011 1.507896432518121 2.807132121248947 2.374119900294267 0.907002995671104 0.6909785110621898 -0.6472724969890027 1.797647281265343 0.9467178409047315 -0.6648896706502492 1.530229712113611 0.7206469850455717 -0.5946598393887041 1.456691133146709 2.248258506550396 0.170271558722038 1.831513988086103 -0.3025803355234244 2.385418663428565 1.146733367404714 2.170143805859415 0.4245995483641243 0.949917907426168 0.3102974373794284 1.715536410523382 3.290316642569182 1.242702262755937 -0.297971103458347 1.369721796531666 -0.1401078494450936 1.074358209119162 -0.248762374937652 0.7741944473547974 0.514507383681623 1.787839763152484 0.4533666972649788 0.3519589345259729 2.036886889178933 0.9996729504490831 2.289381718120088 0.431587546009627 -0.5237002243209679 1.239567952522294 0.6761313996158718 0.9071710058567745 0.8837286842961545 1.035352586794545 2.144365426158453 0.4783708236729683 2.033672969617412 1.801519424099742 -0.8344022765758961 0.5602884656419493 0.72107370598711 -0.1512500608802831 -0.8109086017955673 1.763530583698427 -1.176822767358457 -0.4942355480811056 1.161308989334149 1.222869381750991 2.659789397277532 1.425338756687175 2.404180720178092 0.2384537595927544 2.212845748940844 0.5550363648913921 1.524028774206119 0.5981836136185896 0.5294540428864375 0.8696654600401873 0.7833674134167391 1.508798572309083 0.2534079338862683 0.66946204116125 0.1919327690165185 2.137979101543531 1.811537552181822 2.211414480126396 -0.3667803731405608 1.898889228322251 1.622418177467973 1.947201617655731 1.481833835653183 0.601764812275871 0.6303918701726243 1.51994258186477 -0.1161953258017661 2.024617398746776 1.422339904094752 1.183154710585878 0.9304865918550239 1.071618644818221 0.7033879178281861 1.568412859625416 -0.9750664306465477 0.7493605057505957 1.668917742417485 1.767266362343072 1.722993495766987 1.387423726793022 -0.2828739000600107 0.8299562156869035 1.21918696434744 1.01931148214656 1.310285271474827 -0.9664421259477922 1.024933466312315 1.422898978159024 1.566950347807346 0.9021197631092972 0.6222499071420458 0.1321155481729535 1.288063095180657 1.352656771706024 0.1895061657905364 0.1920373324287531 1.478312968404502 0.6999690645255463 1.718244503797262 1.693025225446638 -0.07355066568929547 1.324113730690591 -0.8910892301905946 0.08014974002211994 1.374078902490329 2.627495388780011 -1.307638481303528 -0.3934299458205341 1.611516609206224 0.9667655239515172 0.597723163117555 0.4400589500081084 0.7293240539552928 1.285119876893175 2.336139525648818 -0.6673365593984886 1.530074485490484 0.2171340933522691 -0.0446354765064918 -0.1038020661746528 -0.1819091484153785 0.8678733854163624 1.232817380535235 0.9167791661759996 1.168013108607712 1.493390515296723 1.536417437632474 0.6797286763461281 0.9035159224287149 -0.8191069059672471 0.6012359607464048 1.344444829105806 0.8183485917139273 1.614778215988197 2.223568213372829 0.3734889189821574 1.422649302463133 0.9888487551982762 -1.496709144298501 0.8646479649906167 1.160002618343647 0.6743586761693734 0.6734460402057004 -0.3863713742678794 -0.9604983793458164 2.461399454123412 1.003077847052847 0.2592339470695195 1.569928793084649 0.6382274482858993 -0.9429592775201137 -0.4581566653767766 0.4315195019696789 2.195549098550692 1.235049479489091 2.840815442383532 0.9410390060244023 2.256518702798669 1.234668114148285 0.9465881972606696 -0.5457540293086707 0.162386254249927 0.534006432843394 2.233921297718268 -0.08698581368483005 1.674617042026495 -0.1836532309896213 3.14122995765425 1.433464763688435 -0.206909665367156 1.444601979938507 1.503067169865708 0.08666750927248379 0.03915375698533374 1.837347527702306 1.064527516722064 2.244462190723304 1.617210464548076 1.267451644905893 0.9200161026441995
AutoHotkey
contributed by Laszlo on the ahk forum <lang AutoHotkey>Loop 40
R .= RandN(1,0.5) "`n" ; mean = 1.0, standard deviation = 0.5
MsgBox %R%
RandN(m,s) { ; Normally distributed random numbers of mean = m, std.dev = s by Box-Muller method
Static i, Y If (i := !i) { ; every other call Random U, 0, 1.0 Random V, 0, 6.2831853071795862 U := sqrt(-2*ln(U))*s Y := m + U*sin(V) Return m + U*cos(V) } Return Y
}</lang>
Avail
<lang Avail>Method "U(_,_)" is [ lower : number, upper : number | divisor ::= ((1<<32)) ÷ (upper - lower)→double; map a pRNG through [i : integer | (i ÷ divisor) + lower] ];
Method "a Marsaglia polar sampler" is [ generator for [ yield : [double]→⊤ | source ::= U(-1, 1); Repeat [ x ::= take 1 from source[1]; y ::= take 1 from source[1]; s ::= x^2 + y^2; If 0 < s < 1 then [ factor ::= ((-2 × ln s) ÷ s) ^ 0.5; yield(x × factor); yield(y × factor); ]; ] ] ];
// the default distribution has mean 0 and std dev 1.0, so we scale the values sampler ::= map a Marsaglia polar sampler through [d : double | d ÷ 2.0 + 1.0]; values ::= take 1000 from sampler;</lang>
AWK
One-liner: <lang awk>$ awk 'func r(){return sqrt(-2*log(rand()))*cos(6.2831853*rand())}BEGIN{for(i=0;i<1000;i++)s=s" "1+0.5*r();print s}'</lang>
Readable version: <lang awk> function r() {
return sqrt( -2*log( rand() ) ) * cos(6.2831853*rand() )
}
BEGIN {
n=1000 for(i=0;i<n;i++) { x = 1 + 0.5*r() s = s" "x } print s
} </lang>
- Output:
first few values only
0.783753 1.16682 1.17989 1.14975 1.34784 0.29296 0.979227 1.04402 0.567835 1.58812 0.465559 1.27186 0.324533 0.725827 -0.0626549 0.632273 1.0145 1.3387 0.861667 1.04147 1.2576 1.02497 0.58453 0.9619 1.26902 0.851048 -0.126259 0.863256
...
BASIC
RANDOMIZE TIMER 'seeds random number generator with the system time pi = 3.141592653589793# DIM a(1 TO 1000) AS DOUBLE CLS FOR i = 1 TO 1000 a(i) = 1 + SQR(-2 * LOG(RND)) * COS(2 * pi * RND) NEXT i
BBC BASIC
<lang bbcbasic> DIM array(999)
FOR number% = 0 TO 999 array(number%) = 1.0 + 0.5 * SQR(-2*LN(RND(1))) * COS(2*PI*RND(1)) NEXT mean = SUM(array()) / (DIM(array(),1) + 1) array() -= mean stdev = MOD(array()) / SQR(DIM(array(),1) + 1) PRINT "Mean = " ; mean PRINT "Standard deviation = " ; stdev</lang>
- Output:
Mean = 1.01848064 Standard deviation = 0.503551814
C
<lang c>#include <stdlib.h>
- include <math.h>
- ifndef M_PI
- define M_PI 3.14159265358979323846
- endif
double drand() /* uniform distribution, (0..1] */ {
return (rand()+1.0)/(RAND_MAX+1.0);
} double random_normal() /* normal distribution, centered on 0, std dev 1 */ {
return sqrt(-2*log(drand())) * cos(2*M_PI*drand());
} int main() {
int i; double rands[1000]; for (i=0; i<1000; i++) rands[i] = 1.0 + 0.5*random_normal(); return 0;
}</lang>
C#
<lang csharp> private static double randomNormal() { return Math.Cos(2 * Math.PI * tRand.NextDouble()) * Math.Sqrt(-2 * Math.Log(tRand.NextDouble())); } </lang>
Then the methods in Random numbers#Metafont are used to calculate the average and the Standard Deviation: <lang csharp> static Random tRand = new Random();
static void Main(string[] args) { double[] a = new double[1000];
double tAvg = 0; for (int x = 0; x < a.Length; x++) { a[x] = randomNormal() / 2 + 1; tAvg += a[x]; }
tAvg /= a.Length; Console.WriteLine("Average: " + tAvg.ToString());
double s = 0; for (int x = 0; x < a.Length; x++) { s += Math.Pow((a[x] - tAvg), 2); } s = Math.Sqrt(s / 1000);
Console.WriteLine("Standard Deviation: " + s.ToString());
Console.ReadLine(); } </lang>
An example result:
Average: 1,00510073053613 Standard Deviation: 0,502540443430955
C++
The new C++ standard looks very similar to the Boost library example below.
<lang cpp>#include <random>
- include <functional>
- include <vector>
- include <algorithm>
using namespace std;
int main() {
random_device seed; mt19937 engine(seed()); normal_distribution<double> dist(1.0, 0.5); auto rnd = bind(dist, engine);
vector<double> v(1000); generate(v.begin(), v.end(), rnd); return 0;
}</lang>
<lang cpp>#include <cstdlib> // for rand
- include <cmath> // for atan, sqrt, log, cos
- include <algorithm> // for generate_n
double const pi = 4*std::atan(1.0);
// simple functor for normal distribution class normal_distribution { public:
normal_distribution(double m, double s): mu(m), sigma(s) {} double operator() const // returns a single normally distributed number { double r1 = (std::rand() + 1.0)/(RAND_MAX + 1.0); // gives equal distribution in (0, 1] double r2 = (std::rand() + 1.0)/(RAND_MAX + 1.0); return mu + sigma * std::sqrt(-2*std::log(r1))*std::cos(2*pi*r2); }
private:
const double mu, sigma;
};
int main() {
double array[1000]; std::generate_n(array, 1000, normal_distribution(1.0, 0.5)); return 0;
}</lang>
This example used Mersenne Twister generator. It can be changed by changing the typedef.
<lang cpp>
- include <vector>
- include "boost/random.hpp"
- include "boost/generator_iterator.hpp"
- include <boost/random/normal_distribution.hpp>
- include <algorithm>
typedef boost::mt19937 RNGType; ///< mersenne twister generator
int main() {
RNGType rng; boost::normal_distribution<> rdist(1.0,0.5); /**< normal distribution with mean of 1.0 and standard deviation of 0.5 */
boost::variate_generator< RNGType, boost::normal_distribution<> > get_rand(rng, rdist);
std::vector<double> v(1000); generate(v.begin(),v.end(),get_rand); return 0;
} </lang>
Clojure
<lang lisp>(import '(java.util Random)) (def normals
(let [r (Random.)] (take 1000 (repeatedly #(-> r .nextGaussian (* 0.5) (+ 1.0))))))</lang>
Common Lisp
<lang lisp>(loop for i from 1 to 1000
collect (1+ (* (sqrt (* -2 (log (random 1.0)))) (cos (* 2 pi (random 1.0))) 0.5)))</lang>
Crystal
<lang ruby>n, mean, sd, tau = 1000, 1, 0.5, (2 * Math::PI) array = Array.new(n) { mean + sd * Math.sqrt(-2 * Math.log(rand)) * Math.cos(tau * rand) }
mean = array.sum / array.size standev = Math.sqrt( array.sum{ |x| (x - mean) ** 2 } / array.size ) puts "mean = #{mean}, standard deviation = #{standev}"</lang>
- Output:
mean = 1.0093442539237896, standard deviation = 0.504694489463623
D
<lang d>import std.stdio, std.random, std.math;
struct NormalRandom {
double mean, stdDev;
// Necessary because it also defines an opCall. this(in double mean_, in double stdDev_) pure nothrow { this.mean = mean_; this.stdDev = stdDev_; }
double opCall() const /*nothrow*/ { immutable r1 = uniform01, r2 = uniform01; // Not nothrow. return mean + stdDev * sqrt(-2 * r1.log) * cos(2 * PI * r2); }
}
void main() {
double[1000] array; auto nRnd = NormalRandom(1.0, 0.5); foreach (ref x; array) //x = nRnd; x = nRnd();
}</lang>
Alternative Version
(Untested)
<lang d>import tango.math.random.Random;
void main() {
double[1000] list; auto r = new Random(); foreach (ref l; list) { r.normalSource!(double)()(l); l = 1.0 + 0.5 * l; }
}</lang>
Delphi
Delphi has RandG function which generates random numbers with normal distribution using Marsaglia-Bray algorithm:
<lang Delphi>program Randoms;
{$APPTYPE CONSOLE}
uses
Math;
var
Values: array[0..999] of Double; I: Integer;
begin // Randomize; Commented to obtain reproducible results
for I:= Low(Values) to High(Values) do Values[I]:= RandG(1.0, 0.5); // Mean = 1.0, StdDev = 0.5 Writeln('Mean = ', Mean(Values):6:4); Writeln('Std Deviation = ', StdDev(Values):6:4); Readln;
end.</lang>
- Output:
Mean = 1.0098 Std deviation = 0.5016
DWScript
<lang delphi>var values : array [0..999] of Float; var i : Integer;
for i := values.Low to values.High do
values[i] := RandG(1, 0.5);</lang>
E
<lang e>accum [] for _ in 1..1000 { _.with(entropy.nextGaussian()) }</lang>
EasyLang
<lang>for i range 1000
a[] &= 1 + 0.5 * sqrt (-2 * logn randomf) * cos (360 * randomf)
. print a[]</lang>
Eiffel
<lang eiffel> class APPLICATION
inherit ARGUMENTS
create make
feature {NONE} -- Initialization
l_time: TIME l_seed: INTEGER math:DOUBLE_MATH rnd:RANDOM Size:INTEGER once Result:= 1000 end
make -- Run application. local ergebnis:ARRAY[DOUBLE] tavg: DOUBLE x: INTEGER tmp: DOUBLE text : STRING
do -- initialize random generator create l_time.make_now
l_seed := l_time.hour l_seed := l_seed * 60 + l_time.minute l_seed := l_seed * 60 + l_time.second l_seed := l_seed * 1000 + l_time.milli_second create rnd.set_seed (l_seed)
-- initialize random number container and math create ergebnis.make_filled (0.0, 1, size) tavg := 0; create math
from x := 1 until x > ergebnis.count loop tmp := randomNormal / 2 + 1 tavg := tavg + tmp ergebnis.enter (tmp , x) x := x + 1 end
tavg := tavg / ergebnis.count text := "Average: " text.append_double (tavg) text.append ("%N") print(text)
tmp := 0 from x:= 1 until x > ergebnis.count loop tmp := tmp + (ergebnis.item (x) - tavg)^2 x := x + 1 end
tmp := math.sqrt (tmp / ergebnis.count) text := "Standard Deviation: " text.append_double (tmp) text.append ("%N") print(text)
end
randomNormal:DOUBLE
local
first: DOUBLE second: DOUBLE
do
rnd.forth first := rnd.double_item rnd.forth second := rnd.double_item
Result := math.cosine (2 * math.pi * first) * math.sqrt (-2 * math.log (second))
end end </lang>
Example Result
Average: 1.0079398405028137 Standard Deviation: 0.49042787564453988
Elena
ELENA 4.1 : <lang elena>import extensions; import extensions'math;
randomNormal() {
^ cos(2 * Pi_value * randomGenerator.nextReal()) * sqrt(-2 * ln(randomGenerator.nextReal()))
}
public program() {
real[] a := new real[](1000); real tAvg := 0; for (int x := 0, x < a.Length, x += 1) { a[x] := (randomNormal()) / 2 + 1; tAvg += a[x] }; tAvg /= a.Length; console.printLine("Average: ", tAvg); real s := 0; for (int x := 0, x < a.Length, x += 1) { s += power(a[x] - tAvg, 2) }; s := sqrt(s / 1000); console.printLine("Standard Deviation: ", s); console.readChar()
}</lang>
- Output:
Average: 0.9842420481571 Standard Deviation: 0.5109070975558
Elixir
<lang elixir>defmodule Random do
def normal(mean, sd) do {a, b} = {:rand.uniform, :rand.uniform} mean + sd * (:math.sqrt(-2 * :math.log(a)) * :math.cos(2 * :math.pi * b)) end
end
std_dev = fn (list) ->
mean = Enum.sum(list) / length(list) sd = Enum.reduce(list, 0, fn x,acc -> acc + (x-mean)*(x-mean) end) / length(list) |> :math.sqrt IO.puts "Mean: #{mean},\tStdDev: #{sd}" end
xs = for _ <- 1..1000, do: Random.normal(1.0, 0.5) std_dev.(xs)</lang>
- Output:
Mean: 1.009079383094275, StdDev: 0.4991894476975088
used Erlang function :rand.normal
<lang elixir>xs = for _ <- 1..1000, do: 1.0 + :rand.normal * 0.5
std_dev.(xs)</lang>
- Output:
Mean: 0.9955701150615597, StdDev: 0.5036412260426065
Erlang
<lang erlang> mean(Values) ->
mean(tl(Values), hd(Values), 1).
mean([], Acc, Length) ->
Acc / Length;
mean(Values, Acc, Length) ->
mean(tl(Values), hd(Values)+Acc, Length+1).
variance(Values) ->
Mean = mean(Values), variance(Values, Mean, 0) / length(Values).
variance([], _, Acc) ->
Acc;
variance(Values, Mean, Acc) ->
Diff = hd(Values) - Mean, DiffSqr = Diff * Diff, variance(tl(Values), Mean, Acc + DiffSqr).
stddev(Values) ->
math:sqrt(variance(Values)).
normal(Mean, StdDev) ->
U = random:uniform(), V = random:uniform(), Mean + StdDev * ( math:sqrt(-2 * math:log(U)) * math:cos(2 * math:pi() * V) ). % Erlang's math:log is the natural logarithm.
main(_) ->
X = [ normal(1.0, 0.5) || _ <- lists:seq(1, 1000) ], io:format("mean = ~w\n", [mean(X)]), io:format("stddev = ~w\n", [stddev(X)]).
</lang>
- Output:
mean = 1.0118289913718608 stddev = 0.5021636849524854
ERRE
<lang> PROGRAM DISTRIBUTION
! ! for rosettacode.org !
! formulas taken from TI-59 Master Library manual
CONST NUM_ITEM=1000
!VAR SUMX#,SUMX2#,R1#,R2#,Z#,I%
DIM A#[1000]
BEGIN ! seeds random number generator with system time
RANDOMIZE(TIMER)
PRINT(CHR$(12);) !CLS SUMX#=0 SUMX2#=0
FOR I%=1 TO NUM_ITEM DO R1#=RND(1) R2#=RND(1) Z#=SQR(-2*LOG(R1#))*COS(2*π*R2#) A#[I%]=Z#/2+1 ! I want a normal distribution with ! mean=1 and std.dev=0.5 SUMX#+=A#[I%] SUMX2#+=A#[I%]*A#[I%] END FOR
Z#=SUMX#/NUM_ITEM
PRINT("Average is";Z#) PRINT("Standard dev. is";SQR(SUMX2#/NUM_ITEM-Z#*Z#))
END PROGRAM </lang>
Euler Math Toolbox
<lang Euler Math Toolbox> >v=normal(1,1000)*0.5+1; >mean(v), dev(v)
1.00291801071 0.498226876528
</lang>
Euphoria
<lang euphoria>include misc.e
function RandomNormal()
atom x1, x2 x1 = rand(999999) / 1000000 x2 = rand(999999) / 1000000 return sqrt(-2*log(x1)) * cos(2*PI*x2)
end function
constant n = 1000 sequence s s = repeat(0,n) for i = 1 to n do
s[i] = 1 + 0.5 * RandomNormal()
end for</lang>
F#
<lang fsharp> let n = MathNet.Numerics.Distributions.Normal(1.0,0.5) List.init 1000 (fun _->n.Sample()) </lang>
- Output:
[0.734433576; 1.54225304; 0.4407528678; 1.177675412; 0.4318617021; 0.6026656337; 0.769764924; 1.104693934; 0.6297500925; 0.9594598077; 1.684736389; 1.160376323; 0.883354356; 0.9513968363; 0.9727698268; 0.5315570949; 0.9599239266; 1.564976755; 0.7232002879; 1.084139442; 1.220914517; 0.3553085946; 1.112549824; 1.989443553; 0.5752307543; 1.156682549; 0.7886670467; 0.02050745923; 1.532060208; 1.18789591; 1.408946777; 1.038812004; 1.724679503; 1.671565045; 1.266831442; 1.363611654; 1.705819067; 0.5772366328; 0.4503488498; 1.496891481; 0.9831877282; 0.3845460366; 0.8253240671; 1.298969969; 0.4265904553; 0.9303696876; 0.445003361; 0.753175816; 0.6143534043; 1.059982235; 0.7143206784; 0.2233328038; 1.005178481; 0.7697392436; 0.5904948577; 0.5127953044; 0.6467346747; 0.7929387604; -0.1501790761; 0.8750780903; 0.941704369; 1.37941579; 0.4739006145; 1.998886344; 1.219428519; 0.06270791476; 1.097739804; 0.7584232803; 1.042177217; 1.166561247; 1.502357164; 1.171525776; 0.1528807432; 0.2289389756; 1.36208422; 0.3714421124; 1.299571092; 1.171553369; 1.317807265; 1.616662281; 1.724223246; 1.059580642; 1.270520918; -0.1827677907; 1.938593232; 1.420362143; 1.888357595; 0.7851629936; 0.7080554899; 0.7747215818; 1.403719877; 0.5765950249; 1.275206565; 0.6292054813; 1.525562798; 0.6224640457; 0.8524078517; 0.7646595627; 0.6799834691; 0.773111053; ...]
==F#==
<lang fsharp>let gaussianRand count =
let o = new System.Random() let pi = System.Math.PI let gaussrnd = (fun _ -> 1. + 0.5 * sqrt(-2. * log(o.NextDouble())) * cos(2. * pi * o.NextDouble())) [ for i in {0 .. (int count)} -> gaussrnd() ]</lang>
Factor
<lang factor>USING: random ; 1000 [ 1.0 0.5 normal-random-float ] replicate</lang>
Falcon
<lang falcon>a = [] for i in [0:1000] : a+= norm_rand_num()
function norm_rand_num()
pi = 2*acos(0) return 1 + (cos(2 * pi * random()) * pow(-2 * log(random()) ,1/2)) /2
end</lang>
Fantom
Two solutions. The first uses Fantom's random-number generator, which produces a uniform distribution. So, convert to a normal distribution using a formula:
<lang fantom> class Main {
static const Float PI := 0.0f.acos * 2 // we need to precompute PI
static Float randomNormal () { return (Float.random * PI * 2).cos * (Float.random.log * -2).sqrt }
public static Void main () { mean := 1.0f sd := 0.5f Float[] values := [,] // this is the collection to fill with random numbers 1000.times { values.add (randomNormal * sd + mean) } }
} </lang>
The second calls out to Java's Gaussian random-number generator:
<lang fantom> using [java] java.util::Random
class Main {
Random generator := Random()
Float randomNormal () { return generator.nextGaussian }
public static Void main () { rnd := Main() // create an instance of Main class, which holds the generator mean := 1.0f sd := 0.5f Float[] values := [,] // this is the collection to fill with random numbers 1000.times { values.add (rnd.randomNormal * sd + mean) } }
} </lang>
Forth
<lang forth>require random.fs here to seed
-1. 1 rshift 2constant MAX-D \ or s" MAX-D" ENVIRONMENT? drop
- frnd ( -- f ) \ uniform distribution 0..1
rnd rnd dabs d>f MAX-D d>f f/ ;
- frnd-normal ( -- f ) \ centered on 0, std dev 1
frnd pi f* 2e f* fcos frnd fln -2e f* fsqrt f* ;
- ,normals ( n -- ) \ store many, centered on 1, std dev 0.5
0 do frnd-normal 0.5e f* 1e f+ f, loop ;
create rnd-array 1000 ,normals</lang>
For newer versions of gforth (tested on 0.7.3), it seems you need to use HERE SEED ! instead of HERE TO SEED, because SEED has been made a variable instead of a value.
<lang>rnd rnd dabs d>f</lang> is necessary, but surprising and definitely not well documented / perhaps not compliant.
Fortran
<lang fortran>PROGRAM Random
INTEGER, PARAMETER :: n = 1000 INTEGER :: i REAL :: array(n), pi, temp, mean = 1.0, sd = 0.5
pi = 4.0*ATAN(1.0) CALL RANDOM_NUMBER(array) ! Uniform distribution
! Now convert to normal distribution
DO i = 1, n-1, 2 temp = sd * SQRT(-2.0*LOG(array(i))) * COS(2*pi*array(i+1)) + mean array(i+1) = sd * SQRT(-2.0*LOG(array(i))) * SIN(2*pi*array(i+1)) + mean array(i) = temp END DO
! Check mean and standard deviation
mean = SUM(array)/n sd = SQRT(SUM((array - mean)**2)/n) WRITE(*, "(A,F8.6)") "Mean = ", mean WRITE(*, "(A,F8.6)") "Standard Deviation = ", sd
END PROGRAM Random</lang>
- Output:
Mean = 0.995112 Standard Deviation = 0.503373
Free Pascal
Free Pascal provides the randg function in the RTL math unit that produces Gaussian-distributed random numbers with the Box-Müller algorithm.
<lang pascal> function randg(mean,stddev: float): float; </lang>
FreeBASIC
<lang freebasic>' FB 1.05.0 Win64
Const pi As Double = 3.141592653589793 Randomize
' Generates normally distributed random numbers with mean 0 and standard deviation 1 Function randomNormal() As Double
Return Cos(2.0 * pi * Rnd) * Sqr(-2.0 * Log(Rnd))
End Function
Dim r(0 To 999) As Double Dim sum As Double = 0.0
' Generate 1000 normally distributed random numbers ' with mean 1 and standard deviation 0.5 ' and calculate their sum For i As Integer = 0 To 999
r(i) = 1.0 + randomNormal/2.0 sum += r(i)
Next
Dim mean As Double = sum / 1000.0
Dim sd As Double sum = 0.0 ' Now calculate their standard deviation For i As Integer = 0 To 999
sum += (r(i) - mean) ^ 2.0
Next sd = Sqr(sum/1000.0)
Print "Mean is "; mean Print "Standard Deviation is"; sd Print Print "Press any key to quit" Sleep</lang> Sample result:
- Output:
Mean is 1.000763573902885 Standard Deviation is 0.500653063426955
Frink
<lang frink> a = new array for i = 1 to 1000
a.push[randomGaussian[1, 0.5]]
</lang>
FutureBasic
Note: To generate the random number, rather than using FB's native "rnd" function, this code wraps C code into the RandomZeroToOne function. <lang futurebasic> include "ConsoleWindow"
local fn RandomZeroToOne as double dim as double result BeginCCode
result = (double)( (rand() % 100000 ) * 0.00001 );
EndC end fn = result
local fn RandomGaussian as double dim as double r
r = fn RandomZeroToOne end fn = 1 + .5 * ( sqr( -2 * log(r) ) * cos( 2 * pi * r ) )
dim as long i dim as double mean, std, a(1000)
for i = 1 to 1000
a(i) = fn RandomGaussian mean += a(i)
next mean = mean / 1000
for i = 1 to 1000
std += ( a(i) - mean )^2
next std = std / 1000
print " Average:"; mean print "Standard Deviation:"; std </lang> Output:
Average: 1.0258434498 Standard Deviation: 0.2771047023
Go
This solution uses math/rand package in the standard library. See also though the subrepository rand package at https://godoc.org/golang.org/x/exp/rand, which also has a NormFloat64 and has a rand source with a number of advantages over the one in standard library. <lang go>package main
import (
"fmt" "math" "math/rand" "strings" "time"
)
const mean = 1.0 const stdv = .5 const n = 1000
func main() {
var list [n]float64 rand.Seed(time.Now().UnixNano()) for i := range list { list[i] = mean + stdv*rand.NormFloat64() } // show computed mean and stdv of list var s, sq float64 for _, v := range list { s += v } cm := s / n for _, v := range list { d := v - cm sq += d * d } fmt.Printf("mean %.3f, stdv %.3f\n", cm, math.Sqrt(sq/(n-1))) // show histogram by hdiv divisions per stdv over +/-hrange stdv const hdiv = 3 const hrange = 2 var h [1 + 2*hrange*hdiv]int for _, v := range list { bin := hrange*hdiv + int(math.Floor((v-mean)/stdv*hdiv+.5)) if bin >= 0 && bin < len(h) { h[bin]++ } } const hscale = 10 for _, c := range h { fmt.Println(strings.Repeat("*", (c+hscale/2)/hscale)) }
}</lang>
- Output:
mean 0.995, stdv 0.503 ** **** ****** ******** ************ ************ ************* ************ ********** ******** ***** *** **
Groovy
<lang groovy>rnd = new Random() result = (1..1000).inject([]) { r, i -> r << rnd.nextGaussian() }</lang>
Haskell
<lang haskell>import System.Random
pairs :: [a] -> [(a,a)] pairs (x:y:zs) = (x,y):pairs zs pairs _ = []
gauss mu sigma (r1,r2) =
mu + sigma * sqrt (-2 * log r1) * cos (2 * pi * r2)
gaussians :: (RandomGen g, Random a, Floating a) => Int -> g -> [a] gaussians n g = take n $ map (gauss 1.0 0.5) $ pairs $ randoms g
result :: IO [Double] result = getStdGen >>= \g -> return $ gaussians 1000 g</lang>
Or using Data.Random from random-fu package: <lang haskell>replicateM 1000 $ normal 1 0.5</lang> To print them: <lang haskell>import Data.Random import Control.Monad
thousandRandomNumbers :: RVar [Double] thousandRandomNumbers = replicateM 1000 $ normal 1 0.5
main = do
x <- sample thousandRandomNumbers print x</lang>
HicEst
<lang hicest>REAL :: n=1000, m=1, s=0.5, array(n)
pi = 4 * ATAN(1) array = s * (-2*LOG(RAN(1)))^0.5 * COS(2*pi*RAN(1)) + m </lang>
Icon and Unicon
The seed &random may be assigned in either language; either to randomly seed or to pick a fixed starting point. ?i is the random number generator, returning an integer from 0 to i - 1 for non-zero integer i. As a special case, ?0 yields a random floating point number from 0.0 <= r < 1.0
Note that Unicon randomly seeds it's generator. <lang icon> procedure main()
local L L := list(1000) every L[1 to 1000] := 1.0 + 0.5 * sqrt(-2.0 * log(?0)) * cos(2.0 * &pi * ?0) every write(!L)
end </lang>
IDL
<lang idl>result = 1.0 + 0.5*randomn(seed,1000)</lang>
J
Solution: <lang j>urand=: ?@$ 0: zrand=: (2 o. 2p1 * urand) * [: %: _2 * [: ^. urand
1 + 0.5 * zrand 100</lang>
Alternative Solution:
Using the normal script from the stats/distribs addon.
<lang j> require 'stats/distribs/normal'
1 0.5 rnorm 1000
1.44868803 1.21548637 0.812460657 1.54295452 1.2470606 ...</lang>
Java
<lang java>double[] list = new double[1000]; double mean = 1.0, std = 0.5; Random rng = new Random(); for(int i = 0;i<list.length;i++) {
list[i] = mean + std * rng.nextGaussian();
}</lang>
JavaScript
<lang javascript>function randomNormal() {
return Math.cos(2 * Math.PI * Math.random()) * Math.sqrt(-2 * Math.log(Math.random()))
}
var a = [] for (var i=0; i < 1000; i++){
a[i] = randomNormal() / 2 + 1
}</lang>
jq
Since jq is a purely functional language, it is convenient to define the pseudo-random number generator functions as filters whose inputs and outputs are arrays containing a "seed".
The following uses the same pseudo-random number generator as the Microsoft C Runtime (see Linear congruential generator).
'A Pseudo-Random Number Generator' <lang jq># 15-bit integers generated using the same formula as rand() from the Microsoft C Runtime.
- The random numbers are in [0 -- 32767] inclusive.
- Input: an array of length at least 2 interpreted as [count, state, ...]
- Output: [count+1, newstate, r] where r is the next pseudo-random number.
def next_rand_Microsoft:
.[0] as $count | .[1] as $state | ( (214013 * $state) + 2531011) % 2147483648 # mod 2^31 | [$count+1 , ., (. / 65536 | floor) ] ;</lang>
'Box-Muller Method' <lang jq># Generate a single number following the normal distribution with mean 0, variance 1,
- using the Box-Muller method: X = sqrt(-2 ln U) * cos(2 pi V) where U and V are uniform on [0,1].
- Input: [n, state]
- Output [n+1, nextstate, r]
def next_rand_normal:
def u: next_rand_Microsoft | .[2] /= 32767; u as $u1 | ($u1 | u) as $u2 | ((( (8*(1|atan)) * $u1[2]) | cos) * ((-2 * (($u2[2]) | log)) | sqrt)) as $r | [ (.[0]+1), $u2[1], $r] ;
- Generate "count" arrays, each containing a random normal variate with the given mean and standard deviation.
- Input: [count, state]
- Output: [updatedcount, updatedstate, rnv]
- where "state" is a seed and "updatedstate" can be used as a seed.
def random_normal_variate(mean; sd; count):
next_rand_normal | recurse( if .[0] < count then next_rand_normal else empty end) | .[2] = (.[2] * sd) + mean;</lang>
Example The task can be completed using: [0,1] | random_normal_variate(1; 0.5; 1000) | .[2]
We show just the sample average and standard deviation: <lang jq>def summary:
length as $l | add as $sum | ($sum/$l) as $a | reduce .[] as $x (0; . + ( ($x - $a) | .*. )) | [ $a, (./$l | sqrt)] ;
[ [0,1] | random_normal_variate(1; 0.5; 1000) | .[2] ] | summary</lang>
- Output:
$ jq -n -c -f Random_numbers.jq [0.9932830741018853,0.4977760644490579]
Julia
Julia's standard library provides a randn
function to generate normally distributed random numbers (with mean 0 and standard deviation 0.5, which can be easily rescaled to any desired values):
<lang julia>randn(1000) * 0.5 + 1</lang>
Kotlin
<lang scala>// version 1.0.6
import java.util.Random
fun main(args: Array<String>) {
val r = Random() val da = DoubleArray(1000) for (i in 0 until 1000) da[i] = 1.0 + 0.5 * r.nextGaussian() // now check actual mean and SD val mean = da.average() val sd = Math.sqrt(da.map { (it - mean) * (it - mean) }.average()) println("Mean is $mean") println("S.D. is $sd")
}</lang> Sample output:
- Output:
Mean is 1.0071784073168768 S.D. is 0.48567118114896807
LabVIEW
Liberty BASIC
<lang lb>dim a(1000) mean =1 sd =0.5 for i = 1 to 1000 ' throw 1000 normal variates
a( i) =mean +sd *( sqr( -2 * log( rnd( 0))) * cos( 2 * pi * rnd( 0)))
next i</lang>
Lingo
<lang Lingo>-- Returns a random float value in range 0..1 on randf ()
n = random(the maxinteger)-1 return n / float(the maxinteger-1)
end</lang>
<lang Lingo>normal = [] repeat with i = 1 to 1000
normal.add(1 + sqrt(-2 * log(randf())) * cos(2 * PI * randf()) / 2)
end repeat</lang>
Lobster
Uses built-in rnd_gaussian
<lang Lobster>
let mean = 1.0
let stdv = 0.5
let count = 1000
// stats computes a running mean and variance // See Knuth TAOCP vol 2, 3rd edition, page 232
def stats(xs: [float]) -> float, float: // variance, mean
var M = xs[0] var S = 0.0 var n = 1.0 for(xs.length - 1) i: let x = xs[i + 1] n = n + 1.0 let mm = (x - M) M += mm / n S += mm * (x - M) return (if n > 0.0: S / n else: 0.0), M
def test_random_normal() -> [float]:
rnd_seed(floor(seconds_elapsed() * 1000000)) let r = vector_reserve(typeof return, count) for (count): r.push(rnd_gaussian() * stdv + mean) let cvar, cmean = stats(r) let cstdv = sqrt(cvar) print concat_string(["Mean: ", string(cmean), ", Std.Deviation: ", string(cstdv)], "")
test_random_normal() </lang>
Logo
The earliest Logos only have a RANDOM function for picking a random non-negative integer. Many modern Logos have floating point random generators built-in. <lang logo>to random.float ; 0..1
localmake "max.int lshift -1 -1 output quotient random :max.int :max.int
end
to random.gaussian
output product cos random 360 sqrt -2 / ln random.float
end
make "randoms cascade 1000 [fput random.gaussian / 2 + 1 ?] []</lang>
Lua
<lang lua>local list = {} for i = 1, 1000 do
list[i] = 1 + math.sqrt(-2 * math.log(math.random())) * math.cos(2 * math.pi * math.random()) / 2
end</lang>
M2000 Interpreter
M2000 use a Wichmann - Hill Pseudo Random Number Generator. <lang M2000 Interpreter> Module CheckIt {
Function StdDev (A()) { \\ A() has a copy of values N=Len(A()) if N<1 then Error "Empty Array" M=Each(A()) k=0 \\ make sum, dev same type as A(k) sum=A(k)-A(k) dev=sum \\ find mean While M { sum+=Array(M) } Mean=sum/N \\ make a pointet to A() P=A() \\ subtruct from each item P-=Mean M=Each(P) While M { dev+=Array(M)*Array(M) } \\ as pointer to arrray =(if(dev>0->Sqrt(dev/N), 0), Mean) } Function randomNormal { \\ by default all numbers are double \\ cos() get degrees =1+Cos(360 * rnd) * Sqrt(-2 * Ln(rnd)) /2 } \\ fill array calling randomNormal() for each item Dim A(1000)<<randomNormal() \\ we can pass a pointer to array and place it to stack of values DisplayMeanAndStdDeviation(A()) ' mean ~ 1 std deviation ~0.5 \\ check M2000 rnd only Dim B(1000)<<rnd DisplayMeanAndStdDeviation(B()) ' mean ~ 0.5 std deviation ~0.28
DisplayMeanAndStdDeviation((0,0,14,14)) ' mean = 7 std deviation = 7 DisplayMeanAndStdDeviation((0,6,8,14)) ' mean = 7 std deviation = 5 DisplayMeanAndStdDeviation((6,6,8,8)) ' mean = 7 std deviation = 1 Sub DisplayMeanAndStdDeviation(A) \\ push to stack all items of an array (need an array pointer) Push ! StdDev(A) \\ read from strack two numbers Print "Mean is "; Number Print "Standard Deviation is "; Number End Sub
} Checkit </lang>
Maple
<lang maple>with(Statistics): Sample(Normal(1, 0.5), 1000);</lang>
or
<lang maple>1+0.5*ArrayTools[RandomArray](1000,1,distribution=normal);</lang>
Mathematica
Built-in function RandomReal with built-in distribution NormalDistribution as an argument: <lang Mathematica>RandomReal[NormalDistribution[1, 1/2], 1000]</lang>
MATLAB
Native support : <lang MATLAB> mu = 1; sd = 0.5;
x = randn(1000,1) * sd + mu;
</lang>
The statistics toolbox provides this function <lang MATLAB> x = normrnd(mu, sd, [1000,1]); </lang>
This script uses the Box-Mueller Transform to transform a number from the uniform distribution to a normal distribution of mean = mu0 and standard deviation = chi2.
<lang MATLAB>function randNum = randNorm(mu0,chi2, sz)
radiusSquared = +Inf;
while (radiusSquared >= 1) u = ( 2 * rand(sz) ) - 1; v = ( 2 * rand(sz) ) - 1;
radiusSquared = u.^2 + v.^2; end
scaleFactor = sqrt( ( -2*log(radiusSquared) )./ radiusSquared ); randNum = (v .* scaleFactor .* chi2) + mu0;
end</lang>
Output: <lang MATLAB>>> randNorm(1,.5, [1000,1])
ans =
0.693984121077029</lang>
Maxima
<lang maxima>load(distrib)$
random_normal(1.0, 0.5, 1000);</lang>
MAXScript
<lang maxscript>arr = #() for i in 1 to 1000 do (
a = random 0.0 1.0 b = random 0.0 1.0 c = 1.0 + 0.5 * sqrt (-2*log a) * cos (360*b) -- Maxscript cos takes degrees append arr c
)</lang>
Metafont
Metafont has normaldeviate
which produces pseudorandom normal distributed numbers with mean 0 and variance one. So the following complete the task:
<lang metafont>numeric col[];
m := 0; % m holds the mean, for testing purposes for i = 1 upto 1000:
col[i] := 1 + .5normaldeviate; m := m + col[i];
endfor
% testing m := m / 1000; % finalize the computation of the mean
s := 0; % in s we compute the standard deviation for i = 1 upto 1000:
s := s + (col[i] - m)**2;
endfor s := sqrt(s / 1000);
show m, s; % and let's show that really they get what we wanted end</lang>
A run gave
>> 0.99947 >> 0.50533
Assigning a value to the special variable randomseed will allow to have always the same sequence of pseudorandom numbers
MiniScript
<lang MiniScript>randNormal = function(mean=0, stddev=1)
return mean + sqrt(-2 * log(rnd,2.7182818284)) * cos(2*pi*rnd) * stddev
end function
x = [] for i in range(1,1000)
x.push randNormal(1, 0.5)
end for</lang>
Mirah
<lang mirah>import java.util.Random
list = double[999] mean = 1.0 std = 0.5 rng = Random.new 0.upto(998) do | i |
list[i] = mean + std * rng.nextGaussian
end </lang>
МК-61/52
<lang>П7 <-> П8 1/x П6 ИП6 П9 СЧ П6 1/x ln ИП8 * 2 * КвКор ИП9 2 * пи
- sin * ИП7 + С/П БП 05</lang>
Input: РY - variance, РX - expectation.
Or:
<lang>3 10^x П0 ПП 13 2 / 1 + С/П L0 03 С/П СЧ lg 2 /-/ * КвКор 2 пи ^ СЧ * * cos * В/О</lang>
to generate 1000 numbers with a mean of 1.0 and a standard deviation of 0.5.
Modula-3
<lang modula3>MODULE Rand EXPORTS Main;
IMPORT Random; FROM Math IMPORT log, cos, sqrt, Pi;
VAR rands: ARRAY [1..1000] OF LONGREAL;
(* Normal distribution. *) PROCEDURE RandNorm(): LONGREAL =
BEGIN WITH rand = NEW(Random.Default).init() DO RETURN sqrt(-2.0D0 * log(rand.longreal())) * cos(2.0D0 * Pi * rand.longreal()); END; END RandNorm;
BEGIN
FOR i := FIRST(rands) TO LAST(rands) DO rands[i] := 1.0D0 + 0.5D0 * RandNorm(); END;
END Rand.</lang>
Nanoquery
<lang nanoquery>list = {0} * 1000 mean = 1.0; std = 0.5 rng = new(Nanoquery.Util.Random)
for i in range(0, len(list) - 1)
list[i] = mean + std * rng.getGaussian()
end</lang>
NetRexx
<lang NetRexx>/* NetRexx */ options replace format comments java crossref symbols nobinary
import java.math.BigDecimal import java.math.MathContext
-- prologue numeric digits 20
-- get input, set defaults parse arg dp mu sigma ec . if mu = | mu = '.' then mean = 1.0; else mean = mu if sigma = | sigma = '.' then stdDeviation = 0.5; else stdDeviation = sigma if dp = | dp = '.' then displayPrecision = 1; else displayPrecision = dp if ec = | ec = '.' then elements = 1000; else elements = ec
-- set up RNG = Random() numberList = java.util.List numberList = ArrayList()
-- generate list of random numbers loop for elements
rn = mean + stdDeviation * RNG.nextGaussian() numberList.add(BigDecimal(rn, MathContext.DECIMAL128)) end
-- report say "Mean: " mean say "Standard Deviation:" stdDeviation say "Precision: " displayPrecision say drawBellCurve(numberList, displayPrecision)
return
-- ----------------------------------------------------------------------------- method drawBellCurve(numberList = java.util.List, precision) static
Collections.sort(numberList) val = BigDecimal lastN = nextN = loop val over numberList nextN = Rexx(val.toPlainString()).format(5, precision) select when lastN = then nop when lastN \= nextN then say lastN otherwise nop end say '*\-' lastN = nextN end val say lastN
return
</lang>
- Output:
Mean: 1.0 Standard Deviation: 0.5 Precision: 1 * 2.7 ** 2.5 * 2.4 *** 2.3 ***** 2.2 ******* 2.1 ************* 2.0 ************* 1.9 ***************************** 1.8 ************************* 1.7 ************************************* 1.6 ****************************************************** 1.5 ******************************************** 1.4 ******************************************************************** 1.3 ***************************************************************** 1.2 ************************************************************************** 1.1 ********************************************************************************************* 1.0 ************************************************************* 0.9 ********************************************************************** 0.8 ************************************************************** 0.7 *********************************************************************** 0.6 ************************************************************** 0.5 ****************************************** 0.4 ******************************* 0.3 *************************** 0.2 *************** 0.1 ********* 0.0 ****** -0.1 *** -0.2 *** -0.3 * -0.4 * -0.6 ** -0.7
NewLISP
<lang NewLISP>(normal 1 .5 1000)</lang>
Nim
<lang nim>import math, strutils
const precisn = 5 var rs: TRunningStat
proc normGauss: float {.inline.} = 1 + 0.76 * cos(2*PI*random(1.0)) * sqrt(-2*log10(random(1.0)))
randomize()
for j in 0..5:
for i in 0..1000: rs.push(normGauss()) echo("mean: ", $formatFloat(rs.mean,ffDecimal,precisn), " stdDev: ", $formatFloat(rs.standardDeviation(),ffDecimal,precisn))</lang>
- Output:
mean: 1.01703 stdDev: 0.50324 mean: 1.01187 stdDev: 0.50060 mean: 1.00216 stdDev: 0.49969 mean: 1.00335 stdDev: 0.50184 mean: 1.00120 stdDev: 0.49830 mean: 1.00217 stdDev: 0.49911
Objeck
<lang objeck>bundle Default {
class RandomNumbers { function : Main(args : String[]) ~ Nil { rands := Float->New[1000]; for(i := 0; i < rands->Size(); i += 1;) { rands[i] := 1.0 + 0.5 * RandomNormal(); };
each(i : rands) { rands[i]->PrintLine(); }; }
function : native : RandomNormal() ~ Float { return (2 * Float->Pi() * Float->Random())->Cos() * (-2 * (Float->Random()->Log()))->SquareRoot(); } }
}</lang>
OCaml
<lang ocaml>let pi = 4. *. atan 1.;; let random_gaussian () =
1. +. sqrt (-2. *. log (Random.float 1.)) *. cos (2. *. pi *. Random.float 1.);;
let a = Array.init 1000 (fun _ -> random_gaussian ());;</lang>
Octave
<lang octave>p = normrnd(1.0, 0.5, 1000, 1); disp(mean(p)); disp(sqrt(sum((p - mean(p)).^2)/numel(p)));</lang>
- Output:
1.0209 0.51048
ooRexx
version 1
<lang oorexx>/*REXX pgm gens 1,000 normally distributed #s: mean=1, standard dev.=0.5*/
pi=RxCalcPi() /* get value of pi */ Parse Arg n seed . /* allow specification of N & seed*/ If n==|n==',' Then n=1000 /* N is the size of the array. */ If seed\== Then Call random,,seed /* use seed for repeatable RANDOM#*/ mean=1 /* desired new mean (arith. avg.) */ sd=1/2 /* desired new standard deviation.*/ Do g=1 For n /* generate N uniform random nums.*/ n.g=random(0,1e5)/1e5 /* REXX gens uniform rand integers*/ End
Say ' old mean=' mean() Say 'old standard deviation=' stddev() Say Do j=1 To n-1 By 2 m=j+1 /*use Box-Muller method */ _=sd*RxCalcPower(-2*RxCalcLog(n.j),.5)*RxCalcCos(2*pi*n.m,,'R')+mean n.m=sd*RxCalcpower(-2*RxCalcLog(n.j),.5)*RxCalcSin(2*pi*n.m,,'R')+, mean /* rand # must be 0???1. */ n.j=_ End /* j */ Say ' new mean=' mean() Say 'new standard deviation=' stddev() Exit
mean:
_=0 Do k=1 For n _=_+n.k End Return _/n
stddev:
_avg=mean() _=0 Do k=1 For n _=_+(n.k-_avg)**2 End Return RxCalcPower(_/n,.5)
- requires rxmath library</lang>
- Output:
old mean= 0.49830002 old standard deviation= 0.283199568 new mean= 1.00377404 new standard deviation= 0.501444536
version 2
Using the nice function names in the algorithm. <lang oorexx>/*REXX pgm gens 1,000 normally distributed #s: mean=1, standard dev.=0.5*/
pi=RxCalcPi() /* get value of pi */ Parse Arg n seed . /* allow specification of N & seed*/ If n==|n==',' Then n=1000 /* N is the size of the array. */ If seed\== Then Call random,,seed /* use seed for repeatable RANDOM#*/ mean=1 /* desired new mean (arith. avg.) */ sd=1/2 /* desired new standard deviation.*/ Do g=1 For n /* generate N uniform random nums.*/ n.g=random(0,1e5)/1e5 /* REXX gens uniform rand integers*/ End
Say ' old mean=' mean() Say 'old standard deviation=' stddev() Say Do j=1 To n-1 By 2 m=j+1 /*use Box-Muller method */ _=sd*sqrt(-2*ln(n.j))*cos(2*pi*n.m)+mean n.m=sd*sqrt(-2*ln(n.j))*sin(2*pi*n.m)+mean n.j=_ End Say ' new mean=' mean() Say 'new standard deviation=' stddev() Exit
mean:
_=0 Do k=1 For n _=_+n.k End Return _/n
stddev:
_avg=mean() _=0 Do k=1 For n _=_+(n.k-_avg)**2 End Return sqrt(_/n)
sqrt: Return RxCalcSqrt(arg(1)) ln: Return RxCalcLog(arg(1)) cos: Return RxCalcCos(arg(1),,'R') sin: Return RxCalcSin(arg(1),,'R')
- requires rxmath library</lang>
PARI/GP
<lang parigp>rnormal()={ my(pr=32*ceil(default(realprecision)*log(10)/log(4294967296)),u1=random(2^pr)*1.>>pr,u2=random(2^pr)*1.>>pr); sqrt(-2*log(u1))*cos(2*Pi*u1) \\ Could easily be extended with a second normal at very little cost. }; vector(1000,unused,rnormal()/2+1)</lang>
Pascal
The following function calculates Gaussian-distributed random numbers with the Box-Müller algorithm: <lang pascal> function rnorm (mean, sd: real): real;
{Calculates Gaussian random numbers according to the Box-Müller approach}
var
u1, u2: real;
begin
u1 := random; u2 := random; rnorm := mean * abs(1 + sqrt(-2 * (ln(u1))) * cos(2 * pi * u2) * sd);
end; </lang>
Delphi and Free Pascal support implement a randg function that delivers Gaussian-distributed random numbers.
Perl
<lang perl>my $PI = 2 * atan2 1, 0;
my @nums = map {
1 + 0.5 * sqrt(-2 * log rand) * cos(2 * $PI * rand)
} 1..1000;</lang>
Phix
<lang Phix>function RandomNormal()
return sqrt(-2*log(rnd())) * cos(2*PI*rnd())
end function
sequence s = repeat(0,1000) for i=1 to length(s) do
s[i] = 1 + 0.5 * RandomNormal()
end for</lang>
Phixmonti
<lang Phixmonti>include ..\Utilitys.pmt
def RandomNormal
drop rand log -2 * sqrt 2 pi * rand * cos * 0.5 * 1 +
enddef
1000 var n 0 n repeat
getid RandomNormal map
dup sum n / var mean "Mean: " print mean print nl
0 swap n for
get mean - 2 power rot + swap
endfor swap n / sqrt "Standard deviation: " print print</lang>
PHP
<lang php>function random() {
return mt_rand() / mt_getrandmax();
}
$pi = pi(); // Set PI
$a = array(); for ($i = 0; $i < 1000; $i++) {
$a[$i] = 1.0 + ((sqrt(-2 * log(random())) * cos(2 * $pi * random())) * 0.5);
}</lang>
PicoLisp
<lang PicoLisp>(load "@lib/math.l")
(de randomNormal () # Normal distribution, centered on 0, std dev 1
(*/ (sqrt (* -2.0 (log (rand 0 1.0)))) (cos (*/ 2.0 pi (rand 0 1.0) `(* 1.0 1.0))) 1.0 ) )
(seed (time)) # Randomize
(let Result
(make # Build list (do 1000 # of 1000 elements (link (+ 1.0 (/ (randomNormal) 2))) ) ) (for N (head 7 Result) # Print first 7 results (prin (format N *Scl) " ") ) )</lang>
- Output:
1.500334 1.212931 1.095283 0.433122 0.459116 1.302446 0.402477
PL/I
<lang PL/I> /* CONVERTED FROM WIKI FORTRAN */ Normal_Random: procedure options (main);
declare (array(1000), pi, temp, mean initial (1.0), sd initial (0.5)) float (18); declare (i, n) fixed binary; n = hbound(array, 1); pi = 4.0*ATAN(1.0); array = random(); /* Uniform distribution */ /* Now convert to normal distribution */ DO i = 1 to n-1 by 2; temp = sd * SQRT(-2.0*LOG(array(i))) * COS(2*pi*array(i+1)) + mean; array(i+1) = sd * SQRT(-2.0*LOG(array(i))) * SIN(2*pi*array(i+1)) + mean; array(i) = temp; END; /* Check mean and standard deviation */ mean = SUM(array)/n; sd = SQRT(SUM((array - mean)**2)/n); put skip edit ( "Mean = ", mean ) (a, F(18,16) ); put skip edit ( "Standard Deviation = ", sd) (a, F(18,16));
END Normal_Random; </lang>
- Output:
Mean = 1.0125630677913652 Standard Deviation = 0.5067289784535284 3 runs with different seeds to random(): Mean = 1.0008390411168471 Standard Deviation = 0.5095810511317908 Mean = 0.9754351286894838 Standard Deviation = 0.4804376530558166 Mean = 1.0177411222687990 Standard Deviation = 0.5165899662493400
PL/SQL
<lang PL/SQL> DECLARE
--The desired collection type t_coll is table of number index by binary_integer; l_coll t_coll;
c_max pls_integer := 1000;
BEGIN
FOR l_counter IN 1 .. c_max LOOP -- dbms_random.normal delivers normal distributed random numbers with a mean of 0 and a variance of 1 -- We just adjust the values and get the desired result: l_coll(l_counter) := DBMS_RANDOM.normal * 0.5 + 1; DBMS_OUTPUT.put_line (l_coll(l_counter)); END LOOP;
END; </lang>
Pop11
<lang pop11>;;; Choose radians as arguments to trigonometic functions true -> popradians;
- procedure generating standard normal distribution
define random_normal() -> result; lvars r1 = random0(1.0), r2 = random0(1.0);
cos(2*pi*r1)*sqrt(-2*log(r2)) -> result
enddefine;
lvars array, i;
- Put numbers on the stack
for i from 1 to 1000 do 1.0+0.5*random_normal() endfor;
- collect them into array
consvector(1000) -> array;</lang>
PowerShell
Equation adapted from Liberty BASIC <lang powershell>function Get-RandomNormal
{ [CmdletBinding()] Param ( [double]$Mean, [double]$StandardDeviation ) $RandomNormal = $Mean + $StandardDeviation * [math]::Sqrt( -2 * [math]::Log( ( Get-Random -Minimum 0.0 -Maximum 1.0 ) ) ) * [math]::Cos( 2 * [math]::PI * ( Get-Random -Minimum 0.0 -Maximum 1.0 ) ) return $RandomNormal }
- Standard deviation function for testing
function Get-StandardDeviation
{ [CmdletBinding()] param ( [double[]]$Numbers ) $Measure = $Numbers | Measure-Object -Average $PopulationDeviation = 0 ForEach ($Number in $Numbers) { $PopulationDeviation += [math]::Pow( ( $Number - $Measure.Average ), 2 ) } $StandardDeviation = [math]::Sqrt( $PopulationDeviation / ( $Measure.Count - 1 ) ) return $StandardDeviation }
- Test
$RandomNormalNumbers = 1..1000 | ForEach { Get-RandomNormal -Mean 1 -StandardDeviation 0.5 }
$Measure = $RandomNormalNumbers | Measure-Object -Average
$Stats = [PSCustomObject]@{
Count = $Measure.Count Average = $Measure.Average StandardDeviation = Get-StandardDeviation -Numbers $RandomNormalNumbers
}
$Stats | Format-List </lang>
- Output:
Count : 1000 Average : 1.01206560135809 StandardDeviation : 0.489099623426272
PureBasic
<lang PureBasic>Procedure.f RandomNormal()
; This procedure can return any real number. Protected.f x1, x2
; random numbers from the open interval ]0, 1[ x1 = (Random(999998)+1) / 1000000 ; must be > 0 because of Log(x1) x2 = (Random(999998)+1) / 1000000
ProcedureReturn Sqr(-2*Log(x1)) * Cos(2*#PI*x2)
EndProcedure
Define i, n=1000
Dim a.q(n-1) For i = 0 To n-1
a(i) = 1 + 0.5 * RandomNormal()
Next</lang>
Python
- Using random.gauss
<lang python>>>> import random >>> values = [random.gauss(1, .5) for i in range(1000)] >>> </lang>
- Quick check of distribution
<lang python>>>> def quick_check(numbers):
count = len(numbers) mean = sum(numbers) / count sdeviation = (sum((i - mean)**2 for i in numbers) / count)**0.5 return mean, sdeviation
>>> quick_check(values) (1.0140373306786599, 0.49943411329234066) >>> </lang>
Note that the random module in the Python standard library supports a number of statistical distribution methods.
- Alternatively using random.normalvariate
<lang python>>>> values = [ random.normalvariate(1, 0.5) for i in range(1000)] >>> quick_check(values) (0.990099111944864, 0.5029847005836282) >>> </lang>
R
<lang r>result <- rnorm(1000, mean=1, sd=0.5)</lang>
Racket
<lang Racket>
- lang racket
(for/list ([i 1000])
(add1 (* (sqrt (* -2 (log (random)))) (cos (* 2 pi (random))) 0.5)))
</lang>
Alternative: <lang Racket>
- lang racket
(require math/distributions) (sample (normal-dist 1.0 0.5) 1000) </lang>
Raku
(formerly Perl 6)
<lang perl6>sub randnorm ($mean, $stddev) {
$mean + $stddev * sqrt(-2 * log rand) * cos(2 * pi * rand)
}
my @nums = randnorm(1, 0.5) xx 1000;
- Checking
say my $mean = @nums R/ [+] @nums; say my $stddev = sqrt $mean**2 R- @nums R/ [+] @nums X** 2; </lang>
Raven
<lang raven>define PI
-1 acos
define rand1
9999999 choose 1 + 10000000.0 /
define randNormal
rand1 PI * 2 * cos rand1 log -2 * sqrt * 2 / 1 +
1000 each drop randNormal "%f\n" print</lang> Quick Check (on linux with code in file rand.rv) <lang raven>raven rand.rv | awk '{sum+=$1; sumsq+=$1*$1;} END {print "stdev = " sqrt(sumsq/NR - (sum/NR)**2); print "mean = " sum/NR}' stdev = 0.497773 mean = 1.01497</lang>
REXX
The REXX language doesn't have any "higher math" functions like SQRT/SIN/COS/LN/LOG/EXP/POW/etc.,
so we hoi polloi REXX programmers have to roll our own.
Programming note: note the range of the random numbers: (0,1]
(that is, random numbers from zero──►unity, excluding zero, including unity).
<lang rexx>/*REXX pgm generates 1,000 normally distributed numbers: mean=1, standard deviation=½.*/
numeric digits 20 /*the default decimal digit precision=9*/
parse arg n seed . /*allow specification of N and the seed*/
if n== | n=="," then n=1000 /*N: is the size of the array. */
if datatype(seed,'W') then call random ,,seed /*SEED: for repeatable random numbers. */
newMean=1 /*the desired new mean (arithmetic avg)*/
sd=1/2 /*the desired new standard deviation. */
do g=1 for n /*generate N uniform random #'s (0,1].*/ #.g = random(1, 1e5) / 1e5 /*REXX's RANDOM BIF generates integers.*/ end /*g*/ /* [↑] random integers ──► fractions. */
say ' old mean=' mean() say 'old standard deviation=' stdDev() call pi; pi2=pi * 2 /*define pi and also 2 * pi. */ say
do j=1 to n-1 by 2; m=j+1 /*step through the iterations by two. */ _=sd * sqrt(ln(#.j) * -2) /*calculate the used-twice expression.*/ #.j=_ * cos(pi2 * #.m) + newMean /*utilize the Box─Muller method. */ #.m=_ * sin(pi2 * #.m) + newMean /*random number must be: (0,1] */ end /*j*/
say ' new mean=' mean() say 'new standard deviation=' stdDev() exit /*stick a fork in it, we're all done. */ /*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ mean: _=0; do k=1 for n; _=_ + #.k; end; return _/n stdDev: _avg=mean(); _=0; do k=1 for n; _=_ + (#.k - _avg)**2; end; return sqrt(_/n) e: e =2.7182818284590452353602874713526624977572470936999595749669676277240766303535; return e /*digs overkill*/ pi: pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862; return pi /* " " */ r2r: return arg(1) // (pi() * 2) /*normalize ang*/ sin: procedure; parse arg x;x=r2r(x);numeric fuzz min(5,digits()-3);if abs(x)=pi then return 0;return .sincos(x,x,1) .sincos:parse arg z,_,i; x=x*x; p=z; do k=2 by 2; _=-_*x/(k*(k+i)); z=z+_; if z=p then leave; p=z; end; return z /*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ ln: procedure; parse arg x,f; call e; ig= x>1.5; is=1 - 2 * (ig\==1); ii=0; xx=x
do while ig&xx>1.5|\ig&xx<.5;_=e;do k=-1;iz=xx*_**-is;if k>=0&(ig&iz<1|\ig&iz>.5) then leave;_=_*_;izz=iz;end xx=izz;ii=ii+is*2**k;end;x=x*e**-ii-1;z=0;_=-1;p=z;do k=1;_=-_*x;z=z+_/k;if z=p then leave;p=z;end; return z+ii
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ cos: procedure; parse arg x; x=r2r(x); a=abs(x); hpi=pi * .5
numeric fuzz min(6, digits() - 3); if a=pi then return -1 if a=hpi | a=hpi*3 then return 0; if a=pi/3 then return .5 if a=pi * 2/3 then return -.5; return .sinCos(1,1,-1)
/*───────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); numeric digits; h=d+6
numeric form; parse value format(x,2,1,,0) 'E0' with g 'E' _ .; g=g * .5'e'_ %2 m.=9; do j=0 while h>9; m.j=h; h=h%2 + 1; end /*j*/ do k=j+5 to 0 by -1; numeric digits m.k; g=(g+x/g)*.5; end /*k*/ numeric digits d; return g/1</lang>
output when using the default inputs:
old mean= 0.5015724 old standard deviation= 0.28652466389342471402 new mean= 0.98807025356443262689 new standard deviation= 0.50002924192766720838
Ring
<lang ring> for i = 1 to 10
see random(i) + nl
next i </lang>
Ruby
<lang ruby>Array.new(1000) { 1 + Math.sqrt(-2 * Math.log(rand)) * Math.cos(2 * Math::PI * rand) }</lang>
Run BASIC
<lang runbasic>dim a(1000) pi = 22/7 for i = 1 to 1000
a( i) = 1 + .5 * (sqr(-2 * log(rnd(0))) * cos(2 * pi * rnd(0)))
next i</lang>
Rust
Using a for-loop: <lang rust>extern crate rand; use rand::distributions::{Normal, IndependentSample};
fn main() {
let mut rands = [0.0; 1000]; let normal = Normal::new(1.0, 0.5); let mut rng = rand::thread_rng(); for num in rands.iter_mut() { *num = normal.ind_sample(&mut rng); }
}</lang>
Using iterators: <lang rust>extern crate rand; use rand::distributions::{Normal, IndependentSample};
fn main() {
let rands: Vec<_> = { let normal = Normal::new(1.0, 0.5); let mut rng = rand::thread_rng(); (0..1000).map(|_| normal.ind_sample(&mut rng)).collect() };
}</lang>
SAS
<lang SAS> /* Generate 1000 random numbers with mean 1 and standard deviation 0.5.
SAS version 9.2 was used to create this code.*/
data norm1000;
call streaminit(123456);
/* Set the starting point, so we can replicate results.
If you want different results each time, comment the above line. */ do i=1 to 1000; r=rand('normal',1,0.5); output; end;
run; </lang> Results:
The MEANS Procedure Analysis Variable : r Mean Std Dev ---------------------------- 0.9907408 0.4844051 ----------------------------
Sather
<lang sather>class MAIN is
main is a:ARRAY{FLTD} := #(1000); i:INT;
RND::seed(2010); loop i := 1.upto!(1000) - 1; a[i] := 1.0d + 0.5d * RND::standard_normal; end;
-- testing the distribution mean ::= a.reduce(bind(_.plus(_))) / a.size.fltd; #OUT + "mean " + mean + "\n"; a.map(bind(_.minus(mean))); a.map(bind(_.pow(2.0d))); dev ::= (a.reduce(bind(_.plus(_))) / a.size.fltd).sqrt; #OUT + "dev " + dev + "\n"; end;
end;</lang>
Scala
One liner
<lang scala>List.fill(1000)(1.0 + 0.5 * scala.util.Random.nextGaussian)</lang>
Academic
<lang scala> object RandomNumbers extends App {
val distribution: LazyList[Double] = { def randomNormal: Double = 1.0 + 0.5 * scala.util.Random.nextGaussian
def normalDistribution(a: Double): LazyList[Double] = a #:: normalDistribution(randomNormal)
normalDistribution(randomNormal) }
/* * Let's test it */ def calcAvgAndStddev[T](ts: Iterable[T])(implicit num: Fractional[T]): (T, Double) = { val mean: T = num.div(ts.sum, num.fromInt(ts.size)) // Leaving with type of function T
// Root of mean diffs val stdDev = Math.sqrt(ts.map { x => val diff = num.toDouble(num.minus(x, mean)) diff * diff }.sum / ts.size)
(mean, stdDev) }
println(calcAvgAndStddev(distribution.take(1000))) // e.g. (1.0061433267806525,0.5291834867560893)
} </lang>
Scheme
<lang scheme>; linear congruential generator given in C99 section 7.20.2.1 (define ((c-rand seed)) (set! seed (remainder (+ (* 1103515245 seed) 12345) 2147483648)) (quotient seed 65536))
- uniform real numbers in open interval (0, 1)
(define (unif-rand seed) (let ((r (c-rand seed))) (lambda () (/ (+ (r) 1) 32769.0))))
- Box-Muller method to generate normal distribution
(define (normal-rand unif m s) (let ((? #t) (! 0.0) (twopi (* 2.0 (acos -1.0)))) (lambda ()
(set! ? (not ?)) (if ? ! (let ((a (sqrt (* -2.0 (log (unif))))) (b (* twopi (unif)))) (set! ! (+ m (* s a (sin b)))) (+ m (* s a (cos b))))))))
(define rnorm (normal-rand (unif-rand 0) 1.0 0.5))
- auxiliary function to get a list of 'n random numbers from generator 'r
(define (rand-list r n) = (if (zero? n) '() (cons (r) (rand-list r (- n 1)))))
(define v (rand-list rnorm 1000))
v
- |
(-0.27965824722565835
-0.8870860825789542 0.6499618744638194 0.31336141955110863 ... 0.5648743998193049 0.8282656735558756 0.6399951934564637 0.7699535302478072)
|#
- check mean and standard deviation
(define (mean-sdev v) (let loop ((v v) (a 0) (b 0) (n 0)) (if (null? v)
(let ((mean (/ a n))) (list mean (sqrt (/ (- b (* n mean mean)) (- n 1))))) (let ((x (car v))) (loop (cdr v) (+ a x) (+ b (* x x)) (+ n 1))))))
(mean-sdev v)
- (0.9562156817697293 0.5097087109575911)</lang>
Seed7
<lang seed7>$ include "seed7_05.s7i";
include "float.s7i"; include "math.s7i";
const func float: frand is func # Uniform distribution, (0..1]
result var float: frand is 0.0; begin repeat frand := rand(0.0, 1.0); until frand <> 0.0; end func;
const func float: randomNormal is # Normal distribution, centered on 0, std dev 1
return sqrt(-2.0 * log(frand)) * cos(2.0 * PI * frand);
const proc: main is func
local var integer: i is 0; var array float: rands is 1000 times 0.0; begin for i range 1 to length(rands) do rands[i] := 1.0 + 0.5 * randomNormal; end for; end func;</lang>
Sidef
<lang ruby>var arr = 1000.of { 1 + (0.5 * sqrt(-2 * 1.rand.log) * cos(Num.tau * 1.rand)) } arr.each { .say }</lang>
Standard ML
SML/NJ has two structures for random numbers:
1) Rand (a linear congruential generator).
You create the generator by calling Rand.mkRandom
with a seed (of word
type).
You can call the generator with ()
repeatedly to get a word in the range [Rand.randMin, Rand.randMax]
.
You can use the Rand.norm
function to transform the output into a real
from 0 to 1, or use the Rand.range (i,j)
function to transform the output into an int
of the given range.
<lang sml>val seed = 0w42;
val gen = Rand.mkRandom seed;
fun random_gaussian () =
1.0 + Math.sqrt (~2.0 * Math.ln (Rand.norm (gen ()))) * Math.cos (2.0 * Math.pi * Rand.norm (gen ()));
val a = List.tabulate (1000, fn _ => random_gaussian ());</lang>
2) Random (a subtract-with-borrow generator). You create the generator by calling Random.rand
with a seed (of a pair of int
s). You can use the Random.randInt
function to generate a random int over its whole range; Random.randNat
to generate a non-negative random int; Random.randReal
to generate a real
between 0 and 1; or Random.randRange (i,j)
to generate an int
in the given range.
<lang sml>val seed = (47,42);
val gen = Random.rand seed;
fun random_gaussian () =
1.0 + Math.sqrt (~2.0 * Math.ln (Random.randReal gen)) * Math.cos (2.0 * Math.pi * Random.randReal gen);
val a = List.tabulate (1000, fn _ => random_gaussian ());</lang>
Other implementations of Standard ML have their own random number generators. For example, Moscow ML has a Random
structure that is different from the one from SML/NJ.
The SML Basis Library does not provide a routine for uniform deviate generation, and PolyML does not have one. Using a routine from "Monte Carlo" by Fishman (Springer), in the function uniformdeviate, and avoiding the slow IntInf's: <lang smlh> val urandomlist = fn seed => fn n => let val uniformdeviate = fn seed => let val in31m = (Real.fromInt o Int32.toInt ) (getOpt (Int32.maxInt,0) ); val in31 = in31m +1.0; val s1 = 41160.0; val s2 = 950665216.0; val v = Real.realFloor seed; val val1 = v*s1; val val2 = v*s2; val next1 = Real.fromLargeInt (Real.toLargeInt IEEEReal.TO_NEGINF (val1/in31)) ; val next2 = Real.rem(Real.realFloor(val2/in31) , in31m ); val valt = val1+val2 - (next1+next2)*in31m; val nextt = Real.realFloor(valt/in31m); val valt = valt - nextt*in31m; in (valt/in31m,valt) end; val store = ref (0.0,0.0); val rec u = fn S => fn 0 => [] | n=> (store:=uniformdeviate S; (#1 (!store)):: (u (#2 (!store)) (n-1))) ; in u seed n end;
local open Math in val bmconv = fn urand => fn vrand => 1.0+0.5*(sqrt(~2.0*ln urand)*cos (2.0*pi*vrand) ) end;
val rec makeNormals = fn once => fn u::v::[] => [once u v] | u::v::rm => (once u v )::(makeNormals once rm );
val anyrealseed=1009.0 ; makeNormals bmconv (urandomlist anyrealseed 2000); </lang>
Stata
<lang stata>clear all set obs 1000 gen x=rnormal(1,0.5)</lang>
Mata
<lang stata>a = rnormal(1000,1,1,0.5)</lang>
Tcl
<lang tcl>package require Tcl 8.5 variable ::pi [expr acos(0)] proc ::tcl::mathfunc::nrand {} {
expr {sqrt(-2*log(rand())) * cos(2*$::pi*rand())}
}
set mean 1.0 set stddev 0.5 for {set i 0} {$i < 1000} {incr i} {
lappend result [expr {$mean + $stddev*nrand()}]
}</lang>
TI-83 BASIC
Builtin function: randNorm()
randNorm(1,.5)
Or by a program:
Calculator symbol translations:
"STO" arrow: →
Square root sign: √
ClrList L1 Radian For(A,1,1000) √(-2*ln(rand))*cos(2*π*A)→L1(A) End
TorqueScript
<lang tqs>for (%i = 0; %i < 1000; %i++) %list[%i] = 1 + mSqrt(-2 * mLog(getRandom())) * mCos(2 * $pi * getRandom());</lang>
Ursala
There are two ways of interpreting the task, either to simulate sampling a population described by the given statistics, or to construct a sample exhibiting the given statistics. Both are illustrated below. The functions parameterized by the mean and standard deviation take a sample size and return a sample of that size, represented as a list of floating point numbers. The Z library function simulates a draw from a standard normal distribution. Mean and standard deviation library functions are also used in this example. <lang Ursala>#import nat
- import flo
pop_stats("mu","sigma") = plus/*"mu"+ times/*"sigma"+ Z*+ iota
sample_stats("mu","sigma") = plus^*D(minus/"mu"+ mean,~&)+ vid^*D(div\"sigma"+ stdev,~&)+ Z*+ iota
- cast %eWL
test =
^(mean,stdev)* <
pop_stats(1.,0.5) 1000, sample_stats(1.,0.5) 1000></lang>
The output shows the mean and standard deviation for both sample vectors, the latter being exact by construction.
< (1.004504e+00,4.915525e-01), (1.000000e+00,5.000000e-01)>
Visual FoxPro
<lang vfp> LOCAL i As Integer, m As Double, n As Integer, sd As Double py = PI() SET TALK OFF SET DECIMALS TO 6 CREATE CURSOR gdev (deviate B(6)) RAND(-1) n = 1000 m = 1 sd = 0.5 CLEAR FOR i = 1 TO n INSERT INTO gdev VALUES (GaussDev(m, 1/sd)) ENDFOR CALCULATE AVG(deviate), STD(deviate) TO m, sd ? "Mean", m, "Std Dev", sd SET TALK ON SET DECIMALS TO
FUNCTION GaussDev(mean As Double, sdev As Double) As Double LOCAL z As Double z = SQRT(-2*LOG(RAND()))*COS(py*RAND()) IF sdev # 0 z = mean + z/sdev ENDIF RETURN z ENDFUNC </lang>
Wren
<lang ecmascript>import "random" for Random
var rand = Random.new()
var randNormal = Fn.new { (-2 * rand.float().log).sqrt * (2 * Num.pi * rand.float()).cos }
var stdDev = Fn.new { |a, m|
var c = a.count return ((a.reduce(0) { |acc, x| acc + x*x } - m*m*c) / c).sqrt
}
var n = 1000 var numbers = List.filled(n, 0) var mu = 1 var sigma = 0.5 var sum = 0 for (i in 0...n) {
numbers[i] = mu + sigma*randNormal.call() sum = sum + numbers[i]
} var mean = sum / n System.print("Actual mean : %(mean)") System.print("Actual std dev: %(stdDev.call(numbers, mean))")</lang>
- Output:
Sample run:
Actual mean : 1.0053988699746 Actual std dev: 0.4961645117026
Yorick
Returns array of count random numbers with mean 0 and standard deviation 1. <lang yorick>func random_normal(count) {
return sqrt(-2*log(random(count))) * cos(2*pi*random(count));
}</lang>
Example of basic use:
> nums = random_normal(1000); // create an array 1000 random numbers > nums(avg); // show the mean 0.00901216 > nums(rms); // show the standard deviation 0.990265
Example with a mean of 1.0 and a standard deviation of 0.5:
> nums = random_normal(1000) * 0.5 + 1; > nums(avg); 1.00612 > nums(rms); 0.496853
zkl
<lang zkl>fcn mkRand(mean,sd){ //normally distributed random w/mean & standard deviation
pi:=(0.0).pi; // using the Box–Muller transform rz1:=fcn{1.0-(0.0).random(1)} // from [0,1) to (0,1] return('wrap(){((-2.0*rz1().log()).sqrt() * (2.0*pi*rz1()).cos())*sd + mean })
}</lang> This creates a new random number generator, now to use it: <lang zkl>var g=mkRand(1,0.5); ns:=(0).pump(1000,List,g); // 1000 rands with mean==1 & sd==1/2 mean:=(ns.sum(0.0)/1000); //-->1.00379
// calc sd of list of numbers:
(ns.reduce('wrap(p,n){p+(n-mean).pow(2)},0.0)/1000).sqrt() //-->0.494844</lang>
ZX Spectrum Basic
Here we have converted the QBasic code to suit the ZX Spectrum:
<lang zxbasic>10 RANDOMIZE 0 : REM seeds random number generator based on uptime 20 DIM a(1000) 30 CLS 40 FOR i = 1 TO 1000 50 LET a(i) = 1 + SQR(-2 * LN(RND)) * COS(2 * PI * RND) 60 NEXT i</lang>
- Programming Tasks
- Basic language learning
- Probability and statistics
- Randomness
- GUISS/Omit
- UNIX Shell/Omit
- Ada
- ALGOL 68
- Arturo
- AutoHotkey
- Avail
- AWK
- BASIC
- BBC BASIC
- C
- C sharp
- C++
- Boost
- Clojure
- Common Lisp
- Crystal
- D
- Tango
- Delphi
- DWScript
- E
- EasyLang
- Eiffel
- Elena
- Elixir
- Erlang
- ERRE
- Euler Math Toolbox
- Euphoria
- F Sharp
- Factor
- Falcon
- Fantom
- Forth
- Fortran
- Free Pascal
- FreeBASIC
- Frink
- FutureBasic
- Go
- Groovy
- Haskell
- HicEst
- Icon
- Unicon
- IDL
- J
- Java
- JavaScript
- Jq
- Julia
- Kotlin
- LabVIEW
- Liberty BASIC
- Lingo
- Lobster
- Logo
- Lua
- M2000 Interpreter
- Maple
- Mathematica
- MATLAB
- Maxima
- MAXScript
- Metafont
- MiniScript
- Mirah
- МК-61/52
- Modula-3
- Nanoquery
- NetRexx
- NewLISP
- Nim
- Objeck
- OCaml
- Octave
- OoRexx
- PARI/GP
- Pascal
- Perl
- Phix
- Phixmonti
- PHP
- PicoLisp
- PL/I
- PL/SQL
- Pop11
- PowerShell
- PureBasic
- Python
- R
- Racket
- Raku
- Raven
- REXX
- Ring
- Ruby
- Run BASIC
- Rust
- Rand
- SAS
- Sather
- Scala
- Scheme
- Seed7
- Sidef
- Standard ML
- Stata
- Tcl
- TI-83 BASIC
- TorqueScript
- Ursala
- Visual FoxPro
- Wren
- Yorick
- Zkl
- ZX Spectrum Basic