# Sample the Pareto Distribution

The Pareto Distribution is a probability distribution that comes up a lot in nature. It has lots of special properties, such as an infinite mean. In this challenge, you will output a number sampled from this distribution.

The Pareto Distribution is defined to be greater than or equal to x with probability 1/x, for all x greater than or equal to 1.

Therefore, a number sampled from this distribution is greater than or equal to 1 with probability 1, greater than or equal to 2 with probability exactly 1/2, greater than or equal to 3 with probability exactly 1/3, greater than or equal to 11.4 with probability exactly 1/11.4, and so on.

Since you will sample this distribution, your program or function will take no input, and output a random number, with the above probabilities. However, if your program doesn't perfectly match the above probabilities due to floating-point impression, that's OK. See the bottom of the challenge for more details.

(This is called the Pareto Distribution with alpha 1 and lower bound 1, to be exact)

Here's 10 example draws from this distribution:

1.1540029602790338
52.86156818209856
3.003306506971116
1.4875532217142287
1.3604286212876546
57.5263129600285
1.3139866916055676
20.25125817471419
2.8105749663695208
1.1528212409680156


Notice how 5 of them are below 2, and 5 are above 2. Since this is the average result, it could have been higher or lower, of course.

Your answer only needs to be correct up to the limits of your floating point type, real number type, or whatever else you use, but you must be able to represent numbers at at least 3 decimal digits of precision, and represent numbers up to 1,000,000. If you're not sure whether something is OK, feel free to ask.

This is code golf.

• For each range [a, b], where 1 <= a < b, the ideal probability that the sample would fall in that range is 1/a - 1/b. The probability that your program produces a number in that range must be with 0.001 of 1/a - 1/b. If X is the output of your program, it is required that |P(a <= X <= b) - (1/a - 1/b)| < 0.001.

• Note that by applying the above rule with a=1 and b sufficiently large, it is the case that your program must output a number greater than or equal to 1 with at least probability 0.999. The rest of the time it may crash, output Infinity, or do whatever else.

I'm fairly certain that the existing submissions of the form 1/1-x or 1/x, where x is a random float in [0, 1) or (0, 1) or [0, 1], all satisfy this requirement.

• Dec 16, 2017 at 5:46
• Note to everyone: issacg has added some rules that allow some imprecisions, therefore most answers here are longer than necessary. [sorry for comment abuse too, but that is what would happen when OP change question significantly] Dec 16, 2017 at 15:16

# MATL, 3 bytes

1r/


Try it online! Or estimate the resulting probabilities by running it 10000 times.

### Explanation

1    % Push 1
r    % Push random number uniformly distributed on the open interval (0,1)
/    % Divide. Implicitly display

• After the clarification and edit in the OP, this answer conforms to the rules of the challenge Dec 17, 2017 at 4:11

# Actually, 4 bytes

G1-ì


Try it online!

Explanation:

G1-ì
G     random()
1-   1-random()
ì  1/(1-random())


# R, 10 bytes

1/runif(1)


Pretty straightforward.

• Note that runif never returns 0 or 1 in the default case so there are no problems with this. Dec 15, 2017 at 12:31
• Yes, thanks. And I didn't thought of it when entering this answer but you can indeed verify the distribution if needed. Dec 15, 2017 at 12:50
• @Mego that is incorrect. The Pareto distribution is absolutely continuous and thus has measure 0 for any number. Dec 15, 2017 at 13:38
• @Mego OK that may be quicksand for me (given i know close to nothing about floating point), but i actually think that while the probability of runif giving 1 is null, the probability of 1/runif to give 1 is not, because of floating point accuracy (i. e. typically 1/0.9999999 returns 1 in R). Dec 15, 2017 at 14:59
• @plannapus Hmm... That's a good point. Floats make this entirely too complicated.
– user45941
Dec 15, 2017 at 15:00

# TI-Basic, 2 bytes

rand^-1      (AB 0C in hex)


