Investigating Power Laws with Mathematica

Fiona Maclachlan

Manhattan College

Abstract

Power-law distributions are found in a broad range of disciplines.  In the late nineteenth century, Vilfredo Pareto identified a power law for the distribution of income.  More recently, power laws have been discovered in the degree distributions of socially constructed networks like the World Wide Web, and have been associated with phenomena characterized by preferential attachment. They have also been linked with the idea of self-organized criticality and have been observed in the size distributions of many natural phenomena, such as sandpiles and earthquakes.  Recent empirical studies of economic data have turned up power law behavior in the return distribution of financial assets, and in the size distributions of firms and market shares.  The latter work has been picked up in the marketing literature and has even found its way into popular business books like The Long Tail.

This presentation illustrates ways in which Mathematica can be used to investigate issues related to power laws.  I show how the non-linear regression capability can be used to investigate the goodness of fit of a power-law distribution to a data set, and to compare the fit with that of other similarly shaped distributions, such as the lognormal.  In addition, I show how simple computational models can be constructed to investigate the extent to which power laws arise from randomness in an underlying process.  For example, the hypothesis that Zipf’s Law in word frequency would be generated by monkeys typing can be investigated with a simulation, as can a model of market shares based on the idea of random division.

Introduction

<< Graphics`Graphics`
<< Statistics`ContinuousDistributions`
<< Statistics`NonLinearFit`

Size versus Rank Distributions

The Pareto distribution is an example of a power law size distribution.  We consider a population (say, of families) and then inquire about the size of a quantity associated with each member of the population (say, family income).  The size of the quantity is the independent variable and the dependent variable is the relative frequency of members of the population associated with that size.  In their cumulative form, size distributions indicate the relative frequency of members of the population associated with a given size and below.  The cumulative distribution function for the Pareto distribution is given by:

CDF[ParetoDistribution[k, α], x]

where k is the minimum size of the quantity in question, and α is the shape parameter for the distribution

With[{k = 1, α =  .9},  
        Plot[CDF[ParetoDistribution[k, α], x],  {x, 1, 100},  
          AxesLabel -> {Size,  CDF},  
          PlotStyle -> Thickness[ .01],  
          ImageSize -> 500]];

With a rank distribution each member of the population is ranked according to the size of some quantity.  Families could be ranked according to their income from richest to poorest, or books could be ranked according to their level of sales from best selling to worst selling.  The rank becomes the independent variable and the dependent variable is the size of the quantity associated with that rank.  A two-parameter power law rank distribution is given by:

(c/(r^β))

where c is the size associated with the highest ranking member of the population, r is the rank, and β is the shape parameter.  Zipf's Law, as it is commonly understood, is a special case of a power law rank distribution in which β is 1. The Zipf-Mandelbrot distribution adds a third parameter to give a more general,  (c/(((r + q))^β)).

In some cases, one might want to consider the relative frequency distribution of ranked data.  This discrete distribution is also know as the zeta distribution and "will be implemented in a future version of Mathematica as ZipfDistribution[rho]".

ListPlot[
        With[{c = 1000, β = 1.1}, Table[(c/(r^β)), {r, 1, 100}]], 
        AxesLabel -> {"rank", "size"},  
        PlotStyle -> PointSize[ .015],  
        ImageSize -> 500];

Log-Log Plots

The command LogLogListPlot can be applied to data to investigate if a power law is present, since a power law distribution, whether in the size or rank form, will show up as a straight line.  For example, let's generate data consistent with the Pareto distribution and apply the command:

Transpose[{Range[12000,  100000,  100], 
          1 - CDF[ParetoDistribution[12000,   .65], 
              Range[12000, 100000,  100]]}];

LogLogListPlot[%,  Ticks -> None];

Here are LogLogListPlots of data from two similar distributions.

Transpose[{Range[12000,  100000,  100], 
          1 - CDF[LogNormalDistribution[10,   .9], 
              Range[12000,  100000,  100]]}];

