# Using a power law to model the number of billionaires

Here’s a brief tutorial about using power laws to numerically model the number of billionaires in a society. We’ll build the model in R.

To get started, let’s install the `poweRlaw` package:

``````# install the powerRlaw package
install.packages("poweRlaw")
``````

The `poweRlaw` package comes with many useful functions for working with power laws. Today, we’ll use the `rplcon` function, which generates random numbers drawn from a continuous power law.

The `rplcon` function takes 3 inputs:

1. `n`: the size of our random sample of numbers
2. `xmin`: the lower threshold for our power law
3. `alpha`: the power-law exponent

For example, here’s how we’d sample 10 numbers from a power-law distribution that has a lower threshold of 100 and a power-law exponent of 3:

``````# example
rplcon(n = 10, xmin = 100, alpha = 3 )
``````

Here’s the type of output we’d get

``````> rplcon(n = 10, xmin = 100, alpha = 3 )

 108.6708 155.1333 112.7486 116.7339 137.0705 230.2303 127.7278 133.8075 122.2766 131.7206
``````

### Counting billionaires

Let’s suppose that our power law describes a distribution of wealth. Given this distribution, we want to know the portion of people who are billionaires. Here’s how we do it.

As, before, we’ll use the `rplcon` function to generate a power-law sample. Let’s set the minimum wealth to `xmin = 1000` and use a power-law exponent of `alpha = 2`. And let’s suppose we have a population of 10 million people (`n = 1e7`). Here’s a sample of the distribution of wealth:

``````# simulate a power-law distribution of wealth
wealth = rplcon( n = 1e7, xmin = 1000, alpha = 2)
``````

Now that we have simulated data for wealth, we can count the number of billionaires. To isolate our billionaires, we’ll use the `which` function. We’ll take our `wealth` vector and ask which elements have values greater than or equal to 1 billion. (In scientific notation, 1 billion is `1e9`.)

``````# isolate billionaires
billionaires = which( wealth >= 1e9 )
``````

Next, we’ll calculate the number of billionaires per capita. To do that, we’ll use the `length` function, which tells us how many elements are in a given vector.

To count the number of billionaires, we’ll use `length(billionaires)`. And to count the whole population, we’ll use `length(wealth)`. The number of billionaires per capita is then the ratio of the two numbers:

``````# count billionaires
n_billionaires = length( billionaires )

# count population
n_population = length( wealth )

# calculate the number of billionaires per capita
billionaires_per_capita = n_billionaires / n_population
``````

Doing the math, I calculate a billionaire density of about 0.8 per million:

``````> billionaires_per_capita
 8e-07
``````

Since we’re working with random numbers, we’ll get a slightly different billionaire density each time we take the sample. Still, there should be a pattern between our power-law parameters and the number of billionaires. To explore the pattern, play around with different input values and see what you find.