1 We address the problem of detecting highly raised crowd density in situations 2 such as indoor dance events. We propose a new method for estimating crowd den- 3 sity by anonymous, non-participatory, indoor Wi-Fi localization of smart phones. 4 Using a probabilistic model inspired by statistical mechanics, and relying only on 5 big data analytics, we tackle three challenges: (1) the ambiguity of Wi-Fi based in- 6 door positioning, which appears regardless of whether the latter is performed with 7 machine learning or with optimization, (2) the MAC address randomization when 8 a device is not connected, and (3) the volatility of packet interarrival times. 9 The main result is that our estimation becomes more { rather than less { accu- 10 rate when the crowd size increases. This property is crucial for detecting dangerous 11 crowd density. 12 Keywords: big data analytics, crowd density estimation, probabilistic modeling, 13 indoor Wi-Fi localization 14 1 Introduction 15 Crowd disasters have taken many human lives. The Love Parade disaster in Duisburg, 16 2010, the Ellis Park Stadium disaster in Johannesburg, 2001, the PhilSports Stadium 17 1 stampede in Manila, 2006, are just a few examples. One of the major factors contributing 18 to crowd disasters are critically dense spots [1, 2, 3], which are dicult to detect due to 19 lack of macroscopic overview of the crowd [1]. In this paper we address the problem of 20 estimating the crowd density distribution in situations such as indoor dance events, to 21 enable prevention of crowd disasters. 22 A lot of research on estimating crowd density concerns processing video records from 23 security cameras [4, 5]. However, this approach does not suce to detect critically raised 24 crowd density. Firstly, as mentioned before, it is dicult to obtain macroscopic overview 25 of the crowd. Secondly, the lighting conditions at a concert might not be sucient 26 for video-based crowd analysis. Finally, the error of counting people increases with the 27 increase of the actual crowd density [6] due to the so-called occlusion e ects. Another way 28 to monitor the crowd density is by using RFID technology [3]. Each participant is asked to 29 wear a tag, and RFID readers are distributed across the venue. This approach, however, 30 requires participation from the crowd and deployment costs. Similar requirements exist 31 for other wireless tracking technologies, like Bluetooth or GPS-based. In addition, GPS 32 is not so suitable for indoor localization. (For more information on crowd monitoring 33 services, we refer the reader to [3].) 34 In our approach, which can complement the video-based analysis, we exploit the 35 ubiquity of smart phones, as it has been done in [7, 8, 9, 10, 11, 12, 13]. More concretely, 36 our approach is non-participatory, that is, it does not require participation from the 37 crowd, and uses the already existing Wi-Fi network at the venue. 38 Despite the recent success in using wireless technologies for indoor positioning and 39 crowd counting, several problems remain open. Firstly, there is the problem of ambiguity 40 when attempting wireless indoor localization [14, 15, 16, 17], which is a major source of 41 localization errors. Secondly, when a phone is not 'connected', its Media Access Control 42 (MAC) address may change (be \randomized") over time [18], complying to privacy 43 policies, thus making it impossible to track the user over time. Finally, since we do not 44 rely on crowd participation, the signals from the phones are quite irregular in time [19], 45 meaning that real-time tracking of a device is also challenging. 46 In this paper we address the aforementioned three problems as follows. To address 47 the ambiguity problem, we apply concepts from statistical mechanics: rather than esti- 48 mating the most likely position of a visitor in real time, we create an evolving probability 49 distribution over all possible positions of the visitor. We rely on the fact that we have 50 a lot of data (or visitors), to estimate the crowd density by aggregating the individual 51 distributions. We use the abundance of data again and the fact that the structure of 52 the MAC address reveals whether it has been randomized, to account for the fact that 53 a portion of the devices are not trackable. Finally, we deploy a time-out based memory 54 model for dealing with the volatile signal rates. 55 Applying primarily the law of large numbers leads to our main result; namely, that 56 our estimation becomes more (rather than less) precise when the crowd size increases, 57 even without requiring crowd participation. This property is crucial for being able to 58 detect critically raised crowd density. 59 The rest of the paper is organized as follows. Section 2 explains brie y our data 60 collection process and the Wi-Fi localization methods that we use. Section 3 introduces 61 in more detail the problems related to crowd density estimation. Section 4 proposes a new 62 method for crowd density estimation that addresses the mentioned problems. In Section 5 63 2 the performances of our method are analyzed and discussed, including comparing the 64 method to related work. Section 6 ends with conclusions and directions for future work. 65 2 Background 66 2.1 Data collection and privacy protection 67 The in-house data and videos used in this paper were collected during the Sensation 2015 68 dance event in the Amsterdam ArenA (today Johan Cruij ArenA) football stadium. 69 More than 30000 visitors were present, and 28847 MAC addresses were detected in the 70 range of the Wi-Fi access points (AP's). We used 30 AP's distributed in the east corner 71 and in the west side of the stadium. (The white dress code of this particular dance event 72 in 2015 made it suitable to evaluate with video data.) We processed the Wi-Fi signals 73 and estimated the coordinates of the Wi-Fi enabled devices using a method similar to 74 trilateration, that we explain in the next subsection. 75 Usage of smart phones to identify the user's locations inevitably raises the question 76 of privacy concerns. The system that we use has been designed from the ground up with 77 privacy in mind - no privacy-sensitive data is ever stored. Only a minimal set of data 78 is collected (timestamp, access point, signal strength and identi er). The unique MAC 79 identi ers of the phones are hashed (anonymized) on reception. The data is not stored on 80 site, but passed on streaming to a trusted third party. The third party maps the hashed 81 identi ers once again. The nal identi er is stored in an environment accessible only to 82 the data scientists participating in this project. As a result, none of the involved parties 83 has sucient information to recover the original MAC address. 84 Following the European Union General Data Protection Regulation (GDPR) and the 85 Dutch law for handling personal data, laid down in the Personal Data Protection Act 86 (Wbp) (Section 12.2) [20], we do not publish any part of the data, and only reveal 87 statistical and aggregated results about the crowd. The data analytics Jupyter Notebook 88 scripts together with the output (aggregated results) that led to the insights and research 89 results presented in this paper can be found in [21]. 90 2.2 Localization of smart phones using Wi-Fi sensors 91 Smart phones transmit Wi-Fi signals which are captured at the Wi-Fi access points. 92 The captured signals contain information about the measured received signal strength 93 (RSS). Widely used methods for positioning using RSS values are multilateration and 94 ngerprinting [22]. In what follows we give a brief overview of the two methods. 95 2.2.1 Localization by multilateration 96 Using the Friis equation for the relationship between an RSS and the distance between 97 a transmitter and a receiver, the distance between a smart phone and an AP can be 98 estimated. When we have the exact distances from a smart phone to at least three AP's, 99 the position of the smart phone can be uniquely determined at the intersection of three 100 circles (Figure 1a). [23]. 101 3 Figure 1: Estimating position with trilateration: a) precise b) rough However, in practice the measured RSS values contain unpredictable variation due to 102 noise and interference such as absorption and re ection by obstacles (e.g. human bodies). 103 As a result, the circles do not intersect at a unique point (see e.g. Figure 1b). In this case, 104 an optimization procedure is undertaken for positioning the smart phone using more than 105 three APs and, in our case, all received RSS values within 500 milliseconds. We use the 106 least-squares optimization method, which in our case has the form of a chi-square data 107 t. An exact description of the method is beyond the scope of this paper, and we refer 108 the reader to [24, 25] for details. We note, however, that the statistical estimation of the 109 position provides us also with the standard deviation, which we will use in Section 4. 110 2.2.2 Localization by ngerprinting 111 Traditional ngerprinting is also RSS-based. There are 'oine' and 'online' phase of the 112 localization. In the 'oine' phase, for a moving client device, the signal strengths from 113 several access points in range are continuously recorded and stored in a database, along 114 with the known coordinates of the device [26, 27]. During the 'online' tracking phase, the 115 current RSS vector of a device at an unknown location is compared to the vectors stored 116 in the ngerprint, and the closest match is returned as an estimation of the location. 117 The closest match is usually determined through probabilistic methods (e.g. expectation 118 maximization, KL-divergence), or through machine learning techniques (e.g. k-nearest 119 neighbors, Support Vector Machines, neural networks). 120 3 Problem statement 121 Under ideal circumstances, the positioning itself would suce to estimate the spatial 122 crowd density distribution: every second we would only need to count the number of 123 detected devices per square meter. However, Wi-Fi based localization comes with the 124 following challenges, which are not related to the mathematical methodology behind the 125 'positioning' step, and that prevent us to apply direct counting. 126 1.Issue 1: Ambiguity of the localization procedures . It has been argued that one of 127 the biggest source of errors, when using RSS values for localization, is the absence 128 of a single global optimum [15]. In the ngerprinting approaches the phenomenon 129 is called " ngerprinting twins" [14, 15], while in the multilateration approaches it 130 is known as " ip ambiguity" [16, 28, 17]. We also sampled randomly 20 MAC 131 addresses from the Sensation data, under various crowd conditions, and plotted 132 4 the estimations of their coordinates through time. We plotted only the estimations 133 with a relatively small (conditional) uncertainty. We observed persistent bi-modal 134 distributions of the estimations through time (as if the device is being tele-ported 135 constantly), an example of which can be seen in Figure 2 (a). This gure shows the 136 estimated x-coordinates (in meters) through time of a static MAC device that was 137 persistent for 24 hours (most probably an AP). 138 (a) (b) Figure 2: (a) Estimation of the x-coordinate in meters of a static MAC device through time; (b) RSS-based positioning can lead to multiple local optima The "twins" phenomenon originates from the fact that the environmental settings 139 and the temporal variations in the RSS create opportunities for multiple local op- 140 tima [17], in case of trilateration, or geographically distant positions to share the 141 same RSS vectors [15], in case of ngerprinting. (Note: we call them 'twins' be- 142 cause the empirical evidence suggests so far that there are no 'triplets'; however, 143 theoretically the latter are not excluded.) To understand the problem, consider 144 Figure 2 (b). In the center of every ring there is an access point with an estimated 145 signal strength to a particular MAC device, with a certain error range. The error 146 range is represented by the width of the ring. Then, there are two possible regions 147 where all three rings overlap (A and B), and that are equally good candidates for 148 positioning the MAC device. Note that this problem can arise regardless of the 149 width of the rings and the number of APs, and that it cannot be alleviated without 150 using additional information or assumption about the visitors or the environment. 151 For example, current solutions require crowd participation [15], or assume crowd 152 mobility [14]. 153 2.Issue 2: Volatility of packet rates . When a MAC device is connected to a Wi-Fi 154 network, it sends packets with a relatively stable and frequent rate. However, when 155 the device is not connected, it is in a \probing" mode, i.e. searching for a network, 156 and in this case the packet rate is quite volatile, ranging from a few seconds to a few 157 minutes [19]. This introduces challenges when trying to estimate the total number 158 of devices present at a certain moment. 159 3.Issue 3: MAC address randomization . Due to the ever-increasing privacy concerns 160 and possibly other business reasons, starting from 2014, Apple has introduced 'ran- 161 5 domization' of the MAC addresses of the phones, when the latter are in probing 162 mode (not connected to the internet) [19]. This means that the devices not con- 163 nected to the internet are continuously changing their MAC addresses and cannot 164 be followed over time. (It is worth noting, however, that from the format of a MAC 165 address itself it can be determined whether the address is authentic or randomized). 166 4 Method 167 In this section, as a main contribution of the paper, we propose solutions to the issues 168 stated in the previous section. 169 4.1 Estimation under localization ambiguity 170 Creating statistical ensembles In order to deal with the localization ambiguity, let 171 us start with the following observation: we are notinterested in the individual locations 172 of the MAC devices, but rather in the density of the crowd. In addition, for prevention 173 of crowd disasters, it is very important to have a precise estimation when the number of 174 people in a stadium is large and dense spots are likely [1, 2]; at the same time, obtaining 175 precision at a low crowd density is of less priority. 176 To this end, we propose a probabilistic model for crowd density estimation. To explain 177 our idea, we use the scenario depicted in Figure 2 (b). We can say that the particular MAC 178 device is located in region A with a probability of 0.5 and in region B with a probability of 179 0.5 (we assign the probabilities in a trivial way for the purpose of illustration). Although 180 this approach does not provide us with very useful information about the location of the 181 MAC device, if we apply the same reasoning for all MAC devices, and we add together the 182 spatial probability distributions of all MAC devices, we end up with a spatial distribution 183 of the crowd density. If we assume that the locations of all MAC devices are mutually 184 independent and identically distributed, we can apply directly the standard law of large 185 numbers and conclude that, for a large crowd, the error of estimation of the density 186 per square meter will vanish. (The di usion animation in [29] provides a nice visual 187 demonstration of the concept.) However, we cannot make those assumptions, because 188 e.g. people tend to go to concerts in groups, that is, their locations are correlated [30, 31]. 189 In this case, the variance of the estimation in the limiting case is equal to the average 190 covariance between the locations. In Section 5 we will show that the average covariance 191 still tends to zero as the crowd size increases, because for a regular concert crowd, the 192 group size is relatively small compared to the entire crowd. 193 Computing individual probability distributions Next, for an arbitrary MAC de- 194 vicemand region R, we proceed with de ning Prob (m2R), that is, the probability 195 thatmis located in R, in order to be able to evaluate the estimated crowd density per 196 region. 197 The localization (" tting") provides us with a series of estimated positions for a mobile 198 device. We estimate the spatial probability distribution for the device along a moving 199 time window. Our rst step at time tis to select the Nestimated positions whose time 200 stamps fall within a speci ed time window [ tt;t] from the data. A natural way to 201 6 construct a two-dimensional probability distribution is to create a histogram, by binning 202 theNpositions and normalizing by N. In case the positions have been estimated with 203 multi-lateration, the optimization procedure provides also for a Gaussian error ( xi,yi) 204 of any estimate ( xi;yi) (see Section 2.2). Therefore, we rst \smooth" each of the ( xi;yi) 205 positions into a bivariate distribution, using a Gaussian kernel with standard deviation 206 (xi,yi). Then, for a MAC device mwe generate a two-dimensional probability density 207 function (pdf) by adding up the separate `bumps' and normalizing by N. (See Figure 3 208 for an example of smoothing a histogram.) 209 Figure 3: Smoothing a histogram with Gaussian kernels. (a) Original histogram. (b) Smoothed histogram. Formally, the implementation of our method is similar to that of kernel density esti- 210 mation [32, 33]. In our case the amount of smoothing is determined by the uncertainty 211 valuesxandy. The pdf for a MAC device mat a location ( x;y) is de ned by 212 ^fm(x;y) =1 NNX i=1K((xxi);xi)K((yyi);yi) (1) where the kernel function Kis given by 213 K(u;) =1 p 2exp u2 22 (2) and (xi;yi) is the result of positioning at step i2f1;Ng. 214 In order to evaluate Prob (m2R), we need to integrate ^fm(x;y) for (x;y)2R. The 215 nal crowd density estimation at point ( x;y) is given by 216 ^fT(x;y) =X m^fm(x;y): (3) Finally, in order to estimate the number of people in region R, we need to integrate 217 ^fT(x;y) over the region R. 218 Remark 1. In our implementation we scale up the individual probability distribution, 219 such that it integrates to one inside the region of the stadium (the concert venue). We 220 assume that if a device is detected by the AP's, it is inside the stadium, and thus the 221 probability that it is in the stadium should be one. (It is very unlikely that the AP's 222 7 have detected devices that are outside, due to the thick walls of the stadium.) In the 223 future we also plan to include the map of the stadium in the calculations, to incorporate 224 the fact that the probability that a visitor is in an inaccessible region is zero. 225 Note that so far we assumed that in every time window there is at least one estimate 226 for every MAC device present in the stadium. In what follows we explain how we capture 227 the cases when this assumption does not hold. 228 4.2 "Conservation of mass" under packet rates volatility 229 To address the second issue raised in Section 3, Volatility of packet rates , we ensure that 230 we do not forget about the MAC devices that were not observed in the last time window. 231 In fact, for every MAC device that was ever observed, until it is observed again, we 232 maintain the old probability distribution. However, we also apply a time-out, that is, if 233 a MAC device has not been observed in a long enough time interval (called 'memory' 234 parameter), it is simply removed from the pool of MAC devices. 235 4.3 Estimation under MAC address randomization 236 The previous discussion assumes that a MAC device does not change its identi er over 237 time. However, as we noted in the last issue in Section 3, a device in a probing mode 238 might randomize its address during a time window, leading to it being counted twice. To 239 address this problem, we rely once again on the fact that we have a lot of data, and on 240 the fact that we can derive from the structure of the MAC address whether it has been 241 randomized or not. Figure 4 (a) shows the time series of numbers of non-randomized and 242 randomized addresses observed per minute from midnight until around 6:00 am during the 243 Sensation concert. We observe that their ratio is stable through time (Figure 4 (b)); the 244 Pearson correlation coecient between the time series of randomized and non-randomized 245 addresses is 0.88. 246 (a) (b) Figure 4: (a) Number of non-randomized and randomized addresses observed per minute; (b) Linear regression plot of randomized and non-randomized addresses observed per minute Therefore, when estimating crowd density, we ignore the MAC devices that have random- 247 ized addresses and at the end we multiply the crowd density by a factor to account for 248 8 the discarded MAC devices. This factor is derived from the slope of the linear regression 249 t of the two time series (Figure 4(b)), which in our case turns out to be 0.2, with a 250 standard error of 0.006. 251 Note also that this proportion should be re-computed periodically, to account for the 252 changing conditions at the smart phones market. In fact, when the crowd is large like in 253 our Sensation scenario, the randomization factor can be updated in real time, during the 254 concert hours, by using all data that arrived in e.g. the last hour. 255 We envision, however, that in the future more people will be connected to the Wi-Fi 256 (and thus the proportion of randomized addresses will become smaller) due to the fact 257 that an increasing number of stadiums across the world o er \smart" services, but also 258 due to the increasing usage of social media to post photos and videos of an event in 259 real-time. 260 Some studies [34, 19] suggest that it is still possible to follow devices despite MAC 261 randomization, and it would be interesting in the future to see if we can improve our 262 methodology taking those studies into account. 263 5 Results and discussion 264 In this section we analyze our method for crowd density estimation in various manners: 265 analytically, with simulations, and using two real-life datasets. (Note: Section 5.1 o ers 266 a formal validation that readers uninterested in technical details can skip without loss.) 267 5.1 Theoretical analysis under correlated groups 268 We show formally that the relative error of the crowd density estimation converges to zero 269 when the crowd size increases, despite having correlated groups of visitors (e.g. friends). 270 Letfmac 1;mac 2;:::mac ngbe allnMAC devices detected in the stadium at time t. 271 Due to the results in Section 4.3, where we show how we can safely discard the randomized 272 addresses from the analysis, we can assume here that all nMAC addresses are xed. In 273 what follows we omit tfrom the notation for clarity. Let Rbe an arbitrary region of the 274 stadium. Denote by mac i2Rthe statement \the device mac iis in region R". LetXi 275 be a random variable de ned by 276 Xi=1 ifmac i2R 0 ifmac i=2R(4) Denote byXthe total number of devices in Rdetected at time t. Clearly,X=Pn i=1Xi. 277 ThenE(X), the expected value of Xis 278 E(X) =E nX i=1Xi! =nX i=1E(Xi) =nX i=1(1Prob (mac i2R) + 0Prob (mac i=2R)) =nX i=1Prob (mac i2R);(5) where byProb (mac i2R) we denote the probability that mac iis in the region R at 279 timet. We will show that the variance of X=n, that is, the variance of the proportion of 280 9 devices detected in region Rout of all detected devices, diminishes when nbecomes large 281 (note that the variance of Xin the limiting case is out of our interest because in this case 282 E(X) is also potentially in nite). This suces to show that our method for estimation 283 of crowd density is theoretically sound, given the probabilities Prob (mac i2R). We have 284 VarX n =1 n2Var nX i=1Xi! =1 n2 nX i=1Var(Xi) +X i6=jCov(Xi;Xj)! (6) Let be an upper limit on the number of people going to a concert together (i.e. 285 whose locations are correlated). Note that, because the random variables fXign i=1take 286 values inf0;1g, the covariances Cov(Xi;Xj) take values in [1;1]. Thus, the covariances 287 are upper-bounded (by 1). Denote by 1 the maximal covariance between any Xiand 288 Xjand let us write ijif and only if the owners of the MAC devices mac iandmac j 289 are in the same group of friends. Then, 290 X i6=jCov(Xi;Xj) = 2X 1i<jnCov(Xi;Xj) = 2X i;j:ijCov(Xi;Xj) 2n ( 1) 2=n ( 1):(7) Here the inequality holds because the maximal number of groups is nand the maximal 291 number of pairs ( i;j) in a group is ( 1)=2. Let us denote by =1 nPn i=1Var(Xi) 292 the average variance of X1;X 2;:::;X n(note that1 from the de nition of fXign i=1). 293 From (6) and (7) we have 294 VarX n 1 n2(n+n ( 1)) =1 n(+ ( 1))1 n(1 + ( 1));(8) which tends to 0 when n!1 . Note that we have greatly overestimated the covariance 295 with the inequality in (7), which means that in practice the variance converges to 0 much 296 faster than presented. 297 With the above, we have proven a version of the law of the large numbers that is 298 generally applicable. We formalize our results in the following proposition: 299 Proposition. LetfX1;X 2;:::X ngbe random variables that always take values in a bounded 300 real interval. Suppose that the set fX1;X 2;:::X ngcan be partitioned into subsets of max- 301 imal size (a xed constant independent of n), such that if XiandXjbelong to di erent 302 subsets, then Cov(Xi;Xj) = 0 . LetSn=1 nPn i=1Xi. Then 303 lim n!1Var(Sn) = 0: Remark 2. In our proof we have assumed that correlation happens only within groups 304 of friends. Note that this is a sucient, however not a necessary condition for the conver- 305 gence of the variance. A necessary condition is that the average correlation in the crowd 306 tends to zero as the crowd size increases. This allows for (positive) correlations outside 307 10 groups of friends. Moreover, in a crowded situation negative correlation is more likely to 308 happen, where people move away from the crowd, looking for empty spots [35], which 309 reduces the variance further. In respect of this discussion, however, it is worth noting 310 that there is onesingularity scenario, a \crowd crush", when the crowd is so dense ( 311 >6 persons per m2) [2] that people cannot move freely anymore and the entire crowd 312 becomes a \group" (as in the Love Parade disaster), implying that all visitors locations 313 are (positively) correlated. In this work we aim to detect high density with our method 314 way before this saturation happens, to be able to react preventively; otherwise, it is too 315 late. 316 Remark 3. In our proof we assumed that all MAC addresses are xed, i.e. that there 317 are no randomized addresses. However, as discussed in Section 4.3, in reality we omit the 318 randomized addresses from the analysis and at the end we multiply the estimation by 319 a so-called randomization factor. Note that, again due to the law of the large numbers, 320 the larger the crowd, the more precise is the estimation of this factor, when it is re- 321 estimated in real time. This means that the above convergence would not be a ected 322 by the randomization factor. However, in the following subsection we also con rm this 323 experimentally. 324 5.2 Analysis with simulations 325 Section 5.1 gives a theoretical validation of the method. We proceed with analyzing 326 the method experimentally, i.e. quantitatively. We are especially interested in how our 327 method performs at dangerously high densities, i.e. more than 4 persons per square 328 meter. Such densities are however dicult to obtain from real life scenarios; moreover, 329 if we use video data as ground truth, the latter is inaccurate at high densities [6]. On 330 the other hand, controlled experiments with thousands of participants and high induced 331 crowd density are beyond the scope of a research paper because of security risks. To 332 validate the method at high densities, we use data-driven stochastic simulations [36]. 333 5.2.1 Simulations setup 334 For linearly increasing crowd size in a football stadium (play eld) we simulate the lo- 335 calization data and apply the proposed method on it. We use the Sensation dataset 336 to derive the probability distributions of the uncertainties of the localization procedure 337 that were discussed in Section 3. Concretely, our analysis of the dataset shows that (i) 338 Median packet inter-arrival time for non-randomized MAC devices is 33 sec at the most 339 overloaded AP, and the time is exponentially distributed (Figure 5); (ii) Mean distance 340 (obtained by random sub-sampling) between the modes of the two Gaussians from Sec- 341 tion 3 for static devices (APs) is circa 20m (see also Figure 2(a)), and the standard 342 deviations of the Gaussians are on average around 3m, regardless of the crowd size; (iii) 343 the average Gaussian error of the positioning process is also 3m and the distribution of 344 all errors is exponential. 345 The crowd simulation is performed as follows. Every person moves in a zig-zag motion, 346 in a random direction, because it has been discussed [35] that under high density visitors 347 try to escape in a sort of zigzag motion, without any preference on direction, but only 348 taking any free space nearby. The velocity is limited by the current crowd density via the 349 11 Figure 5: Distribution of packet inter-arrival times for a non-randomized address. Weidmann's equation [37, 7]. We also introduce correlated positions: for crowd size of n, 350 the number of groups is n=4, and every person is randomly assigned to a group (thus the 351 average group size is 4). 352 After recording the original simulated positions for every second, we create the syn- 353 thetic localization data. The "twins" or "teleportation through time" e ect is introduced 354 by displacing the positions randomly, following the distributions in (ii) above. The po- 355 sitioning errors are sampled following (iii). The packets (and thus location ts) are sub- 356 sampled randomly with inter-arrival time drawn from an exponential distribution with a 357 median 33 s, following (i) above. Each MAC address is randomized with a probability of 358 0.15, following the results in Section 4.3. The randomized addresses change their MAC 359 values in every sent packet. 360 Remark 4. We opted to implement our own crowd simulation instead of using available 361 crowd simulators or models, because the latter are built for di erent purposes. In other 362 words, in crowd simulators it is challenging to implement high crowd density [35] and 363 correlated groups of friends [38]; on the other hand, we are clearly not interested in 364 the exact pattern of movement of the visitors, which is the focus of crowd simulators { 365 any positive e ect from simulating visitors movement on a ne scale is attened by the 366 uncertainties of circa 20m and the signal sparseness introduced by the Wi-Fi data. 367 5.2.2 Deriving the optimal parameters values 368 Our method requires two parameters: the length of the time window, in which the location 369 estimates for a device are counted towards its probability distribution, and the time-out, 370 in terms of number of windows, for keeping the old distributions ('memory'). 371 We argue that the value of the time window  tshould be a compromise between 372 the following requirements: (1) having a window long enough such that we can expect 373 to localize each present non-randomized MAC device at least once in it (which should 374 suce since we also have the memory parameter to make sure that we don't miss devices) 375 and (2) having a window that is not too long, because of the crowd mobility. In order 376 to estimate  t, we combine (1) the median packet inter-arrival time per non-randomized 377 MAC address at the most overloaded access point, which from our data is 33 seconds, 378 and (2) the results from [19], where the author concludes that, while in probing mode, on 379 average smart phones send probing requests 55 times per hour. We choose as intuitively 380 12 good value for  tto be roughly 40 seconds. For high crowd density this value is still small 381 enough to not cause outdated positioning. For example, if the crowd density is 4 p=m2, the 382 maximum distance that a person can travel in 40s according to the Weidmann's equation 383 is 8m[7], which is acceptable given the localization uncertainties. If the crowd density is 384 5p=m2, the maximum distance that a person can travel in 40s is 2 m. 385 Having chosen a window size of 40 seconds, we need to decide how often to update the 386 windows. In practice this would depend on the available computing infrastructure and 387 how often one would want to update the probability distributions. For our experiments 388 we chose to have slightly overlapping windows. The overlap is 10 seconds and thus the 389 'timestep' or 'stride' is in this case 40-10=30 seconds. Thus, probability distributions are 390 updated every 30 seconds. 391 Considering the memory parameter, our experiments show that it depends on the 392 average crowd density, and ranges from 0 windows backwards for expected average density 393 <1p=m2, to 5 windows for densities >4p=m2. Intuitively this is justi ed by the fact 394 that a crowd of low density can move faster and old distributions become soon outdated, 395 while a dense crowd moves slowly and having higher memory enables not missing any 396 visitors, for the reasons of being able to detect highly raised crowd density. On the other 397 hand, it can be checked that the probability that a device will send two signals within 398 5 windows is greater than 99% and it is thus not necessary to keep old estimates for 399 longer. Figure 6 gives examples of the performances of our method w.r.t. the ground 400 truth and tted data for various values of the memory parameter (we will analyze the 401 case memory=5 in Section 5.2.3). The measurements are taken after 9 windows, to allow 402 enough time for the system to 'warm up' [36]. The crowd density estimation in the case 403 of tted locations, where the ts are sparse, has been performed by making a snapshot 404 of all devices that have been detected in the last 40 sec. Note that the latter estimation 405 already incorporates partly our method, regarding the window size; however we present 406 it in the gure as it serves as a reference point. Figure 6: Performances for full simulation, varying the 'memory' parameter Figure 7: Performances for the 'static' crowd scenario 407 5.2.3 Results from simulations analysis 408 To see the advantage of applying our method instead of counting directly the tted 409 positions, we rst perform the following experiment. We simulate a 'static' crowd with 410 13 Figure 8: Full simulation, memory=5 Figure 9: Scatter plot of case in Fig. 8 Figure 10: Performance of our method for very low crowd density Figure 11: An overview of performance through time for one simulation. Figure 12: Performance of our method under high percentage of randomized MAC addresses. regular packet rates and thus regular ts (every second), and no MAC randomization, 411 where by 'static' we mean that nobody moves. In other words, we simulate only the 412 'teleportation' e ect introduced by the bi-modal distributions of localization through 413 time. Our goal is to check the e ect of creating spatial probability distributions in 414 the stadium. We test the performances of our method w.r.t. simple counting of tted 415 positions in a 4mx4m corner of the stadium. Figure 7 shows the performance of our 416 14 method versus simple counting for 10 independent simulations with increasing total crowd 417 size, that ranges between 5000 and 50000. Our method performs well w.r.t. simple 418 counting of last tted locations. This is due to the fact that our method ensures that 419 all detected devices are in the stadium by re-normalizing the corresponding probability 420 distributions. Next, as discussed in Section 5.2.2, we apply our method with window size 421 of 40 sec, time step of 30 sec and memory of 5 windows to independent full simulations 422 as explained in Section 5.2.1 (with correlated crowd movement, packet delays and MAC 423 randomization), increasing linearly the total crowd size in every simulation, from 2000 424 to 64000, which is the capacity of modern stadiums. The results computed in a 4m x 425 4m square in the middle of the play eld are shown in Figures 8 and 9. To contrast the 426 performance for high crowd densities to the performance for low crowd density, we include 427 an example of the latter in Figure 10, which shows that the method can be expected to 428 under-perform for low crowd density. (Note that this estimation is under a worse case 429 scenario, where everybody moves all the time with a normal walking speed). However, 430 we are focused on having precise estimation under high densities. Figure 11 shows the 431 performances through time of our method for a fairly dense crowd in the same 4mx4m 432 region of the play eld. Again, for comparison we include also the plot of the tted data 433 in the last 40s, noting that this calculation partially implements our method. 434 Note that so far we were assuming that the percentage of randomized addresses is 435 15%, as observed from our real data from the Sensation concert. For completeness, we 436 also check the crowd density estimation by our method in the same 4mx4m region when 437 the percentage of randomized addresses is high. The results, after performing again 438 independent full simulations with total crowd size between 2000 and 64000, are shown in 439 Figure 12. We can see that, although the variance of the estimation is higher when the 440 randomization is higher, there is still good agreement between the estimation and the 441 true value. 442 5.3 Analysis under a university campus scenario 443 With simulations we explored worse case scenarios and showed the e ectiveness of the 444 proposed method. In this subsection we explore how the method behaves in normal 445 conditions, with relatively low expected crowd density. We use the publicly available 446 UJIIndoorLoc dataset [39] that has been designed for benchmarking ngerprinting meth- 447 ods for Wi-Fi positioning. The ground truth of the dataset is given by the GPS (Global 448 Positioning System) locations of 25 devices worn by more than 20 participants. The par- 449 ticipants moved across the Jaume university campus (Figure 13). To generate positioning 450 (" tted") data on which to apply our method, we use the software for Wi-Fi positioning 451 based on RSS ngerprinting provided by [40, 41]. More concretely, we apply the anity 452 propagation method for clustering RSS ngerprints and the positioning method provided 453 by [42] to generate the ts. To be able to have more test data and also because we are 454 not interested in obtaining high precision ts, but rather to see the e ect of our method 455 applied to the latter, we use the small test dataset in [39] for training and the bigger 456 train dataset for testing. (This is also in line with the reasoning in the design of the 457 dataset in [41] where the training dataset is smaller than the testing dataset). We obtain 458 tted locations with an average error of 14m. Note that this error is not of Gaussian 459 type (it could be also due to the "twins" e ect), and therefore we do not apply the kernel 460 15 smoothing part of our method to the probability distributions. The memory is 0 in this 461 low density case and the MAC randomization factor does not apply because this was a 462 controlled experiment. Figure 14 shows the results of applying our method to the ' tted' 463 data. We calculated 433 windows with a time step of 30 sec and window 40 sec. The rst 464 window starts at 1804000 seconds after the rst recorded timestamp in the dataset (in 465 this period the number of detected devices was the highest). For counting people based 466 on detected phones sending GPS locations, we also apply a window of 40 sec, that is, we 467 make a snapshot of the last positions of the devices detected in the last 40 sec (because 468 the GPS locations are also not regularly recorded). 469 Figure 13: The location in the Jaume Uni- versity in which measurements were taken (bounded by a red rectangle). Figure 14: Performance of our method w.r.t. ground truth provided by the GPS locations of phones. 5.4 Analysis under a concert crowd scenario 470 To check how our method performs on the Sensation dataset, we used video data posted 471 online by visitors (security cameras are focused on detecting re and are thus very sensi- 472 tive to red light but not good enough for counting people). We have video frames available 473 from the most interesting moment, towards the end of the concert, when the crowd size 474 drops signi cantly in a short period of 15 minutes due to people leaving. From this ma- 475 terial we extracted four time points from which counting could be done (see Figure 15 476 for an overview). 477 Based on the frames selected from the video material, people are counted by manually 478 clicking on their heads, using simple computer scripts to keep track of the number of 479 mouse-clicks. Using multiple spatial reference points the people's locations are projected 480 from the perspective image to the 2D-coordinate system of the Wi-Fi measurements. We 481 compare the manual people counts to our Wi-Fi based estimates on a 122 :5m2regionR 482 that is an intersection of the region covered by Wi-Fi and the region covered by video 483 (the upper triangle in the 15mx15m square of Figure 16). The region Ris located on 484 the football eld. At the center of the eld the DJ stage was positioned. We used time 485 window of 40s and memory 0 as discussed in Section 5.2.2. The time step is 30s. The 486 results of estimation via the method and the video counts through time are given in 487 Figure 17. 488 16 Figure 15: Overview of video frames used for manually counting people (all time points are +2:00 UTC at July 5, 2015). Videos and timestamps are courtesy of Jessey Helslijnen. Figure 16: Visualization of (a) manual people counts, and (b) the WiFi-based crowd density estimate, both within the comparison region (upper triangle), at time point 05:32:04 (+2:00 UTC) (see also Figure 15 ). 5.5 Time complexity 489 An important practical implementation issue for an algorithm is its time complexity with 490 respect to the input size. Note that our algorithm has a linear time complexity w.r.t. the 491 crowd size, because we never consider people in relation to other people. More concretely, 492 let us rst recap the steps needed to compute the crowd density distribution: 1) compute 493 the randomization factor; 2) discard the randomized addresses using the available ags; 494 3) for every remaining MAC device compute the individual probability distribution; 4) 495 aggregate all mprobability distributions and multiply the result by the randomization 496 factor to derive the crowd density distribution. To compute the randomization factor, 497 one needs to keep track of the percentage of randomized addresses detected in e.g. every 498 minute of the last hour, i.e. to update a counter every time a new MAC address is detected 499 in the last minute. Since the number of signals that arrive per minute is linear to the 500 17 Figure 17: Performances of our method w.r.t. video based counting in a 112.5 m2near the main stadium exit, last 15 minutes of the concert. crowd size, this step has linear time complexity. To compute the individual probability 501 distributions, one needs to keep track of the individual ts in the last time window. Note 502 that the computation of the individual distribution does not depend on the other MAC 503 devices, that is, this step is running in linear time as well, which also holds trivially for 504 the last step. 505 5.6 Comparison to previous work 506 In the Introduction we mentioned the bene ts of using smart phones data over video 507 based approaches for density estimation of indoor concert crowd; in fact, it is our opinion 508 that the two approaches complement each other and in the future we plan to integrate 509 both techniques in real time. Thus, in this section we contrast our method to previous 510 work that estimates crowd density using wireless technologies. 511 A number of approaches [43, 44, 45] estimate crowd density based on variations of the 512 received signal strength indicator (RSSI) values of a Wi-Fi network. Under controlled 513 experiments they observe that the more people (obstacles) are present, the greater the 514 RSSI variations. Fadhlullah and Ismail (2016) [43] apply analysis of variance, Yuan et 515 al.(2011) [44] apply k-means clustering to obtain clusters of similar crowd density based 516 on similar RSSI variations, while Yoshida et al. (2016) [45] apply linear and support 517 vector machines regression. Our work also deals with the fact that greater crowd leads 518 to greater RSSI variation (represented by the width of the rings in Figure 2b). However, 519 with our method the relative error of estimation of crowd density decreases as the (indoor) 520 crowd increases, whereas the relative error of estimation of crowd density in the above 521 approaches increases as the crowd increases. Having a decreasing error is essential for 522 being able to detect dangerous crowd density. 523 Other approaches rely on an assumption that pedestrians move from point A to point 524 B in order to estimate crowd density. More concretely, in Wirz et al. (2013) [7] the authors 525 follow a participatory sensing approach in which pedestrians share their GPS locations 526 on a voluntary basis. Since only a fraction of all pedestrians share location information, 527 they infer the crowd density from the walking speed, based on the assumption that the 528 maximal walking speed of pedestrians depends on the crowd density (and thus they 529 assume that the crowd tends to reach a certain destination). A participatory sensing 530 18 approach is also followed by Anzengruber et al. (2013) [8]. They use time series to 531 predict mobility patterns in crowds of spectators, and related to the event agenda over 532 time. Schauer et al. (2014) [9] focus their attention to airport pedestrians, without 533 requiring crowd participation, but exploiting the predetermined direction of movement 534 of passengers. They count unique MAC devices detected with strong signals by two 535 sensors (nodes) at both sides (public and security) of a security check inside a major 536 airport, to estimate pedestrian densities and pedestrian ow. Versichele et al. (2012) 537 [10] use Bluetooth scanners at strategic locations during 10-day festivities, to analyze 538 spatio-temporal dynamics of pedestrians. Their methodology is proximity-based, which 539 is suitable under mobility assumption. Delafontaine et al. (2012) [11] apply sequence 540 alignment methods for the extraction of behavioural patterns within Bluetooth tracking 541 data. Higuchi et al. (2014) [46] use a participation based approach to estimate the 542 number of people and make advantage of the fact that people move in groups to correct 543 the estimations of individual traces. 544 Other approaches use additional devices distributed in the crowd to improve the 545 density estimation. Weppner et al. (2013) [47, 12] estimate crowd densities by distributing 546 volunteers in the crowd, who are carrying smart phones scanning for Bluetooth devices. 547 The authors then use statistics to combine the di erent measurements in space and time. 548 A recent non-participatory approach based on Wi-Fi ngerprinting is given by Tang et al. 549 (2018) [13]. This approach proposes an online training phase ('dynamic ngerprinting') in 550 addition to the oine training phase, which allows for the localization to use the current 551 environmental settings to update the ngerprints database. This requires distributing 552 Wi-Fi devices among the crowd to capture the latest RSS values, which however does 553 not apply to our case of high density crowd in the play eld of a football stadium. 554 Other participation-based approaches include the following. Li et al. (2015) [48] 555 use neural networks to learn the relationship between RSS and the number of people 556 (similarly to using WiFi ngerprinting for estimating individual locations). Martella et al. 557 [49] deploy a participation-based visitors positioning system in a museum. Some visitors 558 wear sensor bracelets, but there are also sensors distributed across the museum. For 559 every visitor a (discrete) probability distribution over all possible locations is maintained, 560 after determining the most likely oor. The distribution is a weighted average of the 561 neighboring sensors distributions, the weights corresponding to the signal strengths. Then 562 the visitor is localized by computing the average of its spatial probability distribution. 563 Note that, unlike the present paper, this approach implicitly assumes uni-modal spatial 564 probability distribution. On another related note, the localization in [49] is more precise 565 when there are more visitors wearing sensors. In our case data is collected from a xed 566 number of sensors; that is, our increase of precision with increased number of people is 567 not related to increased number of sensors. 568 We summarize the added values of our method over previous work for estimating 569 crowd density using wireless technologies as follows. 570 1. Our method does not require participation from the crowd; 571 2. It does not assume that visitors are moving in any particular direction; 572 3. It does not require extra hardware in addition to the available access points; 573 4. It does not assume existence of a global optimum in the localization procedure; 574 19 5. To the best of our knowledge, our method addresses for the rst time the ambiguity 575 problem in Wi-Fi localization in the context of crowd density estimation; 576 6. To the best of our knowledge, we address for the rst time the MAC address ran- 577 domization while estimating crowd size based on Wi-Fi technologies; 578 7. The estimation of crowd density with our approach increases with the increase of the 579 crowd size (without relying on participation from the crowd), which is an essential 580 property for being able to detect critically raised crowd density. To the best of our 581 knowledge, this property has not been addressed or shown before. 582 6 Concluding remarks 583 We proposed a new method for estimating indoor crowd density based on wireless tech- 584 nologies. Our method does not rely on participation or mobility of the crowd but rather 585 uses a big data analytics approach. We addressed three known challenges: 1) the ambigu- 586 ity of localization procedure due to noisy RSS values, 2) the MAC address randomization 587 when a device is in a probing mode, and 3) the irregularity of the packet interarrivial 588 times. We used probabilistic models to address (1) and (2) and a memory-based model 589 to address (3). We showed formally that the error of our estimation tends to zero as 590 the crowd size increases (which is essential for enabling disaster prevention), even in case 591 when the locations of visitors are correlated as in groups of friends. We used data-driven 592 stochastic simulations to evaluate quantitatively the e ectiveness of our method for highly 593 dense crowds, and we used two datasets to evaluate the performances on real scenarios. 594 Finally, we positioned our work in the context of previous related work and showed the 595 added values of our approach. To the best of our knowledge, the approach presented 596 here is the rst to address the ambiguity and MAC randomization problems in context 597 of no-participatory and no-mobility assumption, and also the rst that has e ectiveness 598 at high crowd density as one of its design principles and results. 599 While we propose solutions to several issues related to estimating crowd density based 600 on wireless technologies, it is important to emphasize that methods involving human 601 behaviour cannot be all-encompassing. This is because of the complex and unpredictable 602 nature of the human behaviour itself [50]. For example, in Section 5.2.1 we used the 603 Weidmann's equation to model the velocity of the people; however, it is known that this 604 equation varies across cultures [7, 50]. In addition, the usage of Wi-Fi at public events 605 is likely to change over time, a ecting the method parameters. Thus, it is important 606 to note that anytime a monitoring system is deployed at a di erent venue, a separate 607 calibration process is necessary in principle [7]. Furthermore, one can also include the 608 map of the venue in the calculations, so that the individual probability distributions can 609 be re-normalized to cover only accessible regions. 610 In this paper we discuss estimating crowd density to be able to detect critical density 611 and prevent crowd disasters. We note that planning optimal evacuation and navigation 612 of the crowd are separate research challenges and we refer the reader to [3] for a re- 613 cent overview. Concretely in our case the crowd can be navigated using the large TV 614 screens already present at the stadium, or apps that use the built-in compasses of the 615 smartphones. 616 20 Other possibilities for future work include: integrating Wi-Fi based crowd analysis 617 with video-based analysis, and investigating the performances under edge crowd scenarios, 618 like a crowd crush, and under di erent scenarios, like a crowd of pedestrians instead of 619 concert visitors. A real-time implementation of our prototype model is also important 620 future work. 621 622 7 List of abbreviations 623 RFID Radio-frequency identi cation MAC Media Access Control GDPR General Data Protection Regulation RSS Received signal strength AP Access Point pdf Probability density function GPS Global Positioning System624 8 Declarations 625 8.1 Ethics approval and consent to participate 626 Not applicable. 627 8.2 Consent for publication 628 Not applicable. 629 8.3 Availability of data and materials 630 By the Dutch law for privacy protection [20] we are not allowed to reveal data from 631 Sensation about the individual devices (even if they are anonymized), but we are allowed 632 to reveal aggregated results about the crowd. The dataset used in Section 5.3 is available 633 athttp://indoorlocplatform.uji.es/databases/get/1/ . All software (in Python) 634 related to the research results presented in this paper can be found in [21] ( https: 635 //doi.org/10.5281/zenodo.592479 ). Latest version of the software is available at: 636 https://github.com/ArenA-Crowds/Crowds . 637 8.4 Competing interests 638 The authors declare that they have no competing interests. 639 8.5 Funding 640 The project has been funded by the Netherlands eScience Center. 641 21 8.6 Authors' contributions 642 The authors contributed to the research and manuscript with the order they appear, the 643 last author being also the project principal investigator. All authors discussed the nal 644 results as well as read, improved, and approved the nal manuscript. In particular, SG 645 drafted the main ideas of the proposed solution, and wrote/performed the work presented 646 in sections 1, 2.2.2, 3, 4.1 (creating statistical ensembles), 4.2, 4.3, 5.1, 5.2, 5.3, 6 and 647 7. PR implemented the proposed method and wrote sections 2.2.1, 4.1 (Computing 648 individual probability distributions) and, jointly with SG, 5.4. JA wrote Section 2.1. RB 649 and BdV performed experiments and discussed the simulation strategy jointly with SG. 650 All authors actively participated in discussions also in the research phase. 651 8.7 Acknowledgments 652 We thank the Amsterdam ArenA (presently Johan Cruij ) Innovation Center for their 653 support. We are very thankful to Gijs van den Oord for reviewing the manuscript and 654 for his insights that led to improvement of the manuscript. We are very thankful to 655 Jessey Helslijnen, for his dedication to help with the video data and his unconditional 656 enthusiasm for contributing to science. 657