{"id":732,"date":"2016-05-29T11:03:47","date_gmt":"2016-05-29T15:03:47","guid":{"rendered":"http:\/\/www.clayford.net\/statistics\/?p=732"},"modified":"2020-09-13T10:17:57","modified_gmt":"2020-09-13T14:17:57","slug":"earthquake-data-and-benfords-law","status":"publish","type":"post","link":"https:\/\/www.clayford.net\/statistics\/earthquake-data-and-benfords-law\/","title":{"rendered":"Earthquake data and Benford’s Law"},"content":{"rendered":"

Much<\/a> has been<\/a> written<\/a> about Benford's Law, that weird phenonmenon where if you have a naturally occuring set of numerical data, 30% of the numbers will begin with 1, 18% will begin with 2, 12% will begin with 3, and so on. You might expect the distribution of leading digits to be uniformly distributed, but no, that just isn't the case. 60% of the time the leading digit is 1, 2, or 3. Though it may seem like some curious mathematical anomaly, Benford's Law has been used to detect fraud in elections<\/a> and accounting<\/a>. <\/p>\n

In this post let's use R to verify that earthquake data obtained from the USGS<\/a> does indeed follow Benford's Law. On May 28, 2016, I downloaded data for all earthquakes in the past 30 days<\/a>. The data has a number of fields, including earthquake magnitude, depth of the earthquake, location, time and many others. Benford's Law says those numbers should follow a particular distribution.<\/p>\n

The formula for Benford's Law is as follows:<\/p>\n

\\[P(d) = \\log_{10} \\left(1+\\frac{1}{d}\\right) \\]<\/p>\n

That says the probability that a digit d<\/em> occurs as the first number is equal the log base 10 of 1 + 1\/d<\/em>. We can quickly generate these for d<\/em> = 1 – 9 in R:<\/p>\n

log10(1 + (1\/(1:9)))\r\n<\/code><\/pre>\n
## [1] 0.30103000 0.17609126 0.12493874 0.09691001 0.07918125 0.06694679\r\n## [7] 0.05799195 0.05115252 0.04575749\r\n<\/pre>\n
And we can make a quick plot as follows:<\/p>\n
barplot(log10(1 + (1\/(1:9))), names.arg = 1:9, main = "Benford's Law")\r\n<\/code><\/pre>\n
\"plot<\/p>\n
So according to this law, if we look at the distribution of first digits in our earthquake data, we should see them closely follow this distribution. Let's find out!<\/p>\n
First we need to import the data. Thanks to the USGS, this data comes ready for analysis. All we need to do is read it into R:<\/p>\n
dat <- read.csv("all_month.csv")\r\nnrow(dat)\r\n<\/code><\/pre>\n
## [1] 8584\r\n<\/pre>\n
Over 8500 earthquakes in the past 30 days! A quick look at the magnitude shows us most of them are very small:<\/p>\n
summary(dat$mag)\r\n<\/code><\/pre>\n
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's \r\n##  -1.900   0.790   1.280   1.526   1.930   7.200      20\r\n<\/pre>\n
This also tells us we have some negative numbers (?) as well as some missing data, NA's<\/code>. We also have numbers that start with 0 and have a decimal, such as 0.79. (In this case, Benford would say the number “starts” with a 7, not a 0.) This means when we determine the leading digits, we'll need to ignore negative signs, leading 0s, decimals and missing values. <\/p>\n
Let's investigate at the mag<\/code> column. First we remove the NA's<\/code>. Below is.na(dat$mag)<\/code> generates a logical vector of TRUE and FALSE values. Adding an !<\/code> in front reverses the TRUE and FALSE values. Finally inserting the result in subsetting brackets returns only those values that are TRUE, ie not missing.<\/p>\n
digits <- dat$mag[!is.na(dat$mag)]\r\n<\/code><\/pre>\n
Next we extract the first digits. To help with this I use two R packages, magrittr<\/code> and stringr<\/code>. The magrittr<\/code> package allows us to “chain” commands together with the %>%<\/code> operator. The stringr<\/code> package provides the str_extract<\/code> function that allows us to extract phrases that follow a certain pattern. So what we have below can be translated as follows: <\/p>\n