# 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.

• – isaacg Dec 16 '17 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] – user202729 Dec 16 '17 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 – Luis Mendo Dec 17 '17 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. – Giuseppe Dec 15 '17 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. – plannapus Dec 15 '17 at 12:50
• @Mego that is incorrect. The Pareto distribution is absolutely continuous and thus has measure 0 for any number. – Therkel Dec 15 '17 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). – plannapus Dec 15 '17 at 14:59
• @plannapus Hmm... That's a good point. Floats make this entirely too complicated. – Mego Dec 15 '17 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. – Misha Lavrov Dec 15 '17 at 17:17
• TI-Basic developers knew in advance this challenge will happen...? It seems to win this time. – user202729 Dec 16 '17 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. – Misha Lavrov Dec 16 '17 at 14:55
• @user202729 Yeah, actually I just time traveled a bit to get it all done. Definitely worth it for these upvotes. – Timtech Dec 16 '17 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. – user202729 Dec 15 '17 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. – Giuseppe Dec 15 '17 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. – KSmarts Dec 15 '17 at 15:27
• @KSmarts The "trivial" answer in R wasn't used by anybody actually: actuar::rpareto(1,1,1), because it is longer :) – plannapus Dec 15 '17 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. – plannapus Dec 15 '17 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 '17 at 9:23
• @Neil please wait a moment while I try to understand what your code does... :-) – Charlie Dec 15 '17 at 9:25
• Actually I don't need MapAssignRight any more, 10 bytes! works. – Neil Dec 15 '17 at 9:27
• @Neil assimilation of your code completed! Answer edited. :-D – Charlie Dec 15 '17 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)/) – Laikoni Dec 15 '17 at 11:55
• And as functions are allowed, I'd say you can drop the main=. – Laikoni Dec 15 '17 at 11:57
• Appranetly the range is [0,1) according to this answer – flawr Dec 15 '17 at 13:30
• @flawr Whoops, you're right! I forgot how floats work temporarily. – Mego Dec 15 '17 at 13:32
• Well anyway, thanks for commenting, I would not have had any idea:) – flawr Dec 15 '17 at 13:35

# Excel, 9 bytes

=1/rand()


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

• Also works in LibreOffice Calc :) – ElPedro Dec 16 '17 at 7:59
• You can change this to google sheets for -1 Bytes (=1/Rand() – Taylor Scott Dec 18 '17 at 15:04

# Mathematica, 10 bytes

1/Random[]


-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. – Mego Dec 15 '17 at 11:42
• @Mego where exactly did you find that info? – J42161217 Dec 15 '17 at 11:45
• @Mego what is the probability to get 0? – J42161217 Dec 15 '17 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. – user202729 Dec 15 '17 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. – Kelly Lowder Dec 15 '17 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. – user202729 Dec 16 '17 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. – Erik the Outgolfer Dec 15 '17 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. – user202729 Dec 16 '17 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. – user202729 Dec 16 '17 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. – user202729 Dec 16 '17 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. – user202729 Dec 16 '17 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!) – Shaggy Dec 16 '17 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) – user202729 Dec 17 '17 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. – user202729 Dec 16 '17 at 15:17
• Not sure I could make this much shorter whatever rule changes are made :) – ElPedro Dec 16 '17 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 – Mr. Xcoder Dec 15 '17 at 7:33
• Probably 1/(1-Math.random())? – user202729 Dec 15 '17 at 7:35
• Fixed using u*29's solution – l4m2 Dec 15 '17 at 7:38
• You need _=> at the start to make this a function; snippets aren't allowed. – Shaggy Dec 15 '17 at 8:07
• It's a full program using console running – l4m2 Dec 15 '17 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 '17 at 12:56
• @Dave Doesn’t work, that returns negative values. I need 1-n not n-1 – Mr. Xcoder Dec 15 '17 at 12:58
• Does Pyth not have a constant for 100? – Shaggy Dec 17 '17 at 0:00
• @Shaggy I wish it did. Unfortunately, no constant for 100 AFAIK – Mr. Xcoder Dec 17 '17 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!