LogLogListPlot[%,  Ticks -> None];

Transpose[{Range[12000,  100000,  100], 
          1 - (1/(1 + ((((44 + 27))^3)/(((Range[12000,  100000,  100] + 27))^3))))}];

LogLogListPlot[%,  Ticks -> None];

Nonlinear Curve Fitting

To illustrate nonlinear curve fitting consider the following data which represent the proportion of households falling below various levels of income in 2004.

incomepoints = {{15000, 0.154`}, {25000, 
            0.28300000000000003`}, {35000, 0.402`}, {50000, 0.55`}, {75000, 
            0.7330000000000001`}, {100000, 
            0.8430000000000001`}};

\(\(ListPlot[incomepoints, \ PlotStyle -> PointSize[ .02], \ 
        AxesOrigin -> {0, 0}, 
        AxesLabel -> TraditionalForm /@ \ {"income", \ "cdf"}, \ 
        Ticks -> {{{15000, "$15,000"}, \ {50000, "$50,000"}, \ {100000, \ "$100,000"}}, \ Automatic}, \ 
        ImageSize -> 600];\)\)

Let's fit these data points to different distributions and then see how well the fitted distributions "connect the dots."

Pareto

CDF[ParetoDistribution[k, α], x]

NonlinearRegress[incomepoints,  
        CDF[ParetoDistribution[k, α],  x], x, {k, α}, 
        RegressionReport -> BestFitParameters][[1, 2]]

ParetoCDF = CDF[ParetoDistribution[k, α],  x] /. %

DisplayTogether[
        ListPlot[incomepoints, 
          AxesLabel -> {"income",  "cdf"},  
          PlotStyle -> PointSize[0.02],  AxesOrigin ->  {0, 0},  
          PlotLabel -> "Pareto"], 
        Plot[ParetoCDF,  {x,  5000,  100000}],  
        Ticks -> {None,  Automatic},  
        ImageSize -> 600];

Lognormal

CDF[LogNormalDistribution[μ, σ], x]

NonlinearRegress[incomepoints,  
        CDF[LogNormalDistribution[μ, σ],  x], 
        x, {{μ,  5,  15}, {σ,  .5, 1.5}}, 
        RegressionReport -> BestFitParameters][[1, 2]]

LogNormalCDF = 
      CDF[LogNormalDistribution[μ, σ],  x] /. %

DisplayTogether[
        ListPlot[incomepoints, 
          AxesLabel -> {"income",  "cdf"},  
          PlotStyle -> 
            PointSize[0.02],   PlotLabel -> "LogNormal"], 
        Plot[LogNormalCDF,  {x,  15000,  100000}],  
        Ticks -> {None,  Automatic},  
        ImageSize -> 600];

Beaman

The Beaman distribution was used at duPont in the 1970's to model sales volume of products at various price points.  It was found to give a better fit than the lognormal.  There are two parameters: a represents the minimum price point and b represents the median price point.

Beaman = (1/(1 + ((((b - a))^3)/(((x - a))^3))));

NonlinearRegress[
        incomepoints,  (1 + (b - a)^3/(x - a)^3)^-1,  
        x,   {{b, 25000, 45000}, {a,  -25000, 0}},  
        RegressionReport -> BestFitParameters][[1, 2]]

BeamanCDF = Beaman /. %

DisplayTogether[ListPlot[incomepoints, 
            AxesLabel -> {"income",  "cdf"},  
            PlotStyle -> PointSize[0.02],  
            AxesOrigin ->  {0, 0},  PlotLabel -> "Beaman"],Plot[BeamanCDF,  {x,  5000,  100000}], Ticks -> {None,  Automatic}, ImageSize -> 600];

Two way comparison