For anyone who is wondering, rand returns a random value in (0,1]. "Due to specifics of the random number generating algorithm, the smallest number possible to generate is slightly greater than 0. The largest number possible is actually 1 ..." (source). For example, seeding rand with 196164532 yields 1.

• Strangely, the equivalent code wouldn't work on a TI-89 series calculator. Even though their random number generators are nearly identically implemented, a TI-89 will return 0 whenever a TI-83+ would return 0.99999999999889. Dec 15, 2017 at 17:17
• TI-Basic developers knew in advance this challenge will happen...? It seems to win this time. Dec 16, 2017 at 6:13
• @user202729 Avoiding 0 and 1 makes rand more useful as a subroutine for the calculator's other commands, which is probably why TI made this design decision. For example, randNorm(0,1 returns -7.02129... with seed 196164532. Using the RNG algorithm without the adjustment would give a value of 1e99, which is an unreasonable value for a normally-distributed variable to have. Dec 16, 2017 at 14:55
• @user202729 Yeah, actually I just time traveled a bit to get it all done. Definitely worth it for these upvotes. Dec 16, 2017 at 15:01

# R, 12 bytes

exp(rexp(1))


Try it online!

Verify the distribution

This takes a different approach, exploiting the fact that if Y~exp(alpha), then X=x_m*e^Y is a Pareto with parameters x_m,alpha. Since both parameters are 1, and the default rate parameter for rexp is 1, this results in the appropriate Pareto distribution.

While this answer is a fairly R- specific approach, it's sadly less golfy than plannapus'.

# R, 14 bytes

1/rbeta(1,1,1)


Try it online!

Even less golfy, but another way of getting at the answer.

Another property of the exponential distribution is that if X ~ Exp(λ) then e^−X ~ Beta(λ, 1), hence 1/Beta(1,1) is a Pareto(1,1).

Additionally, a keen observer would recall that if X ~ Beta(a,b) and a=b=1, then X~Unif(0,1), so this truly is 1/runif(1).

• I have no idea. But the reality is, there is a huge confusion about what is allowed and what is not on this challenge. Dec 15, 2017 at 15:22
• @user202729 that's fair, but those who have been raising concerns with regards to that would at least have commented, so the downvote is (in my opinion) unlikely to be related to that. EDIT: mystery downvoter has removed the downvote. Dec 15, 2017 at 15:25
• I downvoted because I thought using R on a challenge like this was trivial, but I was a little trigger-happy. I realize that this uses a different method than most of the other answers, so I removed my downvote. Dec 15, 2017 at 15:27
• @KSmarts The "trivial" answer in R wasn't used by anybody actually: actuar::rpareto(1,1,1), because it is longer :) Dec 15, 2017 at 15:31
• For info, there are ca. 20 distributions hard-coded in base R, but Pareto is not one of them, hence the need to either use a work-around or an additional package. Dec 15, 2017 at 15:34

# Charcoal, 10 bytes

Ｉ∕Ｘφ²⊕‽Ｘφ²


Try it online!

Link is to the verbose version:

Print(Cast(Divide(Power(f, 2), ++(Random(Power(f, 2))))));


• Charcoal only has methods to get random integer numbers, so in order to get a random floating-point number between 0 and 1 we have to get a random integer between 0 and N and divide by N.
• Previous version of this answer that used the 1/(1-R) formula: In this case, N is set to 1000000 as the OP asks it to be the minimum. To get this number Charcoal provides a preset variable f=1000. So just calculating f^2 we get 1000000. In the event that the random number is 999999 (the maximum), 1/(1-0.999999)=1000000.
• Neil's tip (saving 3 bytes): If I have 1/(1-R/N) where R is a random number between 0 and N, it is the same as just calculate N/(N-R). But considering that the random integers N-R and R have the same probability to occur, that is the same as just calculating N/R (being R in this last case a number between 1 and N inclusive to avoid division by zero).
• 10 bytes
– Neil
Dec 15, 2017 at 9:23
• @Neil please wait a moment while I try to understand what your code does... :-) Dec 15, 2017 at 9:25
• Actually I don't need MapAssignRight any more, 10 bytes! works.
– Neil
Dec 15, 2017 at 9:27
• @Neil assimilation of your code completed! Answer edited. :-D Dec 15, 2017 at 9:43

