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

`n`

: the size of our random sample of numbers`xmin`

: the lower threshold for our power law`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 )
[1] 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
[1] 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.