I am attempting to identify possible outliers in data which is skewed to the right and I assume it is Poisson distributed. I am a novice in all things statistics, and the following may be utterly erroneous. However, I am eager to learn.
I have scoured Cross Validated and Stack Overflow for ideas in detecting outliers in situations like mine, but I couldn't find any instances where someone attempted to write an R script for their project to find out outliers in skewed, Poisson distributed data.
My actual data is shown below as the vector parktimes, (n=5222). It is the results of a survey where respondents answered how long it took them (in minutes) to park their car in a postal code area in Helsinki, Finland. Respondents could answer for multiple postal code areas at the same time, leaving the data with some identical timestamps with different values for different postal code areas. Most people reported finding a parking place almost instantaneously, making the data skewed to the right. The allowed sequence here was 0-99, but 99 minutes to find a parking place in Helsinki seems unbelivable and yet someone answered with that value for multiple postal code areas. I would like to find a statistical solution to remove these improbable values if they are indeed outliers. The code below does not provide the exact timestamps to be more concise, substituting to the index.
Here is a histogram of parktime values with ggplot: ggplot(thesisdata, aes(parktime)) + geom_histogram(color = "black", binwidth = 5)
Using this source and this source I have written a R script that I think detects outliers in my data. Simplified, my attempted outlier detection is as follows:
- Import data
- Apply Anscombe transform to data.table column
parktimelike so:anscombe_parktime <- 2.0 * sqrt(parktime + 3.0 / 8.0) - Calculate the probability of observing a point under a Poisson distribution:
ppois(anscombe_parktime, mean(anscombe_parktime) - Plot results
With Anscombe transform (y axis is parktime):
Without Anscombe transform:
Is this a legitimate way to search for outliers in my data? Can Anscombe transform be used in this way to wrangle the data? Is my data even applicable for this kind of analysis?
My code:
library(ggplot2)
library(data.table)
library(outliers)
parktimes <- c(99,5,0,1,10,99,99,1,1,3,1,1,2,5,2,2,2,5,10,5,2,2,0,1,1,1,5,3,5,5,
1,0,0,5,1,0,0,2,2,0,5,10,1,1,1,5,5,3,10,1,1,1,1,0,10,2,10,7,10,7,
3,3,13,1,3,1,1,1,4,4,1,2,3,1,1,1,1,1,1,2,1,1,2,3,0,7,8,3,3,3,5,4,
25,0,10,0,10,6,3,0,0,1,2,1,0,0,0,0,0,0,3,1,0,1,2,1,0,1,5,5,5,3,0,
0,0,0,2,1,3,0,1,5,5,5,2,0,2,0,5,15,3,4,3,4,2,5,1,10,10,2,0,1,1,1,
0,0,1,0,10,5,15,1,10,0,0,2,1,5,1,1,2,2,3,1,1,1,1,4,4,1,3,3,1,3,1,
2,1,0,1,2,2,5,1,2,1,3,5,1,1,1,1,5,4,5,2,15,15,2,5,2,5,8,2,8,5,5,2,
0,1,3,2,1,1,1,1,1,1,1,1,10,3,1,8,10,10,12,5,5,3,6,4,2,1,3,2,0,0,1,
3,1,1,1,1,2,1,3,1,1,2,1,1,3,1,1,1,3,2,1,1,2,2,1,4,1,1,1,1,2,1,2,3,
4,1,2,1,2,10,1,0,0,3,3,10,1,4,0,2,5,5,1,4,0,5,1,1,1,3,0,1,5,1,1,1,
1,1,1,5,5,5,5,5,10,20,1,1,1,0,0,0,0,1,0,2,0,2,2,2,0,1,1,1,2,2,2,0,
1,0,1,2,1,5,0,0,10,1,2,1,2,2,3,2,3,1,1,2,5,2,1,5,5,2,10,2,4,0,5,0,
1,1,5,1,2,5,1,1,3,4,1,6,6,5,2,10,10,10,60,7,1,15,10,0,5,15,1,0,2,
0,0,0,2,1,2,3,3,2,2,3,3,2,3,1,3,5,1,2,1,3,10,1,1,1,1,5,3,1,6,12,5,
7,6,5,2,0,3,1,5,10,30,45,45,30,45,0,0,0,0,5,5,0,3,5,2,5,10,10,2,5,
10,2,1,30,5,2,2,7,1,1,2,4,5,5,1,1,1,5,2,2,2,2,1,5,0,1,3,5,5,1,2,
15,10,0,1,10,8,10,25,5,10,5,12,20,7,12,2,5,2,10,3,10,5,5,5,5,5,7,
3,7,3,6,9,7,1,1,10,10,1,1,1,1,2,1,15,30,1,10,5,20,1,10,1,35,10,0,
5,25,35,10,1,5,5,10,20,5,5,5,10,10,15,2,2,1,1,1,1,1,3,5,5,5,1,1,5,
10,10,15,15,25,20,5,15,5,0,5,5,2,5,3,10,2,5,5,1,15,8,4,6,5,15,20,
20,20,15,15,15,30,15,10,5,5,10,10,10,10,5,5,0,10,1,5,1,2,0,2,2,5,
10,15,3,15,3,4,3,2,1,3,4,5,4,2,10,1,1,1,1,5,1,10,5,5,10,5,1,5,7,
10,10,5,10,5,1,2,15,10,1,10,10,15,10,10,5,2,2,2,5,5,10,5,5,2,5,5,
2,5,10,10,20,5,1,2,2,5,2,5,1,1,15,10,20,15,4,15,15,5,15,5,0,5,1,0,
0,5,6,7,1,3,2,3,2,0,10,15,10,10,3,30,10,30,5,10,20,10,0,1,10,1,2,
2,1,1,0,1,10,0,10,15,5,5,10,5,8,4,10,10,3,3,5,5,1,4,0,15,2,10,10,
2,2,10,2,5,10,1,1,1,1,1,2,2,1,1,1,2,1,1,2,2,8,4,5,1,3,5,10,1,2,1,
2,1,0,1,0,8,10,3,15,0,0,0,1,2,0,1,0,5,2,10,5,2,10,5,1,1,0,2,5,1,1,
1,3,2,3,2,2,6,9,9,9,8,2,9,10,5,10,1,15,10,4,5,5,5,1,7,1,10,2,2,8,
2,2,7,1,1,10,2,8,10,2,5,5,4,3,5,5,8,6,8,4,2,10,15,4,8,3,6,5,5,6,0,
1,10,15,10,3,5,1,8,10,7,1,1,2,5,10,10,15,0,2,5,5,5,10,3,5,1,4,1,1,
14,24,5,5,15,3,0,5,0,5,5,6,0,1,2,1,1,4,1,10,2,5,1,1,5,8,5,10,19,0,
3,5,2,5,0,2,2,5,1,2,2,5,1,2,2,1,5,2,2,1,1,5,15,1,1,1,5,1,1,7,5,3,
5,1,10,1,1,2,4,1,1,2,4,2,1,0,1,2,1,10,5,10,3,15,1,1,15,5,10,1,1,
1,10,20,20,5,1,10,15,1,10,5,1,5,5,5,5,5,20,20,5,1,5,5,10,5,5,20,
5,15,15,10,2,0,0,3,2,5,1,2,1,0,3,0,5,1,1,1,5,1,1,5,10,10,0,1,1,1,
1,5,5,10,5,5,1,8,10,10,10,2,3,5,3,15,3,5,0,0,0,1,1,1,1,0,1,1,1,1,
1,1,1,1,0,1,2,1,1,1,1,0,1,1,1,10,15,10,10,10,20,5,3,1,7,7,5,20,1,
2,5,5,5,5,0,7,1,5,1,1,1,1,1,1,5,1,3,1,3,2,2,5,0,45,5,10,10,5,10,5,
1,2,5,2,5,2,1,1,5,2,15,20,10,35,5,5,5,5,10,20,15,15,1,2,5,5,2,2,3,
5,1,1,10,10,1,1,1,0,2,3,7,2,1,2,2,1,2,3,4,2,1,28,20,1,5,5,8,2,0,0,
3,8,1,3,2,15,15,15,8,4,20,0,2,2,5,1,1,5,7,5,0,5,1,15,2,2,12,10,6,
15,0,2,4,5,5,10,1,1,1,1,2,6,2,1,0,1,3,3,5,3,6,8,2,60,90,15,3,10,1,
5,3,1,6,1,2,2,7,3,3,15,25,10,5,10,8,7,1,1,1,5,3,5,1,2,5,0,1,2,1,2,
1,1,1,1,5,2,25,20,0,0,4,1,5,1,15,10,1,1,3,1,1,5,6,5,1,14,15,6,15,
8,7,1,4,8,5,2,1,0,1,1,1,2,6,3,5,5,2,8,4,1,10,5,4,8,3,3,3,1,3,2,1,
2,3,1,2,6,3,4,6,2,8,1,5,5,1,2,6,1,3,1,2,0,1,5,3,1,3,5,3,5,7,2,5,
15,2,2,5,1,3,5,7,10,5,5,10,10,10,5,2,10,7,20,2,5,10,5,2,2,4,3,5,
2,1,10,2,5,20,5,20,5,1,0,0,2,2,1,5,30,99,10,1,5,10,10,5,2,10,1,5,
3,2,10,4,1,5,5,2,10,5,1,2,10,4,5,3,2,2,1,0,2,55,0,3,10,3,20,5,20,
5,5,3,5,5,5,3,1,5,10,10,5,1,10,0,2,5,1,2,20,5,2,10,5,5,8,1,5,10,2,
5,1,3,1,2,3,5,1,1,5,5,20,5,5,15,1,5,1,5,1,5,99,99,20,99,99,99,99,
2,2,2,1,2,3,1,2,2,1,2,1,2,1,1,2,2,2,1,2,1,1,1,1,1,1,1,1,4,1,1,1,
2,2,3,2,3,2,1,2,3,2,2,2,2,5,2,5,5,3,2,3,2,3,3,5,2,5,5,1,1,1,1,3,2,
2,3,3,2,10,5,1,3,3,0,2,10,5,2,2,3,2,5,3,2,15,5,7,10,1,5,5,2,2,3,2,
2,10,10,15,2,5,15,5,10,6,3,5,2,5,5,5,8,4,4,5,5,4,2,2,5,2,5,5,0,5,
2,5,5,0,0,0,5,10,5,10,1,5,5,1,1,3,20,20,0,0,3,0,2,1,2,1,1,2,1,1,8,
2,2,5,5,0,3,20,6,1,2,4,1,15,2,4,5,5,2,5,10,5,1,1,1,3,2,1,2,3,4,6,
5,10,5,5,2,10,10,10,10,10,10,0,10,10,5,10,10,5,5,5,10,10,10,5,1,1,
3,10,5,5,1,1,0,0,2,10,10,5,5,5,2,2,5,2,10,5,10,1,10,3,2,1,3,2,3,3,
5,1,1,2,6,3,5,5,10,5,3,5,5,10,5,4,5,3,3,1,2,1,3,5,1,1,1,1,1,2,2,5,
6,2,4,2,2,2,5,10,2,2,3,3,2,1,2,2,4,2,1,5,10,5,1,1,3,0,5,3,5,5,1,2,
2,5,3,1,10,2,5,3,10,10,3,10,5,2,3,10,0,2,3,2,1,0,10,2,0,1,2,4,2,2,
5,2,7,0,0,5,7,7,5,1,5,10,5,1,3,4,6,5,2,15,5,4,10,3,2,10,3,3,4,10,
2,8,5,0,2,1,1,3,3,1,1,1,1,1,1,2,1,3,1,1,10,2,1,1,0,1,0,10,30,5,15,
5,5,10,5,5,5,5,1,0,0,0,7,1,5,5,2,1,2,5,20,30,15,15,1,0,0,0,0,2,5,
0,0,0,3,0,0,2,5,0,0,4,0,1,2,3,0,4,3,1,1,3,20,5,5,10,10,15,15,10,5,
3,1,4,10,10,2,10,2,1,5,5,2,2,2,1,1,1,1,1,3,2,2,3,1,7,1,1,3,1,1,3,
3,2,5,2,2,5,5,2,1,3,1,1,1,2,5,5,1,10,2,3,5,1,5,10,0,5,5,0,0,3,3,1,
1,1,15,3,15,2,2,5,1,5,0,1,1,2,2,1,4,5,1,3,2,10,3,5,7,10,3,3,3,4,3,
2,2,0,0,1,1,4,1,3,1,1,3,5,1,10,15,3,3,1,1,5,5,2,10,2,5,5,7,5,8,7,
6,4,5,4,4,2,8,10,9,15,8,5,0,0,2,5,0,5,1,3,2,5,20,10,30,10,30,15,
10,15,15,10,10,10,10,5,15,1,1,2,0,1,4,5,5,0,2,5,4,1,2,0,0,1,2,1,5,
6,1,1,3,1,1,1,1,3,5,10,5,5,2,5,0,1,3,0,3,5,5,15,10,10,0,5,10,5,2,
10,5,2,10,5,2,5,10,5,1,20,5,15,5,5,5,5,5,5,5,10,10,5,5,5,5,5,10,5,
0,0,10,10,5,5,1,25,5,1,1,5,1,2,1,1,1,2,3,10,1,30,10,20,10,20,5,15,
10,10,15,25,15,1,0,7,2,1,0,3,3,4,15,5,15,10,3,10,5,3,2,1,1,3,1,3,
25,0,10,5,7,5,20,10,18,20,5,2,1,1,1,1,1,1,2,2,5,2,2,5,5,10,5,10,10,
3,2,1,1,8,5,2,2,5,5,5,1,5,5,2,15,0,0,2,10,5,1,1,2,0,5,1,5,5,5,2,10,
5,0,5,5,1,4,1,0,4,0,3,4,1,1,0,0,3,5,1,2,1,10,5,5,2,2,3,0,20,2,5,1,0,
3,1,5,5,15,5,5,5,2,0,3,3,0,0,5,5,5,1,2,3,1,10,10,1,1,3,1,0,5,0,10,5,
10,10,10,0,2,3,2,0,10,2,15,2,6,2,10,5,2,3,10,3,5,3,3,5,3,5,4,3,10,5,
5,5,10,2,4,5,6,8,5,5,4,2,15,4,15,5,10,5,5,2,1,1,1,2,3,2,3,4,5,0,10,
15,5,5,1,3,15,1,10,3,1,10,5,5,5,3,7,8,1,10,3,3,0,0,7,15,15,5,3,15,
2,10,1,7,5,20,2,10,5,1,1,1,2,1,5,15,15,5,1,5,7,9,3,2,5,5,15,10,20,
0,20,25,5,15,10,2,3,2,2,5,2,1,5,5,6,6,1,1,3,1,1,3,3,10,2,20,20,5,5,
4,0,30,20,5,15,0,10,10,1,6,3,1,2,2,10,2,1,1,1,0,10,2,2,5,5,4,5,16,
2,1,10,30,15,5,3,2,10,10,1,3,1,3,2,2,10,2,1,3,1,1,1,1,3,3,5,7,5,3,
10,5,1,10,2,2,1,1,5,1,2,3,2,2,2,5,1,1,1,10,2,1,1,1,3,1,6,1,3,5,1,
3,10,10,0,0,0,0,0,15,10,10,15,1,7,3,5,5,1,5,10,6,2,4,2,2,1,1,4,2,
1,2,4,1,3,3,1,1,1,2,1,2,2,2,4,1,1,1,2,2,1,2,1,2,4,4,2,1,8,3,1,3,2,
5,5,2,2,4,3,3,1,1,1,2,1,2,2,1,2,3,2,2,5,0,0,0,3,5,1,1,1,1,2,2,5,5,
5,0,4,1,1,5,10,5,5,3,1,3,3,4,5,1,3,2,3,3,3,2,3,2,4,5,3,5,2,5,5,6,1,
3,7,4,30,3,1,1,3,15,10,2,1,5,1,1,2,1,3,1,1,2,3,1,1,1,1,1,2,1,1,10,
2,2,2,2,5,1,25,30,10,3,15,5,5,30,20,20,40,35,20,10,5,0,5,2,15,20,
2,7,10,2,2,1,15,5,0,20,10,0,10,15,1,3,1,0,1,2,1,0,3,5,2,4,7,6,7,4,
2,2,1,2,2,2,2,6,1,8,6,5,2,5,4,2,5,2,3,3,1,2,1,1,3,2,3,15,2,2,1,4,
1,2,1,1,1,2,1,2,1,1,2,2,1,2,1,1,1,1,1,2,10,2,5,10,20,10,5,10,10,5,
20,15,10,5,20,20,15,10,25,15,20,15,10,15,2,15,5,5,3,1,5,1,5,2,1,0,
5,4,1,2,1,3,5,5,5,5,10,8,1,5,10,5,5,2,10,2,2,10,1,5,5,1,1,10,5,2,
5,1,3,2,5,10,10,5,10,1,10,3,15,1,10,5,2,3,5,10,3,15,30,5,20,1,2,2,
1,3,7,8,10,5,7,5,9,6,5,8,9,7,6,5,5,7,6,2,3,10,10,15,5,1,2,5,2,1,3,
10,1,5,1,10,1,5,1,2,15,5,1,15,1,5,5,10,15,5,2,10,0,0,5,6,0,1,2,0,3,
0,1,5,7,2,5,1,2,1,10,2,2,2,5,5,10,5,0,5,2,10,1,1,3,10,3,1,4,2,0,1,
5,1,8,5,5,1,3,5,5,2,1,5,5,5,5,0,5,0,13,10,2,9,2,0,0,5,5,5,5,5,0,1,
0,2,1,5,4,2,5,4,1,1,5,1,1,15,10,5,0,15,15,0,0,4,5,2,15,5,15,3,3,
10,10,5,3,7,13,0,0,2,4,1,2,4,1,5,3,8,10,10,5,10,2,5,10,7,10,8,2,5,
7,6,7,5,2,5,1,2,1,8,4,10,5,15,10,5,3,1,5,2,5,1,2,5,1,1,5,2,1,5,0,
10,20,5,5,2,2,10,5,2,0,1,1,2,1,1,1,1,1,1,1,1,2,1,3,1,1,5,2,3,1,2,
0,1,1,5,1,5,2,2,2,5,5,5,15,15,5,10,5,5,15,5,10,5,10,5,7,5,1,5,7,5,
10,1,2,3,2,1,2,1,3,5,3,5,3,2,4,5,2,1,5,5,20,5,10,10,10,10,5,3,5,2,
10,4,1,3,5,5,4,7,5,3,5,2,2,10,4,0,8,2,4,3,15,5,2,8,3,10,5,20,2,0,
0,10,1,1,1,1,1,1,0,0,2,0,10,20,2,10,2,1,3,2,2,5,3,4,1,5,3,1,1,7,2,
4,5,4,5,5,5,10,1,1,3,5,5,0,0,1,1,1,5,0,0,0,0,1,1,2,0,3,0,10,1,2,1,
1,10,0,2,2,5,1,5,3,5,1,3,3,10,0,0,0,5,5,1,2,1,1,2,3,10,10,5,4,1,5,
5,2,3,1,1,5,1,2,25,0,5,5,2,3,1,1,2,1,2,1,5,5,5,5,15,5,5,1,3,2,5,2,
4,2,10,1,7,10,20,5,10,5,1,3,10,2,20,10,15,1,10,1,5,1,3,2,5,6,3,10,
3,15,7,5,10,1,1,1,1,1,1,4,1,10,0,0,0,0,0,2,0,0,2,0,0,0,10,5,2,2,3,
3,4,1,2,2,10,8,1,3,1,4,15,5,1,5,0,2,0,3,2,3,0,1,5,2,1,0,1,3,1,10,0,
3,3,1,1,1,5,1,1,1,1,1,1,3,1,3,2,10,0,10,2,10,1,1,1,1,1,1,1,0,3,0,1,
3,0,1,4,3,5,1,10,5,2,5,10,2,2,3,15,10,10,5,10,5,2,5,5,10,2,1,2,0,5,
5,2,2,2,2,2,10,10,10,3,10,2,1,1,2,3,1,5,2,1,1,3,4,1,2,1,3,2,1,1,2,
1,2,0,1,3,5,1,3,3,2,1,2,3,2,5,3,2,3,1,3,8,1,4,2,2,4,5,11,1,6,2,10,
3,0,0,0,20,10,15,5,15,7,7,10,3,5,2,3,1,0,0,0,0,5,1,3,2,1,1,1,2,1,2,
2,5,2,1,1,2,1,2,0,0,3,0,0,0,2,2,5,5,5,1,60,15,2,0,3,5,5,1,2,10,2,0,
2,15,5,1,20,3,0,10,0,5,10,0,0,10,0,0,5,0,5,2,2,10,1,1,5,1,5,2,5,2,
15,20,15,5,5,5,15,5,2,10,20,1,1,2,1,1,5,1,5,3,3,1,3,15,6,15,10,10,
15,20,10,1,1,1,3,3,4,4,15,1,10,5,5,4,0,1,2,2,2,2,3,2,3,5,2,1,1,2,
3,2,5,15,4,3,1,5,0,1,2,1,3,0,1,5,1,1,0,5,0,0,0,10,5,5,5,5,10,0,1,
1,2,15,10,30,1,1,0,2,3,2,4,5,10,3,10,1,1,1,7,3,1,3,3,3,10,5,3,2,7,
0,5,2,0,30,20,10,10,10,10,10,10,10,10,10,5,5,5,5,10,2,5,5,2,20,5,
30,15,10,5,6,5,20,1,10,10,1,1,5,5,1,5,5,10,15,15,5,10,10,5,3,3,5,
10,5,0,5,5,1,5,5,15,20,5,5,5,1,15,5,20,1,2,10,1,2,0,1,5,5,10,1,5,
1,1,1,1,1,2,2,10,10,3,5,0,3,1,1,1,0,1,3,1,1,5,0,10,5,0,0,3,3,5,0,
0,1,10,5,5,3,10,10,10,2,35,20,25,15,5,5,2,2,5,2,5,0,3,3,1,30,10,
15,5,20,5,10,10,20,15,5,10,5,5,15,20,15,5,0,1,4,10,3,4,26,5,10,10,
1,5,0,0,5,5,5,5,10,30,2,2,5,1,3,3,1,1,1,3,1,3,7,3,15,20,0,15,5,25,
3,25,0,30,0,5,1,1,2,1,1,5,10,5,0,0,20,1,0,15,5,5,15,15,15,15,15,10,
10,15,10,30,30,20,20,5,5,1,4,4,5,5,10,2,0,5,1,1,15,15,5,4,1,1,3,3,
1,0,15,0,10,20,15,5,4,0,0,2,1,0,2,0,2,1,1,2,2,1,0,5,4,3,3,5,5,2,1,
5,4,2,10,2,2,10,3,3,5,10,1,0,10,5,0,10,5,10,5,10,10,60,30,30,99,0,
2,1,0,1,1,2,1,2,1,5,1,1,1,5,5,5,1,0,1,0,0,0,0,3,3,10,2,5,2,2,1,5,3,
6,2,3,7,5,3,1,1,1,1,1,5,5,5,5,7,2,5,5,10,2,2,5,5,5,10,5,5,5,5,5,5,
10,15,5,5,5,5,0,2,10,0,2,5,0,1,10,2,1,1,2,4,5,1,2,2,0,5,2,2,3,3,1,
1,10,0,3,0,1,10,12,3,2,6,9,3,5,2,1,1,1,3,4,5,10,5,10,15,20,6,5,5,
5,1,5,15,5,5,10,8,3,15,12,0,5,2,5,5,3,5,4,1,1,3,1,5,2,10,20,1,15,
15,10,3,1,3,2,0,5,0,1,0,1,2,2,1,1,0,1,10,1,5,1,1,1,4,0,5,1,1,15,10,
1,5,5,5,1,10,0,10,2,1,99,99,99,99,99,5,1,10,30,3,5,5,10,10,0,10,0,
4,1,12,5,1,4,1,3,0,15,3,10,5,1,2,1,1,1,2,1,0,1,1,3,5,2,25,15,20,1,
5,2,10,3,3,4,1,3,2,1,5,3,10,1,10,5,1,25,5,20,10,20,15,15,10,10,18,
0,5,1,0,5,2,10,5,5,2,5,5,3,1,3,2,0,2,1,5,99,99,99,99,99,99,99,99,
99,99,2,5,1,3,5,5,0,2,5,7,10,2,15,3,30,20,2,1,0,1,0,1,2,5,4,1,1,1,
2,2,0,2,2,2,2,2,1,3,10,20,15,10,2,3,5,10,5,0,10,10,10,15,1,1,9,2,
1,7,5,5,5,3,2,2,1,2,1,1,5,1,20,2,5,15,5,5,3,5,2,3,15,1,5,3,5,0,5,5,
10,5,7,1,1,1,3,20,1,3,0,5,1,1,1,15,30,5,35,15,5,5,5,2,2,1,1,15,1,
4,3,2,3,1,5,3,1,3,3,2,10,1,5,1,5,1,2,7,30,20,15,5,30,10,10,5,10,10,
10,5,5,0,5,10,10,10,10,10,5,15,10,15,15,15,10,15,20,15,20,20,5,5,
20,10,10,5,1,0,2,5,2,5,5,1,2,2,2,10,1,2,7,2,15,15,15,5,15,5,10,1,
20,2,1,99,0,2,0,5,2,5,1,10,5,5,5,1,5,2,2,5,5,5,3,5,1,0,5,15,7,2,4,
5,5,10,2,10,10,10,3,3,10,5,5,15,5,10,10,2,5,20,5,5,1,5,10,15,1,3,
2,1,3,1,1,1,1,1,1,1,2,1,1,1,1,2,1,1,1,2,2,1,1,1,1,1,3,3,1,5,7,10,
2,5,10,15,2,5,2,2,3,4,3,2,5,4,10,5,3,2,2,2,5,1,1,5,2,5,5,10,5,15,
1,1,1,1,15,2,5,2,10,3,5,2,1,6,5,1,5,5,1,3,5,3,1,4,5,3,5,4,1,8,5,1,
5,5,9,5,5,9,4,3,4,2,5,2,1,5,10,10,5,1,10,1,5,1,1,3,2,1,5,3,3,5,1,
5,1,2,2,0,7,7,2,0,1,3,10,1,2,1,1,5,5,1,5,1,1,2,0,5,15,5,15,5,5,15,
2,2,1,1,10,1,5,10,1,1,1,1,15,1,4,1,1,1,2,1,10,1,5,15,5,10,15,3,1,
1,1,0,5,5,5,0,5,7,1,7,9,2,1,6,5,10,2,2,5,2,8,1,1,1,1,2,5,10,1,10,
1,7,5,4,5,5,5,10,10,15,5,0,10,15,99,99,99,99,5,1,1,2,5,1,5,1,5,5,
10,10,5,10,5,5,10,2,15,0,1,0,7,5,0,1,0,0,5,5,5,3,10,5,3,1,10,15,3,
6,6,1,3,2,0,15,2,20,10,0,1,0,2,5,15,5,2,1,1,5,5,1,5,1,20,15,15,1,
1,2,1,3,0,5,3,0,0,5,6,3,5,6,4,1,2,4,1,10,5,6,3,7,10,5,10,10,5,2,5,
1,1,5,1,2,5,2,5,2,2,2,5,1,8,1,1,1,1,1,4,7,0,3,3,1,3,2,1,6,1,0,2,1,
0,5,1,1,6,1,5,1,3,3,3,3,7,2,10,4,3,5,5,7,3,5,3,6,1,5,1,4,4,3,2,1,
1,2,1,2,15,18,5,0,1,5,0,3,5,0,0,0,1,1,1,3,0,0,1,2,0,2,20,2,4,2,2,
34,0,1,0,4,10,0,7)
thesisdata <- data.table(id = seq(1:length(parktimes)),
parktime = parktimes)
Anscombe <- function(x) {
# https://github.com/broxtronix/pymultiscale/blob/master/pymultiscale/anscombe.py
# Compute the Anscombe variance stabilizing transform.
# the input x is noisy Poisson-distributed data
# the output fx has variance approximately equal to 1.
# Reference: Anscombe, F. J. (1948), "The transformation of Poisson,
# binomial and negative-binomial data", Biometrika 35 (3-4): 246-254
return (2.0 * sqrt(x + 3.0 / 8.0))
}
CalculatePoissonDist <- function(thesisdata, colnam) {
# According to:
# https://www.sqlservercentral.com/articles/scoring-outliers-in-non-normal-data-with-r
# We're going to use the ppois() function to calculate an "outlier score" for
# every observation in our dataset. The intuitive way to think about this
# score is the "likelihood of observing a point this large". This is a
# somewhat loose interpretation of a p-value, but suitable for detecting
# outliers.
# This function fails if input dataframe is not a data.table dataframe.
# Calculate Poisson distribution for parktime or walktime. Creates two new
# columns, Score (double) and Outlier (boolean). Explicitly prints results
# and returns the inputted dataframe with updates.
# Try Anscombe transform for the parameter column
anscombe_col <- paste0("anscombe_", colnam)
thesisdata[, (anscombe_col) := Anscombe(thesisdata[, get(colnam)])]
# Calculate a "p-value" for outliers, based on the poisson probabilities.
# Use get() to enable string column names in data.table syntax
thesisdata[, Score := 1 - ppois(q = get(anscombe_col),
lambda = mean(get(anscombe_col)))]
# Apply a Bonferroni correction factor to the p-value, to control the long-run
# error rate
thesisdata[, Outlier := Score < 0.05 / 1000]
# Add a Method column with all values "Poisson"
thesisdata[, Method := "Poisson"]
# Visualise the results
p <- ggplot(thesisdata, aes(x = id, y = !!sym(colnam))) +
geom_point(aes(colour = Outlier), size = 3, alpha = 0.7) +
scale_colour_manual(values = c("darkgrey", "red")) +
scale_y_continuous(breaks = scales::pretty_breaks(n = 10)) +
theme_minimal()
print(p)
return(thesisdata)
}
# Outliers in count data?
thesisdata <- CalculatePoissonDist(thesisdata, "parktime")
Best Answer
The data is just too noisy to analyze by direct inspection, so I can understand why the question of outliers arose. However, identifying outliers requires that a logical, physical reason for that status to be at least postulated. The only outliers herein are the 99+ answers, which literally lie outside the range of the data. What is occurring with the human responses can be seen using a more precise histogram.
As seen in the minute by minute histogram the responses to your question as to how long it takes to park are responded to with human time estimates that increases at certain time intervals, 1, 5, 10, 15, 20, 25, 30... min. Which are clock face interval estimates. That is we are postulating is that it would be more frequent to say (approximately) 15 min rather than 14 or 16 min. Consequently, it is hard to find a distribution that fits the data as raw data. However, I did a Gaussian kernel smooth on the data (in Mathematica) just to get some idea of what it looks like and got.
Now, without smoothing I got
About that, if one ignores the mixture distributions, which are attempting to model the noise, and not very successfully, one is left with a geometric distribution or a negative binomial distribution.
After smoothing the candidates are a gamma distribution or a beta distribution. I noticed that in the raw data the maximum value of 99 is populated several times, which is likely why the beta distribution was identified after smoothing.
Thinking physically about this problem, there are no whole number wait times. That is, no one parks at 1 min elapsed time exactly and actual times might be more closely 5341 milliseconds or 3 min 34.453 s. So a gamma distribution wait time model might be more appropriate. This is related to a Poisson process, and is a continuous model for it. I would suggest you fit a gamma CDF to the observed CDF, as that will damp the noise without falsifying the model.
To create the CDF, truncate the 99+ entries so that the CDF data for fitting stops at 0.994064, which is $1-\dfrac{31}{5222}$, where 31 is the number of 99+ answers, and 5222 the total number of realizations.
So, just for fun, I did that. The CDF gamma distribution is:
$$\begin{array}{cc} \Bigg\{ & \begin{array}{cc} Q\left(a,0,\frac{x}{b}\right) & x>0 \\ 0 & \text{Elsewhere} \\ \end{array} \\ \end{array}\text{ },$$
where $Q(\cdot,\cdot,\cdot)$ is the generalized regularized incomplete gamma function, and careful as Mathematica might parametrize as b or 1/b compared to other implementations. The coefficients I got from ordinary least squares regression were $a=0.6618887062, b=6.679277804$ and the fit plot was this:
I note that it works a bit more realistically if I shift the data one min to the right. In that case $a=1.113789864, b=4.648996063$. Then as $a>1$, the pdf gamma distribution assigns 0 probability of parking in 0 time (which is physical because human reaction time is not zero, it can be within the first minute, which is <1 but not zero. Same confusion for birthdays, the first birthday is when the first year has finished.) and has a peak at 0.529008630 min, as below
Which has the following density formula:
$$\frac{b^{-a} t^{a-1} e^{-\frac{t}{b}}}{\Gamma (a)},$$ where $t$ is time in minutes, and where $a=1.11379, b=4.64900$-min, and $a$ has no units (dimensionless). That is, $$0.190915 e^{-0.215100 t} t^{0.113790}.$$
BTW, the median wait estimate is 3-min from the raw data.