The function randomIO :: IO Float produces random numbers in the interval [0,1), so transforming them using x -> 1/(1-x) will produce pareto realizations.

import System.Random
randomIO>>=print.(1/).((1::Float)-)


Try it online!

• Moving the type annotation saves a few bytes: randomIO>>=print.((1::Float)/) Dec 15, 2017 at 11:55
• And as functions are allowed, I'd say you can drop the main=. Dec 15, 2017 at 11:57
• Appranetly the range is [0,1) according to this answer Dec 15, 2017 at 13:30
• @flawr Whoops, you're right! I forgot how floats work temporarily.
– user45941
Dec 15, 2017 at 13:32
• Well anyway, thanks for commenting, I would not have had any idea:) Dec 15, 2017 at 13:35

# Excel, 9 bytes

=1/rand()


Yay, Excel is (semi-)competitive for a change!

• Also works in LibreOffice Calc :) Dec 16, 2017 at 7:59
• You can change this to google sheets for -1 Bytes (=1/Rand() Dec 18, 2017 at 15:04

# Mathematica, 10 bytes

1/Random[]


Try it online!

-4 bytes from M.Stern

• This has the potential to fail, because RandomReal outputs a real number in the closed range [0, 1]. Thus, division by 0 is possible. You’ll need to manipulate the random value in order to remove that possibility.
– user45941
Dec 15, 2017 at 11:42
• @Mego where exactly did you find that info? Dec 15, 2017 at 11:45
• @Mego what is the probability to get 0? Dec 15, 2017 at 12:20
• Jenny_mathy: According to the proposal on meta, the burden of proof should be on the person claiming to have a valid answer - it's your work to prove that it's valid, not to ask @Mego to provide an invalid test case. Also because float are discrete the probability to get 0 is nonzero. Dec 15, 2017 at 12:52
• Back to the topic, I do not believe there is a possibility of getting a zero using this function. Mathematica will in fact produce numbers less than $MinMachineNumber. Try this: Table[RandomReal[{0,$MinMachineNumber}], 100]. Turns out Mathematica is smart enough to abandon machine numbers and switch over to arbitrary precision numbers. LOL. Dec 15, 2017 at 15:11

# Ruby, 14 8 bytes

p 1/rand


Trivial program, I don't think it can get any shorter.

• Note to everyone: issacg has added some rules that allow some imprecisions, therefore most answers here are longer than necessary. Dec 16, 2017 at 15:20

# Excel VBA, 6 Bytes

Anonymous VBE immediate window function that takes no input and outputs to the VBE immediate window

?1/Rnd


# Python, 41 bytes

lambda:1/(1-random())
from random import*


Try it online!

Using the builtin is actually longer:

# Python, 43 bytes

lambda:paretovariate(1)
from random import*


Try it online!

Both solutions work in both Python 2 and Python 3.

• Full programs are shorter for tasks which don't use input, using print saves a byte. Dec 15, 2017 at 13:57

# J, 5 bytes

%-.?0


## How ot works:

?0 generates a random value greater than 0 and less than 1

-. subtract from 1

% reciprocal

Try it online!

• Note to everyone: issacg has added some rules that allow some imprecisions, therefore most answers here are longer than necessary. Dec 16, 2017 at 15:18

# Red, 19 bytes

1 /(1 - random 1.0)


Try it online!

• Note to everyone: issacg has added some rules that allow some imprecisions, therefore most answers here are longer than necessary. Dec 16, 2017 at 15:18

# APL (Dyalog), 5 bytes

÷1-?0


Try it online!

How?

 ÷   1-     ?0
1÷  (1-  random 0..1)

• Note to everyone: issacg has added some rules that allow some imprecisions, therefore most answers here are longer than necessary. Dec 16, 2017 at 15:18

# Japt, 6 bytes

1/1-Mr is the same length but this felt a little less boring!

°T/aMr


Try it

## Explanation

Increment (°) zero (T) and divide by (/) its absolute difference (a) with Math.random().

• Note to everyone: issacg has added some rules that allow some imprecisions, therefore most answers here are longer than necessary. Dec 16, 2017 at 15:17

# Jelly, 5 bytes

Jelly also doesn't have random float, so this uses x/n where x is an random integer in range [1, n] (inclusive) to emulate a random float in range (0, 1]. In this program n is set to be 108.

ȷ8µ÷X


Try it online!

## Explanation

ȷ8     Literal 10^8.
÷   Divide by
X  random integer.


# Enlist, 3 bytes

ØXİ


Try it online!

Enlist beats Jelly! (TI-Basic not yet)

## Explanation

  İ    The inverse of...
ØX     a random float in [0, 1)


Of course this has nonzero probability of take the inverse of 0.

• Would the Enlist solution not fail if ØX returned 0? (Disclaimer: I don't know Enlist at all!) Dec 16, 2017 at 23:34
• @Shaggy your program must output a number greater than or equal to 1 with at least probability 0.999. The rest of the time it may crash (from the challenge rules) Dec 17, 2017 at 2:49

# IBM/Lotus Notes Formula, 13 bytes

1/(1-@Random)


Sample (10 runs) • Note to everyone: issacg has added some rules that allow some imprecisions, therefore most answers here are longer than necessary. Dec 16, 2017 at 15:17
• Not sure I could make this much shorter whatever rule changes are made :) Dec 16, 2017 at 22:41

# Java 8, 22 18 bytes

v->1/Math.random()


(Old answer before the rules changed: v->1/(1-Math.random()))

Try it here.

# JavaScript REPL, 15 19 bytes

1/Math.random()

• This will not yield correct results if Math.random() returns 0 Dec 15, 2017 at 7:33
• Probably 1/(1-Math.random())? Dec 15, 2017 at 7:35
• Fixed using u*29's solution
– l4m2
Dec 15, 2017 at 7:38
• You need _=> at the start to make this a function; snippets aren't allowed. Dec 15, 2017 at 8:07
• It's a full program using console running
– l4m2
Dec 15, 2017 at 8:19

# Pyt, 2 bytes

ṛ⅟


Explanation:

ṛ           Random number in [0,1)
⅟          Multiplicative inverse
Implicit print


Try it online!

# J, 9 Bytes

p=:%@?@0:


I couldn't figure out how to make it take no input, since p=:%?0 would evaluate immediately and remain fixed. Because of this its sort of long.

### How it works:

p=:        | Define the verb p
0:  | Constant function. Returns 0 regardless of input.
?@    | When applied to 0, returns a random float in the range (0,1)
%@      | Reciprocal


Evaluated 20 times:

    p"0 i.20
1.27056 1.86233 1.05387 16.8991 5.77882 3.42535 12.8681 17.4852 2.09133 1.82233 2.28139 1.58133 1.79701 1.09794 1.18695 1.07028 3.38721 2.88339 2.06632 2.0793


# Pyth, 4 bytes

c1O0


Try it here!

Alternative: c1h_O0.

• c1tOZ is 5, does it not work?
– Dave
Dec 15, 2017 at 12:56
• @Dave Doesn’t work, that returns negative values. I need 1-n not n-1 Dec 15, 2017 at 12:58
• Does Pyth not have a constant for 100? Dec 17, 2017 at 0:00
• @Shaggy I wish it did. Unfortunately, no constant for 100 AFAIK Dec 17, 2017 at 6:20

# Clean, 91 bytes

import StdEnv,Math.Random,System.Time
Start w=1.0/(1.0-hd(genRandReal(toInt(fst(time w)))))


Clean doesn't like random numbers.

Because the random generator (a Mersenne Twister) needs to be given a seed, I have to take the system timestamp to get something that differs passively per-run, and to do anything IO-related I need to use a whole Start declaration because it's the only place to obtain a World.

Try it online!