DisplayTogether[
        ListPlot[incomepoints, AxesLabel -> {"income",  None},  
          PlotStyle -> PointSize[0.01],  AxesOrigin ->  {0, 0},  
          PlotLabel -> "Beaman in Blue & Lognormal in Red"],  
        Plot[LogNormalCDF,  {x,  5000,  100000},  
          PlotStyle -> {Thickness[ .003], 
              RGBColor[1, 0, 0]}], 
Plot[
          BeamanCDF,  {x,  5000,  100000},  
          PlotStyle -> {Thickness[ .003], 
              Hue[0.7]}], 
Ticks -> {None,  
            Automatic},  ImageSize -> 600];

Leading Digits

An interesting consequence of [Zipf's Law] is easily verified by perusing any of the lists mentioned  ... There is an overwhelming preponderance of 1s, followed by 2s, and so forth.  The percentage of population figures with initial 9s, is extremely small.  That behavior of the first digit is predicted by the rule, which, if exactly obeyed, would give a proportion of initial 1s to intial 9s of 45 to 1.
	Murray Gell-Mann, The Quark and the Jaguar, pp. 93-94.

LDP[list_List] := 
      N[Table[ Count[First /@ First /@ RealDigits[list], d], {d, 9}]/
          Length[list]]

N[Table[1/x,  {x,  10000}]];

LDP[%]

%[[1]]/%[[9]]

Table[N[Table[i/x,  {x,  10000}]],  {i, 15}];

LDP /@ %

#[[1]]/#[[9]] & /@ %

Monkeys Typing

The Theoretical Distribution

Let a monkey choose randomly from three letters, A, B & C, and a space.
What is the probability at any point in his random typing that he will type an m-letter word?
There are two possibilities to consider:
(1) He has just typed a space.  That will be true 25% of the time on average.
(2) He has just typed a letter.  That will be true 75% of the time on average.
 
If he has just typed a space, then the probability of his typing a particular m-letter word (1/(4^(m + 1))).  
If he has just typed a letter, then the probability of his typing a particular m-letter word is (1/(4^(m + 2))).  
So the probabilty of a particular m-letter word is (1/4)×(1/(4^(m + 1))) + (3/4)×(1/(4^(m + 2))) .
But there are (3^m) different m-letter words.
So the probability of the monkey typing any m-letter word is:

((1/4)×(1/(4^(m + 1))) + 
          (3/4)×(1/(4^(m + 2)))) (3^m)

Simplify[%]

probABC[m_] := 7 (3^m) (4^(-3 - m))

Here are the probabilities of words of sizes 1 through 20:

Table[probABC[m], {m, 20}] // N

Here are the cumulative probabilities:

FoldList[Plus, 0,  %]

ListPlot[%];

Ranking  Ranks 1 to 3 are taken up by the one letter words.  Then the two-letter words occupy from rank 4 up to

3 + (3^2)

The three-letter words occupy ranks 13 up to

3 + (3^2) + (3^3)

When different entities tie for the same rank, it's not clear the best way to represent the ranking graphically. The one letters "A", "B" and "C" tie for first place. One possibility is to give them all a rank of 3, meaning they are all in the top three.  The 9 two-letter words would then be given a rank of 12, the 27 three-letters words a rank of 39, and so on.

\(\(BarChart[
        Transpose[{Table[7\ 3\^m\ 4\^\(\(-3\) - m\)\ , \ {m, \ 1, \ 20}], \ 
            FoldList[Plus, 3, 3^Range[2, 20]]}], \ 
        Ticks -> {{{1, 3}, \ {2, 12}, {3, 
                39}, \ {4, "..."}, \ {5, \ "and"}, \ {6, "so"}, {7, \ "on"}, \ {8, "..."}}, \ Automatic}, \ 
        AxesLabel -> {"Word Rank", \ "Probability"}, \ 
        ImageSize -> 600];\)\)

A loglog plot gives the expected straight line.

LogLogListPlot[
        Transpose[{FoldList[Plus, 3, 3^Range[2, 100]], 
            Table[7 (3^m) (4^(-3 - m)) ,  {m,  1,  100}]}],  
        AxesLabel -> {"Log Rank", "Log Probability"}, 
        PlotJoined -> True,  PlotStyle -> Thickness[ .005],  
        ImageSize -> 600];

Distribution from a Simulation

A simulation confirms the above results.

MonkeysTypingABC[n_] := 
      Module[{keys = Append[CharacterRange["a", "c"], "_"], 
          typing,  spaces,  words}, 
typing = 
          Table[keys[[Random[
                  Integer,  {1, 4}]]],  {n}]; 
spaces = 
          Flatten[Position[typing, "_"]]; 
words = 
          Sort[Append[
              Prepend[Table[
                  typing[[Range[spaces[[i]] + 1, 
                        spaces[[i + 1]] - 1]]],  {i, 1,  
                    Length[spaces] - 1}],  
                typing[[Range[1,  spaces[[1]] - 1]]]],  
              typing[[Range[Last[spaces] + 1,  
                    Length[typing]]]]]]]

dataABC = MonkeysTypingABC[100000];

sizesABC = Length /@ dataABC;

Length[sizesABC]

First[Position[Table[Count[sizesABC,  i], {i, 50}], 0]]

wABC = Table[
          Count[sizesABC,  i], {i, 50}][[Range[1, 26]]]

\(\(BarChart[
        Transpose[{wABC[\([Range[15]]\)], 
            FoldList[Plus, 3, 3^Range[2, 15]]}], 
        Ticks -> {{{1, 3}, \ {2, 12}, {3, 
                39}, \ {4, "..."}, \ {5, \ "and"}, \ {6, "so"}, {7, \ "on"}, \ {8, "..."}, \ {15, "21,523,359"}}, \ Automatic}, \ 
        AxesLabel -> {"Word Rank", \ "Frequency"}, \ 
        ImageSize -> 600];\)\)

LogLogListPlot[
        Transpose[{FoldList[Plus, 3, 3^Range[2, 15]], 
            wABC[[Range[15]]]}],  
        AxesLabel -> {"Log Rank", "Log Frequency"},  
        PlotStyle -> Thickness[ .005],  PlotJoined -> True,  
        ImageSize -> 600];

The Broken Stick Rule and the Long Tail

The broken stick rule gives the expected proportional sizes of segments of a line broken randomly into n parts.  Ranked by size largest to smallest, the size of the ith segment is given by

(((∑_(j = 0)^(n - i))(1/(n - j)))/n)

More generally, the expected sizes of the parts of a line broken into n parts is given by:

BrokenStickParts[n_] := 
      Table[1/n Sum[1/(n - j), {j, 0, n - i}], {i, n}]

For n = 10,

BrokenStickParts[10] // N

ListPlot[BrokenStickParts[10],  AxesOrigin -> {0, 0},  
        PlotStyle -> PointSize[ .03],  
        AxesLabel -> {"Rank",  "Proportional Size"}];

If we think of the line segments as the market shares of different types of products and display them in a bar chart format, we can generate something that looks like the Long Tail logo.  The red portion represents the head and is defined as those products sold through traditonal distribution outlets, while the orange part is the long tail of products distributed through the Internet.

bsp1K = BrokenStickParts[1000];

BarChart[bsp1K,  Ticks -> None,  BarEdges -> False,  
        BarStyle -> (If[# < bsp1K[[200]],  RGBColor[1,  .5, 0], 
                RGBColor[1, 0, 0]] &),  
        AxesLabel -> {"Products",  "Sales"}];

A log-log plot reveals that the distribution does not appear to follow a power law at the tail end.

LogLogListPlot[bsp1K];

But by varying n we see that the larger the number of products, the greater is the range over which a power law is apparent.

loglogplot[n_] := 
      LogLogListPlot[Transpose[{Range[n], BrokenStickParts[n]}],  
        PlotLabel  -> n Products]

Table[loglogplot[n],  {n,  200,  1400,  200}];

From www.longtail.com :