Analysis of Sparse Trajectory Features Based on Mobile Device Location for User Group Classification Using Gaussian Mixture Model

Abstract

Location data collected from mobile devices via global positioning system often lack semantic information and can form sparse trajectories in space and time. This study investigates whether user age groups can be accurately classified solely from such sparse spatial–temporal trajectories. We propose a feature extraction method based on a Gaussian mixture model (GMM), which assigns representative points (RPs) by clustering the location data and aggregating user trajectories into these RPs. We then construct three machine learning (ML) models—support vector classifier (SVC), random forest (RF), and deep neural network (DNN)—using the GMM-based features and compare their performance with that of the improved DNN (IDNN), which is an existing feature extraction approach. In our experiments, we introduced a missing value ratio ${\theta }_{\mathrm{th}}$ to quantify trajectory sparsity and analyzed the effect of trajectory sparsity on the classification accuracy and generalizability performance of the ML models. The results indicate that GMM-based features outperform IDNN-based features in both classification accuracy and generalization performance. Notably, the RF model achieved the highest accuracy, whereas the SVC model displayed stable generalizability. As the missing value ratio ${\theta }_{\mathrm{th}}$ increases, the IDNN becomes more susceptible to overfitting, whereas the GMM-based approach preserves accuracy and robustness. These findings suggest that sparse trajectories can still offer meaningful classification performance with appropriate feature design and model selection even without semantic information. This approach holds promise for domains where large-scale, sparse trajectory data are common, including urban planning, marketing analysis, and public policy.

1. Introduction

Over the first two decades of the twenty-first century, the use of mobile devices, such as smartphones with global positioning system (GPS) functions, has spread rapidly. Consequently, location data from these mobile devices via the GPS and the attribute data of users attached to these location data are now being collected on a grand scale, with the expectation that these data would be utilized for urban planning and the analysis of human activity, among other applications [1]. Location data and their associated data are regarded as the trajectories of mobile device users. There has been much research on classification and prediction models for user trajectories [2,3,4,5]. This study mainly focuses on classifying user groups based on location data and their associated data, as obtained from mobile devices.

Many studies have aimed at improving the accuracy of trajectory classification using location information consisting of spatiotemporal dimensions obtained from mobile terminals and associated semantic information [6,7,8,9,10,11]. In particular, Ferrero et al. proposed a MasterMovelet that characterizes trajectories with semantic dimensions using a series of uniform procedures [12]. However, whereas these studies assume that the semantic dimension is attached to the trajectory, the semantic information associated with the location data acquired from mobile devices is often missing. In this context, methods for trajectory similarity computation that do not require semantic dimensionality information other than spatial–temporal dimensionality have been proposed by Han et al. and Zhou et al. [13,14]. In addition, Endo et al. proposed a stacked denoising autoencoder (SDA) [15] based method for extracting trajectory features for transportation mode classification [16]. However, these methods do not assume that trajectories are sparse in the spatiotemporal dimension. Because irregular events obtain mobile terminal-based location information, the trajectories of such information could be sparser in spatial–temporal space.

Many studies have focused on trajectory classification using dense or semantically rich data. However, location data generated by mobile devices often lack semantic information and exhibit spatial–temporal sparsity. This raises the challenge of effectively utilizing sparse trajectories for applications such as user group classification and mobility analysis. Addressing this issue is crucial for harnessing the potential of trajectory data in fields such as urban planning, marketing analysis, and public policy. In addition, increasing data integrity can significantly enhance the usability of location data for these applications.

For these reasons, in our past work, we developed a sparse trajectory feature based on spatial–temporal information for classifying user age groups [17]; however, we conducted only a superficial analysis for limited binary classification and did not significantly examine feature availability. Furthermore, this past research focused solely on analyses based on the support vector classifier (SVC) [18,19] model and did not include comparisons with other machine learning (ML) models and trajectory features based on existing methods. Therefore, this study provides a detailed analysis of whether user groups can be estimated using only a sparse spatial–temporal trajectory without considering prior knowledge and past historical data associated with mobile device users. In particular, we analyze the feature based on [17], which is extracted from users’ location data using point-type dynamic population (PTDP) data [20] obtained from the GPS function of mobile devices. Moreover, we adopt multiple types of machine learning (ML) models that learn the features for six-class classification and conduct a detailed analysis. Additionally, by evaluating the user age group classification performance of the ML models, we verify the effectiveness of features based on sparse spatial–temporal trajectories in user age group estimation.

Each PTDP data point has an ID for identifying mobile devices by day, and those linked to this ID are the associated attributes registered by the user of the device and the location data from the device. Although location data are acquired in one-minute units, most of the data for a single day are missing because the data acquisition timing is determined by an event, such as the user opening a specific application of the device. Thus, we use a simple method to impute the missing location data and regard these imputed data as the user trajectory to make the data length uniform. Furthermore, for this study, we select age groups as the specific user groups that are subjected to classification and assume that the age groups can be categorized into six classes: 0 to 19, 20 s, 30 s, 40 s, 50 s, and 60 s and over.

Typically, mobile users’ location data, which consist of longitude and latitude, are expressed in the geographic coordinate system. In other words, the location data that are part of a trajectory are represented as points in the region [180° W–180° E, 90° S–90° N]. Because using such trajectories as features is not realistic, we need to summarize a trajectory into feature vectors consisting of a finite number of dimensions. Hence, we generate n-dimensional feature vectors by assigning n representative points (RPs) to a target region and aggregating the user locations into RPs. A RP is generated using the Gaussian mixture model (GMM), which assumes a Gaussian mixture distribution (GMD) for the population distribution in a target region. In particular, the GMD is generated from the raw location data of all users in a target period, and we regard the centroids of each Gaussian distribution as the RPs for that distribution. Thus, the RPs are selected in locations where the population is concentrated. In this way, each user trajectory is aggregated into RPs, and n-dimensional feature vectors for classifying user age groups are generated. We employed ML models that use the n-dimensional feature vector as input and output an age group linked to the mobile users. For the ML models, we employed SVC, random forest (RF) [21], and deep neural network (DNN) [22] models, and for the hyperparameter tuning of these ML models, we used fivefold grid search cross-validation (GS+CV). With feature vectors generated using GMD in this manner, it is possible to construct a classification model that uses only spatial–temporal information without complicated preprocessing, which differs for each target region. In addition, this study aims to verify the effectiveness of features based on the GMM for classifying sparse trajectories without semantic information. To achieve this goal, a comparison is conducted with existing feature extraction methods for trajectory data. In particular, the comparison focuses on features extracted using the improved DNN (IDNN), developed by Endo et al. [16], which is based on an SDA [15]. The IDNN does not assume that the trajectory is sparse, but part of this feature extraction method can be applied to our problem settings. This comparison demonstrates the superiority of the GMM-based features for sparse trajectories proposed in our past work, highlighting their effectiveness in this context. This allows for an evaluation of the relative performance and applicability of the GMM-based features.

In our experiments, we demonstrate the performance of three ML models trained with GMM-based features via PTDP data acquired between January and June 2021 in Narashino, Chiba, Japan. Because the number of GMD components that determine the feature vector’s dimensions can significantly influence the classification accuracies of the ML models, we conduct a sensitivity analysis on the number of components. In addition, we evaluate the accuracy and generalization performance of the combinations of three ML models and GMM- and IDNN-based features with respect to the level of missing value ratio, which quantifies the sparsity of trajectories. The trajectory classification ability of GMM-based features consisting of the sparse trajectory and the spatial–temporal sparsity of the user trajectory are examined in terms of their effectiveness for classification accuracy via a comparison with IDNN-based features. As a result, this study demonstrates that GMM-based features exhibit greater robustness to trajectory sparsity and higher generalization performance than the features extracted by the IDNN. These findings highlight the superior capability of GMM-based features in effectively classifying user age groups from sparse trajectories, compared to existing trajectory-based features, making them a promising approach for scenarios where only spatial–temporal information is available.

2. Literature Review

This section reviews relevant research studies to highlight the originality and practicality of our GMM-based feature extraction method. In this work, we treat age group classification as an example of user group classification using a sparse trajectory.

Traditionally, age group classification tasks have relied heavily on human biometric data. For instance, various age estimation techniques based on facial images have been developed over the years [23,24,25]. Additionally, Saxena et al. introduced an age group estimation method leveraging fingerprint data, while Nabila et al. proposed using gait patterns for the same purpose [26,27]. Beyond biometric-based approaches, alternative methods have also emerged. For example, Guimaraes et al. developed a strategy for classifying age attributes based on user-generated text on social media platforms [28].

Despite this wealth of research on age group classification, direct applications of trajectory data to age group classification are rare. Thus, we focus on research studies that, although not explicitly designed for user age group classification, can be adapted for this purpose. This section is organized into two key aspects: independence from semantic information and applicability for sparse trajectories, with potential applications and limitations of each method being discussed in detail.

2.1. Trajectory Analysis Using Semantic Information

Many studies have proposed trajectory classification methods that utilize data with both spatial–temporal and other semantic dimensions to increase accuracy. Wang et al. proposed an estimation method for the demographics of user groups based on trajectories obtained by analyzing wireless access points and prior knowledge of mobile users [29]. Montasser et al. also proposed a method for predicting the demographics of geographic units by analyzing geotagged tweets [30]. Many other methods have also been proposed, corresponding to classification or prediction models for user groups based on additional geographic and social network information [6,7,8,9,10,11]. However, these methods must process a considerable amount of prior knowledge data linked to mobile users and perform text analyses on posts in social network services.

In this context, Ferrero et al. proposed MasterMovelets, which designs features for trajectories with semantic dimensions, excluding spatial–temporal dimensions [12]. MasterMovelets can be used to examine the features of spatial–temporal trajectories with rich information regarding other semantic dimensions via a series procedure. This method increases classification accuracy by focusing on local subtrajectories in a target trajectory. Most of the research introduced so far assumes that the trajectories have semantic dimensions in addition to spatial–temporal dimensions and aims to increase the classification accuracy for such trajectories. However, user information linked to mobile device location data tends to have many missing values because of the nonuniformity of users’ registration statuses. Therefore, given that homogeneous information is not always added to all user trajectories, these methods cannot be applied to trajectories, which tend to have missing information.

2.2. Trajectory Analysis Using Only Spatial–Temporal Information

2.2.1. Similarity-Based Methods

DTW [31], a traditional similarity measurement method for sequential data, can measure the similarity between trajectories on the basis of spatiotemporal information in the context of trajectory classification. The similarity of trajectories can also be utilized for classifying user attributes associated with trajectories. Although a variety of extended DTWs have been proposed [32,33], these methods have a high computational cost. By contrast, Han et al. proposed a highly accurate prediction method for the similarity of trajectory patterns composed of only spatial–temporal information based on a graph-based approach [13]. Zhou et al. developed a noise-robust method for calculating the similarity of trajectories and solved data sparsity for the search space for the same problem setting [14]. However, these methods do not assume spatial–temporal sparsity of the trajectory based on mobile device location data. In contrast, location data from mobile devices are obtained almost entirely through irregular events, such as the user opening an application. In addition, with regard to anonymity, location data often do not account for microscopic movement patterns associated with users’ residences. The trajectories consisting of such location data are very sparse.

2.2.2. Trajectory Prediction and Recovery for Sparse Trajectory

Wang et al. proposed a model for predicting the future mobility of mobile users using such sparse trajectories based only on spatial–temporal dimensions [34]. Their model learns the trajectories of heterogeneous users accumulated in the past and predicts the future movements of users on the basis of their similarity with past users’ data. However, that study focused only on predicting spatial–temporal information about how users move from one point to another; it did not aim to estimate a user group without spatial–temporal information.

If it is possible to recover a complete trajectory for a sparse trajectory that lacks information, conventional methods could be applied to user-group classification. Various methods for such trajectory recovery have been proposed [35,36,37,38]. Xia et al. developed an attentional neural network-based model that recovers sparse trajectories by learning past accumulated historical trajectory data of users [35]. Sun et al. also proposed a trajectory recovery model for similar problem settings using a graph neural network, and both models demonstrate high performance; however, these models do not assume a new user who does not store past historical data [36]. For this problem, Si et al. proposed a new method that mitigates the cold start problem of such new users by using a transformer [37]. However, these methods require large amounts of historical data linked to specific users. Conversely, since unique IDs of mobile users’ location data covered by this study are reset on a daily basis, these methods cannot be applied to our problem settings.

2.3. Feature Extraction Applicable to Sparse Trajectories Without Semantic Information: A Comparative Baseline

Endo et al. proposed a feature extraction method for classifying transportation modes of trajectories using a SDA [15,16]. In this method, trajectories are converted into images and flattened into vectors, which are then used as inputs for a DNN constructed via SDA. Fine-tuning is performed using annotated labels, and the resulting improved DNN (IDNN) is employed to extract features from the trajectories. The semantic information associated with the trajectories is subsequently added, enabling high-accuracy transportation mode classification. In this study, we target trajectories that lack semantic information; however, the feature extraction method using the IDNN can still be applied. On the other hand, because IDNN-based feature extraction relies heavily on the shape of trajectories, it assumes dense trajectory data. Therefore, its performance on sparse trajectories remains uncertain. In our experiments, we compare the performance on sparse trajectories using features extracted by our proposed method and those extracted by the IDNN.

3. Trajectory Data

3.1. PTDP Data

The PTDP data used in this study are obtained from the GPS function of mobile devices and are provided by Agoop Corp. (Tokyo, Japan) [20]. Agoop Corp., the provider of these data, is a company that offers location-based services and is 100% owned by Softbank, one of the major telecommunications carriers in Japan. PTDP data are used for various purposes, such as urban planning, marketing analysis, and public policy development [39,40,41,42]. Accurately estimating user groups, e.g., age groups, can significantly enhance the usability of PTDP data in these research areas, enabling more precise analyses and evidence-based decision making. The columns of the PTDP data table, which are used in this research study, consist of the following attributes:

• Daily ID: User ID associated with the mobile device. This ID is reset on a daily basis.
• Timestamp: timestamp recorded in yyyy-mm-dd HH:MM for obtaining location information.
• Longitude: longitude of the location data.
• Latitude: latitude of the location data.
• Age group: age groups composed of those aged under 15, 70 and over and other age intervals divided into five-year increments, i.e., 0–14, 15–19, 20–24, …, 65–69, 70 and over.

Each row in the PTDP data table records the above attributes, which represent a user’s location information, demographic attributes, and time. For each user, daily location information is recorded across multiple rows. These point data are acquired through specific events, such as when a user opens an application on a smartphone. Each row of data is logged at a minimum interval of one minute. Therefore, by arranging these location data in a time series order, trajectories associated with individual users can be constructed.

Notably, there are missing values for the age group attribute; however, in this study, we excluded data with missing age group information from the experiments. Using the age groups in their original form could result in an excessive number of classes, which might lead to an extremely small number of data points in each class. Moreover, the original classification does not clearly indicate the attributes of individuals in each class, such as whether they are students or workers. For these reasons, this study reclassifies the age groups into six categories: 0–19, 20 s, 30 s, 40 s, 50 s, and 60 s and over. These groups can generally be associated with standard attributes, such as students, early-career workers, young workers, middle-career workers, senior workers, and retirees. Notably, the behavior patterns of the 0–19 age group may overlap with those of their caregivers.

3.2. Data Assumption

As mentioned in Section 3.1, this study assumes that PTDP data are the location data that constitute the trajectories. The set of IDs linked to the mobile devices (where each ID is linked to exactly one mobile device) is denoted by K, and the set of location data corresponding to ID $k\in K$ is denoted by $L_k$. Note that all the IDs are reset on a daily basis. If the maximum length of the location data in a day is denoted by $\mathcal{T}$, the location data size for a day is $|L_k|\le \mathcal{T}$, i.e.,

$$\begin{array}{c}\hfill L_k=\{{\mathit{x}}_{k,t}\mid t\in {\mathbb{N}}_{\le \mathcal{T}}\},\left|L_k\right|\le \mathcal{T},\end{array}$$

where ${\mathbb{N}}_{\le \mathcal{T}}$ is a natural number that is less than or equal to $\mathcal{T}$. Note that ${\mathit{x}}_{k,t}$ is a point of ID k on the geographic coordinate system at time t, and ${\mathit{x}}_{k,t}={(x_{k,t}^{\mathrm{lon}},x_{k,t}^{\mathrm{lat}})}^{\top },x_{k,t}^{\mathrm{lon}}\in$ [180° W–180° E], $x_{k,t}^{\mathrm{lat}}\in$ [90° S–90° N].

Although the acquisition timing of location data ${\mathit{x}}_{k,t}$ is at a unit time $t_{\mathrm{unit}}$ interval, if a mobile device with ID k does not cause any event, location data are not acquired. Hence, for almost all sets $L_k$ for $\mathcal{T}$, a large portion of the location data acquired by mobile devices are missing, i.e., $|L_k|\ll \mathcal{T}$.

3.3. Imputation for Missing Location Data

As mentioned in Section 3.2, a large portion of the location data for any $L_k$ is missing. Accordingly, generating features is difficult because the data length of $L_k$ differs from that of ID k. Thus, we must impute a location dataset $L_k$ to a trajectory $T_k=\{{\mathit{\tau }}_{k,1},{\mathit{\tau }}_{k,2},\dots ,{\mathit{\tau }}_{k,\mathcal{T}}\}$ with aligned data length, where ${\mathit{\tau }}_{k,t}$ is a point of ID k on the geographic coordinate system at time t.

As we mentioned in Section 2.2, although various trajectory recovery methods have been proposed, these models require large amounts of historical trajectory data linked to unique user IDs. However, because the mobile location data covered by this study do not store historical data for user anonymization, we cannot use such historical data. Thus, in this study, we use simple linear interpolation for trajectory recovery.

Let us consider a scenario in which an ID is obtained at Site A at time $t_1$. Then, later, the same ID is acquired at Site B at time $t_2$, where $t_1<t_2$ and $t_1,t_2\in \mathbb{N}$. When $t_2-t_1>1$, we assume that the ID contains missing location data between $t_1$ and $t_2$. An assumption needs to be made regarding the missing data to be imputed if there is no information about the ID location. Hence, we assume that a mobile user departs from Site A at $t_1$ and then arrives at Site B at $t_2$. In other words, on the basis of the unit time $t_{\mathrm{unit}}$ and the acquisition times $t_1$ and $t_2$, the path can be obtained by dividing the interval between Sites A and B by $(t_1-t_2)/t_{\mathrm{unit}}$. By this assumption, we obtain trajectories $T_k,\forall k\in K$.

The missing value ratio of $T_k$, i.e., the ratio of imputed data for the entire $T_k$, is defined by the function $\Theta (·)\mapsto [0,1)$ with location dataset $L_k$, which is the source of $T_k$ as an input, as follows:

$$\begin{array}{c}\hfill \Theta (\mathcal{T},L_k)={\theta }_k={\displaystyle \frac{\mathcal{T}-|L_k|}{\mathcal{T}}}.\end{array}$$

(1)

Note that $|L_k|>0$. Afterward, denoting the missing value ratio by ${\theta }_k$, we use the function $\Theta$ with $L_k$ as the input. The missing value ratio ${\theta }_k$ refers to the sparsity of trajectory $T_k$ such that the closer the value of ${\theta }_k$ is to 1, the sparser the trajectory $T_k$.

Figure 1 shows the distribution of the missing value ratios of the trajectories based on PTDP data obtained at Narashino, Chiba, from January to June 2021. Figure 1 indicates that almost all the location data are ${\theta }_k\ge 0.95$. Notably, the data with ${\theta }_k\ge 0.99$ comprise 63% of the data. Figure 1 thus implies that almost all the trajectories based on PTDP data are sparse. Figure 2 visualizes the imputation.

4. Feature Extraction Method for Sparse Trajectories

4.1. Feature Extraction by GMM

In our previous study, we proposed a method that determines and examines the features of trajectories by combining the imputed trajectory in Section 3.3 and the GMM in [17]. In this study, we analyze the features of trajectories via the same method; however, the properties of the method are related to those of the experiment presented in this paper and are therefore introduced in detail below.

The PTDP data consist of points denoted by $\mathit{x}={(x_{\mathrm{lon}},x_{\mathrm{lat}})}^{\top }$ in the geographic coordinate system. Because any $\mathit{x}$ is a set of coordinates in the region [180° W–180° E, 90° S–90° N], it is not realistic to regard its coordinates as features. Hence, n RPs are determined for the given region, and by allowing the PTDP data to aggregate into the RPs, we generate n-dimensional feature vectors ${\mathit{f}}_k^{\mathrm{GMM}}$ for each ID k.

4.1.1. Determining RPs via GMM

Generally, it is desirable to select densely populated locations, such as stations and commercial facilities, as RPs; however, this method requires an analysis of the regional characteristics of a target area and thus would not result in a general framework. Instead, we use the GMM to aggregate a target area into a number of RPs. The GMM is a clustering method that represents data points scattered in space as overlapping n-Gaussian distributions.

In this study, because the target space is two-dimensional, the GMD is represented by an overlapping bivariate Gaussian distribution as follows:

$$\begin{array}{c}\hfill N(\mathit{x}\mid \mathit{\mu },\mathbf{\Sigma })={\displaystyle \frac{1}{{2\pi \left|\mathbf{\Sigma }\right|}^{1/2}}}exp\left\{-{\displaystyle \frac{1}{2}}{(\mathit{x}-\mathit{\mu })}^{\top }{\mathbf{\Sigma }}^{-1}(\mathit{x}-\mathit{\mu })\right\},\end{array}$$

where $\mathit{\mu }$ and $\mathbf{\Sigma }$ are the mean vector and variance–covariance matrix, respectively. Thus, the GMD is defined as follows:

$$\begin{array}{cc}\hfill p\left(\mathit{x}\right)& =\sum _{i\in I}{\pi }_{i}N(\mathit{x}\mid {\mathit{\mu }}_i,{\mathbf{\Sigma }}_i),\hfill \\ \hfill N(\mathit{x}\mid {\mathit{\mu }}_i,{\mathbf{\Sigma }}_i)& ={\displaystyle \frac{1}{2\pi |{\mathbf{\Sigma }}_i{|}^{1/2}}}exp\left\{-{\displaystyle \frac{1}{2}}{(\mathit{x}-{\mathit{\mu }}_i)}^{\top }{\mathbf{\Sigma }}_i^{-1}(\mathit{x}-{\mathit{\mu }}_i)\right\},\hfill \\ \hfill \sum _{i\in I}{\pi }_i& =1,\hfill \end{array}$$

where I is the set of bivariate Gaussian distributions and ${\mathit{\mu }}_i$ and ${\mathbf{\Sigma }}_i$ are the i-th mean vector and variance–covariance matrix, respectively. In addition, ${\pi }_i$ is the mixing coefficient that expresses the weight of the i-th Gaussian distribution. The obtained means ${\mathit{\mu }}_i,i\in I$ of n bivariate Gaussian distributions are regarded as the coordinates of the RPs with which we approximate the target space. The number of overlapping Gaussian distributions, i.e., the number of components n, can be determined arbitrarily. In this context, because the lengths of the location data are not necessarily the same for any ID k, we use the raw location data $\{\mathit{x}\mid \mathit{x}\in L_k,\forall k\in K\}$ as samples for the GMD.

4.1.2. Aggregation of User Trajectories into RPs

As mentioned at the beginning of this section, using trajectory $T_k$ of ID k imputed based on Section 3.3 as feature vectors is not realistic. Accordingly, we define the feature vector to aggregate trajectory $T_k$ of data length $\mathcal{T}$ into RPs obtained via the GMM. In particular, the nearest RP for each coordinate in trajectory $T_k$ is determined, and the obtained RP is assigned as the representative for that coordinate. The number of times that a coordinate is assigned to a certain RP is stored in the RP and is regarded as the stay time of ID k at that RP. This operation generates feature vectors ${\mathit{f}}_k^{\mathrm{GMM}},\forall k\in K$ with n features based on trajectory $T_k,\forall k\in K$.

The specific procedures are indicated below. Let the coordinates of ID k obtained at time t be ${\mathit{x}}_{k,t}=(x_{k,t}^{\mathrm{lat}}$, $x_{k,t}^{\mathrm{lon}})$. The cluster variable is denoted by $z_i$, which equals 1 if sample ${\mathit{x}}_{k,t}$ belongs to the Gaussian distribution $i\in I$; otherwise, it equals 0. In this context, the probability that $z_i=1$, given ${\mathit{x}}_{k,t}$, is $p(z_i=1\mid {\mathit{x}}_{k,t})$, and from Bayes’ theorem,

$$\begin{array}{cc}\hfill p(z_i=1\mid {\mathit{x}}_{k,t})& ={\displaystyle \frac{p(z_i=1)p({\mathit{x}}_{k,t}\mid z_i=1)}{{\sum }_{i\in I}p(z_i=1)p({\mathit{x}}_{k,t}\mid z_i=1)}}\hfill \\ \hfill & ={\displaystyle \frac{{\pi }_{i}p\left({\mathit{x}}_{k,t}\right)}{{\sum }_{i\in I}{\pi }_{i}p({\mathit{x}}_{k,t}\mid z_i=1)}}\hfill \end{array}$$

is obtained. Here, because $p({\mathit{x}}_{k,t}\mid z_i=1)$ is the probability of sample ${\mathit{x}}_{k,t}$ on the i-th Gaussian discussion of the GMD,

$$\begin{array}{cc}\hfill p(z_i=1\mid {\mathit{x}}_{k,t})& ={\displaystyle \frac{{\pi }_{i}N({\mathit{x}}_{k,t}\mid {\mathit{\mu }}_i,{\mathbf{\Sigma }}_i)}{{\sum }_{i\in I}{\pi }_{i}N({\mathit{x}}_{k,t}\mid {\mathit{\mu }}_i,{\mathbf{\Sigma }}_i)}}\hfill \end{array}$$

(2)

is obtained based on

$$\begin{array}{c}\hfill p({\mathit{x}}_{k,t}\mid z_i=1)=N({\mathit{x}}_{k,t}\mid {\mathit{\mu }}_i,{\mathbf{\Sigma }}_i).\end{array}$$

Equation (2) is the conditional probability that $z_i=1$ given sample ${\mathit{x}}_{k,t}$; thus, the Gaussian distribution i that maximizes Equation (2), i.e.,

$$\begin{array}{c}\hfill {\mathrm{argmax}}_{i\in I}p(z_i=1\mid {\mathit{x}}_{k,t}),\end{array}$$

is the distribution to which sample ${\mathit{x}}_{k,t}$ belongs. By applying this operation to ${\mathit{x}}_{k,t},\forall t\in {\mathbb{N}}_{\le \mathcal{T}}$, we can transform $T_k$ of ID k to an n-dimensional feature vector ${\mathit{f}}_k^{\mathrm{GMM}}$ and generate a training dataset for the ML model. Note that the components of vector ${\mathit{f}}_k^{\mathrm{GMM}}$ are assumed to be $f_i,i\in I$, $f_i\in \{x\in {\mathbb{R}}_{\ge 0}\mid x\le \mathcal{T}\}$ and ${\sum }_{i\in I}{f}_i=\mathcal{T}$. Because this n-dimensional feature vector is obtained by aggregating the target area into n RPs, the size of n is regarded as the expressiveness for the source area and may affect the classification accuracy.

A diagram of this procedure is presented in Figure 3.

4.2. Deep Feature Extraction: A Comparative Baseline

Endo et al. proposed a feature extraction method for classifying the transportation modes of trajectories obtained via GPS data [16]. In this method, the target trajectories for classification are converted into images, and fine-tuning is performed on an nueral network (NN) constructed via SDA with annotated labels. The IDNN obtained through this process is then used to extract features from the trajectories. Furthermore, the extracted features are combined with semantic information associated with the trajectories to generate new features.

Alternatively, we aim to examine whether user age group classification can be performed using sparse trajectory data without semantic information. Although the entire method proposed by Endo et al. [16] requires semantic information, the IDNN, which is a part of their method, does not rely on such information. Therefore, the IDNN is considered applicable to the nonsemantic and sparse trajectory data targeted in this study. In this study, we compare the classification accuracy for sparse trajectories using feature extraction via the IDNN with that via the GMM-based method.

The feature extraction method using the IDNN generates features through the following three phases:

• Trajectory image generation
  The target area is divided into rectangular regions of a specified size $W_m\times H_m$, and each region is represented as a pixel. This allows the target area to be treated as an image, where each element of the trajectory $T_k$ is converted into a pixel value $v(i,j)$ for the region $(i,j),i\in \{0,1,\dots ,W_m\},j\in \{0,1,\dots ,H_m\}$. The pixel value $v(i,j)$ is incremented by one for every minute the user k spends in the region. Thus, $v(i,j)$ represents the total time that user k spends in region $(i,j)$, satisfying the following equation:

  $$\begin{array}{c}\hfill \sum _i^{W_m}\sum _j^{H_m}v(i,j)=\mathcal{T}.\end{array}$$
  Finally, the pixel values $v(i,j)$ are normalized such that they fall within the range $\widehat{v}(i,j)\in [0,1]$. This process generates a grayscale image $I_k=\left(\widehat{v}(i,j)\right)\in {[0,1]}^{W_m\times H_m}$ associated with each ID k.
• Pretraining by SDA
  The autoencoder (AE) is a type of NN designed to compress and reconstruct input vectors, and it is primarily used for feature extraction and dimensionality reduction. For example, when a flattened vector of an image $I_k$ is denoted as ${\mathit{x}}_k^{\mathrm{img}}\in {[0,1]}^{W_m{H}_m}$, an AE constructs an NN that minimizes the loss function $\mathcal{L}({\mathit{x}}_k^{\mathrm{img}},\mathit{z})$ between ${\mathit{x}}_k^{\mathrm{img}}$ and the reconstructed vector $\mathit{z}\in {\mathbb{R}}^{W_m{H}_m}$. The denoising autoencoder (DAE) is an extension of the AE that enhances robustness by adding noise, sampled from a distribution $q_D$, to the input vector ${\mathit{x}}_k^{\mathrm{img}}$ before constructing the AE. SDA is a type of DNN built by stacking multiple DAEs as middle layers, where each layer is trained independently by its respective DAE [15]. Specifically, the l-th middle layer of an SDA is trained by minimizing the loss function $\mathcal{L}({\mathit{x}}_{l-1},{\mathit{z}}_l)$ between the features ${\mathit{x}}_{l-1}$ learned by the $(l-1)$-th layer DAE and the output ${\mathit{z}}_l$ of the l-th layer. The number of middle layers $\lambda$ in a DNN constructed by an SDA equals the number of stacked DAEs. In the l-th DAE, the input vector ${\mathit{x}}_{l-1}$ is first passed through the encoding function $\phi (·)$ after noise is added according to $q_D$ and then reconstructed as a vector ${\mathit{z}}_l$ through the decoding function ${\phi }^{\prime }(·)$. Typically, $q_D$ employs masking noise, which randomly forces selected elements to zero. The conceptual diagram of an SDA is shown in Figure 4a.
• Fine tuning and feature extraction
  Fine-tuning is performed by adding an output layer to the top of the pretrained DNN constructed by the SDA and minimizing the loss function ${\mathcal{L}}_{\mathrm{supv}}(\mathit{y},{\mathit{z}}_{\mathrm{supv}})$ between the output ${\mathit{z}}_{\mathrm{supv}}$ and the target labels $\mathit{y}$ (Figure 4b). In the fine-tuning process, the target labels are defined as binary vectors representing the classes. Finally, the fine-tuned IDNN takes the vector ${\mathit{x}}_k^{\mathrm{img}}$ as input, and the feature vector ${\mathit{f}}_k^{\mathrm{IDNN}}$ is extracted from the $\lambda$-th intermediate layer of the IDNN. The dimensionality of the extracted feature vector ${\mathit{f}}_k^{\mathrm{IDNN}}$ is equal to the number of neurons in the $\lambda$-th middle layer.

Through this series of procedures, the feature vectors ${\mathit{f}}_k^{\mathrm{IDNN}},\forall k\in K$ are obtained. Notably, the IDNN proposed by Endo et al. [16] utilizes not only the feature vectors derived from trajectories but also other semantic information as features. By contrast, as this study aims to classify trajectories without semantic information, the classification performance is evaluated using only the feature vectors ${\mathit{f}}_k^{\mathrm{IDNN}}$ obtained through the IDNN.

5. Experimental Setup

5.1. Data Summary

PTDP data from January to June 2021 from Narashino, Chiba, Japan, provided by Agoop Corp. [20], were used for this experiment. The target periods of data extraction comprise only workdays because people’s behaviors on holidays are commonly different from those on workdays, which could affect classification accuracy. Narashino is in the region [139.59–140.50° E], [35.38–35.42° N], covering an area of 20.27 [km²], and as of January 2021, its population was 176,197 [43,44]. The location data acquisition timing is unit time $t_{\mathrm{unit}}=1$ [min]. Because IDs linked to mobile devices are reset at 12:00 a.m., each user’s behavior can be identified only for a maximum of one day. In other words, the data length of trajectories $T_k,\forall k\in K$ is $\mathcal{T}=60\times 24=1440$. Each ID k is associated with one age group among five-year age divisions. Although an age group is also linked to each mobile device, not all IDs have age group information, and some data are missing. Therefore, in this study, we construct and evaluate the ML model using only data for IDs with age group information.

This study aims to determine whether the features generated by the sparse trajectory $T_k$ can lead to accurate classifications. To facilitate this evaluation, we simplify the problem into a six-class classification problem, with each class comprising ages 0 to 19, 20 s, 30 s, 40 s, 50 s, and 60 s and over.

Because people’s behavior can vary depending on the season, we divided the data from January to June 2021 by month, i.e., $1,2,\dots ,6$. Figure 5 shows the distribution of the data size by month and class; the figure indicates that the dataset is imbalanced. Additionally, Figure 6 shows the distribution of the data size for the missing value rate ${\theta }_k,\forall k\in K$ for each month. Figure 6 indicates that the distributions for all the periods are almost the same shape.

Note that during the target period, the target area was under the outbreak of the novel oronavirus (COVID-19), and declarations of a state of emergency and a semistate of emergency were enforced from 8 January to 21 March and 20 April to 1 August 2021, respectively [45]. This is noted in the experimental results and discussion.

5.2. Under Sampling for Imbalanced Data

Typically, an ML model that learns imbalanced data has problems in that the estimation accuracy for each class is biased. Many methods for sampling and evaluating imbalanced data have been proposed [46]. In the experiments, the imbalanced dataset was balanced by undersampling to simplify the evaluation of the classification models.

Undersampling is the process of sampling data from all classes according to the class with the smallest data size. In other words, if datasets whose elements are pairs $(\mathit{f},y)$ of inputs $\mathit{f}$ and outputs y are denoted by $\mathcal{D}$, the data size of any month that should be sampled is as follows:

$$\begin{array}{c}\hfill S_m=\left|C\right|(\underset{c\in C}{min}\left|\{(\mathit{f},y)\mid y=c,\forall (\mathit{f},y)\in \mathcal{D}\}\right|),\forall m\in M,\end{array}$$

(3)

where C and M are the sets of classes and months, respectively.

Here, we must consider the differences in the data size by month, as shown in Figure 5. Because the data sizes $S_m$ obtained using Equation (3) are not necessarily equal for all $m\in M$, the sizes are imbalanced among months. The larger the data size is compared with the feature dimension, the better the robustness and accuracy of an ML model in general [47,48]. Each target period’s reliability gap in terms of classification accuracy may be significant if the data sizes differ. For this reason, the datasets for each month are undersampled to the class with the minimum data size to facilitate comparisons among the classification models for each month. That is, letting ${\mathcal{D}}_m,\left|{\mathcal{D}}_m\right|=S_m,m\in M$ be a dataset for month m, the data sizes for the months are calculated as follows:

$$\begin{array}{c}\hfill S=|C|\left(\underset{m\in M}{min}\underset{c\in C}{min}\right|\{(\mathit{f},y)\mid y=c,\forall (\mathit{f},y)\in {\mathcal{D}}_m\}\left|\right).\end{array}$$

Note that in this study, $M=\{1,2,\dots ,6\}$ and C = {‘0 to 19’, ‘20 s’, ‘30 s’, ‘40 s’, ‘50 s’, ‘60 s and over’}.

To construct and evaluate the ML models, we divided the dataset into training and test datasets. Because dividing a dataset into 70% training and 30% test datasets was shown to be desirable in many instances, the same was adopted for the training datasets in this study. In contrast, many test datasets can increase the validation accuracy of an ML model; thus, all extra datasets after training datasets were extracted and used as test datasets. Therefore, the number of training datasets was constant for any month, but the test datasets differed by month. If the training and test data sizes are denoted by $S^{\mathrm{train}}$ and $S_m^{\mathrm{test}}$, respectively, then their sizes are defined as follows:

$$\begin{array}{cc}\hfill S^{\mathrm{train}}& =⌊0.7S⌋,\hfill \\ \hfill S_m^{\mathrm{test}}& =S_m-S^{\mathrm{train}},\forall m\in M.\hfill \end{array}$$

The training and test data are sampled such that there are no duplications from the dataset of size $S_m$. Because randomness is present when the training and test data are extracted from the source dataset, we conduct undersampling for 20 random seeds in the range of 0–19. In this way, 20 datasets for each month are prepared.

The specific sizes S and $S_m$, which are the basis for the training and test data for each month, respectively, are shown in Figure 5.

5.3. ML Models for Classification

To evaluate the proposed features for sparse trajectories (Section 4.1) and the features obtained via the method described in Section 4.2, we employed multiple ML models, including SVC with an RBF kernel, RF, and DNN. The models were implemented and evaluated using scikit-learn version 1.5.2 [49].

Although the components of the feature vector ${\mathit{f}}_k^{\mathrm{GMM}},\forall k\in K$ mentioned in Section 4 are given by $f_i\in \{x\in {\mathbb{R}}_{\ge 0}\mid x\le \mathcal{T}\},\forall i\in {\mathbb{N}}_{\le \mathcal{T}}$, for calculation convenience, we use the feature vector to normalize the components into values between 0 and 1, as follows:

$$\begin{array}{c}\hfill {\mathit{f}}_k^{\prime \mathrm{GMM}}={\displaystyle \frac{1}{\mathcal{T}}}{\mathit{f}}_k^{\mathrm{GMM}}.\end{array}$$

Note that ${\sum }_{k\in K}{\mathit{f}}_k^{\prime \mathrm{GMM}}=1$. Similarly, each element of ${\mathit{f}}^{\mathrm{IDNN}}$ is normalized to the range ${\mathit{f}}^{\prime \mathrm{IDNN}}\in {[0,1]}^n$ via min–max scaling. Because each ML model has its own hyperparameters, a GS is applied to tune them by providing a specified search range for each model.

The SVC with an RBF kernel has two hyperparameters: the regularization parameter C and the kernel coefficient $\gamma$. These parameters are tuned over the following ranges:

$$C\in \{{10}^1,{10}^2,{10}^3\},\gamma \in \{{10}^{-3},{10}^{-2},{10}^{-1}\}.$$

The RF has several hyperparameters, but tuning all of them is impractical. Therefore, we focus on two key hyperparameters: the maximum depth $d_{max}$ and the number of trees $n_{\mathrm{trees}}$. These parameters are tuned over the following ranges:

$$d_{max}\in \{4,6,8\},n_{\mathrm{trees}}\in \{100,300,500\}.$$

The other hyperparameters of the RF are set to the default values provided by scikit-learn (version 1.5.2).

The DNN also has numerous hyperparameters. To simplify tuning, we select two of them: the hidden layer size $\mathit{h}$ and the regularization parameter $\alpha$. These parameters are tuned over the following ranges:

$$\mathit{h}\in \{\left(64\right),\left(128\right),(64,64)\},\alpha \in \{{10}^{-4},{10}^{-3},{10}^{-2}\}.$$

The DNN uses the Adam optimizer for parameter updates and ReLU as the activation function. All other hyperparameters are set to the default values provided by scikit-learn (version 1.5.2).

Moreover, we conducted fivefold CV for each parameter pair to prevent overfitting. This method is called GS+CV and is visualized in Figure 7.

5.4. Metrics

The classification accuracies of the learned ML models are evaluated on the test data. The metrics employed for the analysis are the following macro averages of accuracy, precision, recall, and F-measure, which are assumed to be $A_{\mathrm{M}},P_{\mathrm{M}},R_{\mathrm{M}}$, and $F{1}_{\mathrm{M}}$, respectively, and are defined as follows:

$$\begin{array}{cc}\hfill A_{\mathrm{M}}& ={\displaystyle \frac{1}{\left|C\right|}}\sum _{c\in C}{\displaystyle \frac{{\mathrm{TP}}_c+{\mathrm{TN}}_c}{{\mathrm{TP}}_c+{\mathrm{FP}}_c+{\mathrm{TN}}_c+{\mathrm{FN}}_c}},\hfill \\ \hfill P_{\mathrm{M}}& ={\displaystyle \frac{1}{\left|C\right|}}\sum _{c\in C}{P}_c,P_c={\displaystyle \frac{{\mathrm{TP}}_c}{{\mathrm{TP}}_c+{\mathrm{FP}}_c}},\hfill \\ \hfill R_{\mathrm{M}}& ={\displaystyle \frac{1}{\left|C\right|}}\sum _{c\in C}{R}_c,R_c={\displaystyle \frac{{\mathrm{TP}}_c}{{\mathrm{TP}}_c+{\mathrm{FN}}_c}},\hfill \\ \hfill F{1}_{\mathrm{M}}& ={\displaystyle \frac{1}{\left|C\right|}}\sum _{c\in C}{\displaystyle \frac{2P_c{R}_c}{P_c+R_c}},\hfill \end{array}$$

(4)

where ${\mathrm{TN}}_c,{\mathrm{FN}}_c,{\mathrm{FP}}_c$, and ${\mathrm{TP}}_c$ are the true negative, false negative, false positive, and true positive of each class c, respectively.

Because the dataset is a balanced six-class dataset, the accuracy can be expressed simply as follows:

$$\begin{array}{c}\hfill A_{\mathrm{M}}\left\{\begin{array}{cc}>{\displaystyle \frac{1}{6}}:\hfill & \mathrm{Performance}\mathrm{is}\mathrm{higher}\mathrm{than}\hfill \\ \hfill & \mathrm{that}\mathrm{of}\mathrm{random}\mathrm{model},\hfill \\ \le {\displaystyle \frac{1}{6}}:\hfill & \mathrm{Performance}\mathrm{is}\mathrm{equivalent}\mathrm{to}\mathrm{or}\hfill \\ \hfill & \mathrm{lower}\mathrm{than}\mathrm{that}\mathrm{of}\mathrm{random}\mathrm{model},\hfill \end{array}\right.\end{array}$$

That is, the approximate chance level of this problem is 17%. The dataset utilized in this study was sampled to ensure balanced data across all classes, as described in Section 5.2. For balanced data, there is no requirement to adjust for class imbalance when evaluating accuracy, which allows the metric defined in Equation (4) directly comparable to the chance level. We thus evaluate the model performance on the basis of accuracy. In addition, as mentioned in Section 5.2, given that 20 datasets were prepared to decrease randomness, 20 models are also constructed. In other words, 20 values of each metric are obtained for each parameter setting. Thus, the average and standard deviation can be obtained for the metrics of all the ML models.

6. Numerical Experiment

6.1. Methods

6.1.1. Tuning GMD Components n for Model Performance

As shown in Section 4.1.1, the dimension of the feature vector ${\mathit{f}}_k^{\mathrm{GMM}}$ generated based on trajectory $T_k$ is equal to the number n of GMD components. In general, too large the feature vector dimension for a sample size has a negative effect on the robustness and regularization of the model [47,48]. By contrast, for the feature vector obtained by aggregating a trajectory to RPs, the dimension is regarded as the expressiveness for a trajectory; thus, a larger dimension n might increase the classification accuracy. Therefore, the trends of the model performance by components n of the GMD are verified, and a search for the components n that provide the highest accuracy is performed. Let the search range of components n be $n\in \left\{2\right\}\cup \{10,20,\dots ,200\}$. The minimum dimensionality required for the feature vector is $n=2$ to ensure distinctiveness. This is because setting $n=1$ results in all feature vectors being ${\mathit{f}}_k^{\mathrm{GMM}}={\left(\mathcal{T}\right)}^{\top }$, eliminating any distinctiveness of trajectories. Note that the GMD is computed via the expectation–maximization algorithm [50].

The metric employed for the GMD is often either Akaike’s information criterion or the Bayesian information criterion. These metrics provide the log-likelihood of the GMD. However, in this study, the performance of the classification model trained with features based on the GMD is essential, whereas the likelihood of the GMD for the actual population distribution is not essential. Thus, we directly use the model performance, i.e., the metrics in Section 5.4 as evaluation indices for the GMD components n. By contrast, only the accuracy in Section 5.4 is used in the experiment because the objective is to determine the trends of model performance by components n and not the detailed behavior of the models. In particular, the average and standard deviation of the accuracy are calculated for the results of 20 seeds.

Because the numbers of datasets and component candidates are 20 and 21, respectively, the number of constructed models is $20\times 21=420$ for each ML method. The hyperparameter tuning described in Section 5.3 is conducted for all 420 models for each ML method.

6.1.2. Comparison of GMM and IDNN

The performances of the ML models that use the feature vectors ${\mathit{f}}_k^{\mathrm{GMM}}$ and ${\mathit{f}}_k^{\mathrm{IDNN}}$ are compared. The IDNN is trained via the same training and test datasets with ${\theta }_{\mathrm{th}}=1.00$ as the ML models in Section 6.1.1.

In [16], a clipping process is applied to center the trajectory centroid within the image, as their objective is to classify transportation modes, focusing primarily on the shape of the trajectory. By contrast, this study aims at attribute classification, where the relative position of the trajectory within the region, in addition to its shape, is crucial. Therefore, no clipping process is applied.

The IDNN has several hyperparameters that need to be tuned. It is not realistic to tune all the hyperparameters in terms of the computational cost; thus, we focus on the number of hidden layers $\lambda$ and the number of neurons in the l-th layer $n^{\left(l\right)}$ to tune these parameters via the GS+CV.

To simplify the process, we assume that all layers have the same number of neurons, i.e., $n^{\left(1\right)}=n^{\left(2\right)}=\cdots =n^{\left(\lambda \right)}$. We suppose the search ranges for these hyperparameters with $\lambda \in \{2,3,4\}$ and $n^{\left(l\right)}\in \{100,150,200\},\forall l\in \{1,2,\dots ,\lambda \}$. The other parameters and settings for the IDNN are as follows:

• Image size: $W_m=H_m=32$
• Noise distribution $q_D$: masking noise (rate = 0.2)
• Activation functions:

  −

  $\phi (·)$, ${\phi }^{\prime }(·)$: ReLU

  −

  ${\phi }_{\mathrm{supv}}(·)$: Softmax

• Loss functions:

  −

  $\mathcal{L}$, ${\mathcal{L}}_{\mathrm{supv}}$: mean squared error

• Optimizer for each DEA and fine tuning of the DNN: Adam

The image size determines the expressiveness of the trajectory representation. A larger image size is anticipated to enhance the representational capacity of the trajectory; however, it may also increase susceptibility to noise. Conversely, a smaller image size is more robust to noise but may compromise the expressive power of the trajectory representation. In [16], an image size of $W_m=H_m=40$ was reported to achieve the best accuracy and demonstrate good generalization performance. In this study, the task involves a higher degree of trajectory sparsity and noise compared to [16]. To address these challenges while avoiding significant loss in representational capability due to an overly small image size, the image size was set to $W_m=H_m=32$.

Masking noise is used during the pre-training phase of SDA to deliberately introduce noise, with the goal of improving the robustness of extracted features. In [16], it was reported that a masking noise rate of 0.4 achieved the highest accuracy and improved generalization performance. However, because the trajectories in this study inherently contain noise, using a high masking noise rate would significantly reduce accuracy. Consequently, the masking noise rate was set to 0.2 in this study.

In [16], sigmoid functions were used as activation functions. However, sigmoid functions are known to cause the vanishing gradient problem. In recent years, ReLU functions have gained popularity in many tasks [51]. Furthermore, since the output layer during fine-tuning is compared with binary vectors representing the classes, the Softmax function is more suitable than the sigmoid function. For these reasons, this study adopted ReLU for $\phi (·)$ and ${\phi }^{\prime }(·)$, and the Softmax function for ${\phi }_{\mathrm{supv}}(·)$. For the loss function, although Cross-Entropy is commonly used for classification tasks, some studies report that mean squared error (MSE) performs better for certain problems [52]. Since there is no definitive conclusion about the superiority of one over the other, this study follows [15,16] and adopts MSE as the loss function.

Regarding the optimizer, SGD was utilized in [15,16]. However, SGD requires careful tuning of the learning rate and is prone to becoming stuck in local optima if not configured correctly. On the other hand, Adam overcomes these limitations by automatically adjusting the learning rate. Considering the balance between computational cost and improved accuracy, this study selected Adam as the optimizer.

Based on these parameter settings, the average and standard deviation of the accuracy are calculated for the results of 20 seeds and compared with the results obtained in Section 6.1.1.

6.1.3. Performance by the Missing Value Rate of the User Trajectory

As in Section 3.3, the missing value rate ${\theta }_k$ of trajectory $T_k$ indicates its sparsity. Hence, the classification performance increases if only trajectories with low sparsity are learned by providing the threshold ${\theta }_{\mathrm{th}}$ for the missing value range ${\theta }_k,\forall k\in K$. Let $T_k^{{\theta }_{\mathrm{th}}}$ be a trajectory generated from element $L_k^{{\theta }_{\mathrm{th}}}$ of

$$\begin{array}{c}\hfill L^{{\theta }_{\mathrm{th}}}=\{L_k\mid \Theta (\mathcal{T},L_k)<{\theta }_{\mathrm{th}},\forall k\in K\},\end{array}$$

(5)

obtained using any ${\theta }_{\mathrm{th}}$. Then, we use the function $\Theta$ defined in Equation (1) into Equation (5). The dataset for month m consists of the pair $({\mathit{f}}_k^{{\theta }_{\mathrm{th}}},c)$ of feature ${\mathit{f}}_k^{{\theta }_{\mathrm{th}}}$ generated by $T_k^{{\theta }_{\mathrm{th}}}$, and the corresponding class c is denoted by ${\mathcal{D}}_m^{{\theta }_{\mathrm{th}}}$. If we conduct a sensitivity analysis of the threshold ${\theta }_{\mathrm{th}}$ for ${\mathcal{D}}_m^{{\theta }_{\mathrm{th}}}$, we can verify the effect of trajectory sparsity on the classification accuracy.

In this experiment, a sensitivity analysis is conducted in the range ${\theta }_{\mathrm{th}}\in \{0.95,0.96,\dots ,1.00\}$. Then, ${\mathcal{D}}_m^{{\theta }_{\mathrm{th}}=1.00}$ indicates that the entire trajectory $T_k\forall k\in K$ is to be used for training and testing. When sampling data using any random seeds, we ensure that the datasets become

$$\begin{array}{c}\hfill {\mathcal{D}}_m^{{\theta }_{\mathrm{th}}=0.95}\subseteq {\mathcal{D}}_m^{{\theta }_{\mathrm{th}}=0.96}\subseteq \cdots \subseteq {\mathcal{D}}_m^{{\theta }_{\mathrm{th}}=1.00},\forall m\in M.\end{array}$$

The sizes of the training and test data in the datasets ${\mathcal{D}}_m^{{\theta }_{\mathrm{th}}},\forall m\in M$ for each ${\theta }_{\mathrm{th}}$ value are listed in

Table 1. Sizes of training and test data in datasets ${\mathcal{D}}_m^{{\theta }_{\mathrm{th}}},\forall m\in M$ by ${\theta }_{\mathrm{th}}$ setting.

Data	Month	0.95	0.96	0.97	0.98	0.99	1.00
$S^{\mathrm{train}}$	-	96	198	300	510	1044	2286
$S_m^{\mathrm{test}}$	January	42	84	126	216	444	972
	February	102	180	330	612	1176	2538
	March	834	1338	2004	3222	6342	16,752
	April	1146	1686	2592	4194	7848	21,660
	May	1044	1476	2556	4206	7680	20,952
	June	1638	2280	3996	6810	12,888	38,136

. Additionally, for the models for each month, we adopted the GMD components and the IDNN hyperparameters with the maximum average accuracy in Section 6.1.1 and Section 6.1.2.

Figure 8 summarizes the flow of all the experimental procedures.

6.2. Results and Discussion: Tuning GMD Components n for Model Performance

Based on Section 6.1.1, the trends in the classification performance of each ML model for GMD components n are verified. Figure 9 shows the test accuracy values obtained using components $n\in \left\{2\right\}\cup \{10,20,\dots ,200\}$ and months $\forall m\in M$. The results are shown for the SVC with an RBF kernel, RF, and DNN, which are represented by red, blue, and green lines, respectively. In particular, the averages and standard deviations (1SD) of each model by dataset ${\mathcal{D}}_m$ for the 20 seeds are indicated. The vertical and horizontal axes represent the accuracy and GMD components n, respectively. The title presents the model that achieved the highest average accuracy, along with the GMD component number n, which yielded the highest average accuracy. In addition, we conducted the Shapiro–Wilk test on the distribution of 20 accuracy values when the average was maximal, confirming the normality of the distribution. According to the test results, $p\ge 0.050$, $p<0.050$, and $p<0.010$ are indicated as ‘n.s.’, ‘*’, and ‘**’, respectively. The purpose of the Shapiro–Wilk test is to define the confidence interval for the classification accuracy of the ML models and determine the significance of the ML model for the random six-class classification model. Figure 10 shows the training accuracy values, with the vertical and horizontal axes providing the same information as in Figure 9.

Figure 9 shows that for all months, the accuracy tends to increase as n increases, particularly in the SVC and RF. This is because if the GMD components n are tiny, the expressiveness of the RPs in the target area might be rough, and the difference in the similarity of each ID feature might be less likely to occur. In particular, models with minimal components, such as $n=2$, significantly reduce the classification performance. Conversely, components n that are too large may cause the expressiveness to peak and are less likely to contribute to the classification performance. The classification performance varies by the period of data acquisition; focusing on the maximal average by each month, the highest performance is achieved in January by SVC, with an accuracy of =0.336 (0.013), and the lowest performance is achieved in June by RF, with an accuracy of =0.238 (0.040) and a chance level of approximately $0.170$. Nonetheless, in any month, the normality of the accuracy distribution is confirmed to be $p>0.050$, and the classification performance, excluding the results from DNN, is higher than that of random classification ($\mathrm{accuracy}\approx 0.170$) in the range of 3SD, i.e., a confidence interval of 99.7%. As a result, features consisting of only sparse trajectories would include sufficient information to classify the user age groups.

Here, we focus on the learning trends of each model. The SVC achieves the highest performance in January and February, whereas the RF achieves superior accuracy from March to June. Notably, during the period from March to June, the RF consistently outperforms the other two models across all the components n. In contrast, the DNN achieves significantly lower accuracy than both the SVC and the RF do in all months. This may be because the sample size of training dataset is small for DNN. Therefore, the following discussion focuses primarily on the SVC and RF.

Focusing on the transitions in accuracy with respect to the GMD components n, RF tends to reach accuracy saturation earlier than SVC across all months for the test data. This trend is also observed in the transitions in accuracy for each model on the training data, as shown in Figure 10, although the training accuracy of SVC tends to increase more gradually than that of the test data. From this observation, it can be seen that although SVC tends to have lower accuracy than RF does, it appropriately handles the complexity associated with the increase in the dimensionality of feature vectors, effectively suppressing overfitting. This is particularly evident during the period from March to June. On the other hand, RF demonstrates high accuracy even with small values of n, that is, low-dimensional feature vectors, but it also tends to overfit. From these results, it can be concluded that SVC exhibits higher generalization performance than RF does. Moreover, the RF achieves high accuracy even with low-dimensional feature vectors. In the computation of the GMD, an increase in the number of components n leads to an increase in the computational cost; thus, the RF is appropriate for scenarios where the computational cost is a priority. Conversely, when the generalization performance is emphasized, SVC is preferred. These findings are based on the accuracy obtained within the hyperparameter space described in Section 5.3. Hence, different results may emerge when a broader hyperparameter space is considered.

Note that because the GMD components n obtained in this experiment are the results for only one area, i.e., Narashino, Chiba, Japan, the components are not necessarily practically applicable to any area and could depend on the scale and shape of the target area.

6.3. Results and Discussion: Comparison of the GMM and IDNN

The performance of the ML models that use the feature vectors ${\mathit{f}}_k^{\mathrm{GMM}}$ and ${\mathit{f}}_k^{\mathrm{IDNN}}$ is compared on the basis of Section 6.1.2.

Table 2. Average and standard deviation of 20 accuracies on training and test data in ${\theta }_{\mathrm{th}}=1.00$ by the combinations of feature extraction methods and ML models for each month. The highest average accuracy and its standard deviation are indicated in bold. The parameters of the GMM and IDNN that yielded the highest average accuracy on test data are indicated.

Month	Data	GMM-SVC	GMM-RF	GMM-DNN	IDNN-SVC	IDNN-RF	IDNN-DNN
January	Train	0.631 (0.008)	0.597 (0.011)	0.415 (0.027)	0.907 (0.019)	0.794 (0.017)	0.350 (0.039)
	Test	0.336 (0.013)	0.331 (0.011)	0.310 (0.012)	0.323 (0.014)	0.298 (0.012)	0.268 (0.018)
	Param.	(200)	(60)	(180)	(2, 200)	(2, 200)	(2, 200)
February	Train	0.575 (0.033)	0.599 (0.016)	0.390 (0.030)	0.907 (0.017)	0.785 (0.024)	0.313 (0.033)
	Test	0.298 (0.010)	0.296 (0.008)	0.276 (0.009)	0.296 (0.009)	0.270 (0.009)	0.236 (0.012)
	Param.	(160)	(50)	(200)	(2, 200)	(2, 150)	(2, 200)
March	Train	0.482 (0.065)	0.635 (0.035)	0.310 (0.020)	0.905 (0.058)	0.812 (0.029)	0.285 (0.036)
	Test	0.234 (0.006)	0.243 (0.004)	0.221 (0.007)	0.226 (0.007)	0.224 (0.006)	0.205 (0.011)
	Param.	(120)	(180)	(120)	(2, 200)	(3, 200)	(2, 200)
April	Train	0.453 (0.080)	0.590 (0.035)	0.326 (0.028)	0.712 (0.203)	0.743 (0.059)	0.290 (0.037)
	Test	0.233 (0.005)	0.244 (0.005)	0.227 (0.007)	0.226 (0.004)	0.229 (0.006)	0.216 (0.006)
	Param.	(130)	(180)	(190)	(2, 200)	(2, 150)	(2, 200)
May	Train	0.500 (0.095)	0.599 (0.038)	0.331 (0.025)	0.715 (0.222)	0.756 (0.054)	0.294 (0.033)
	Test	0.235 (0.004)	0.250 (0.004)	0.230 (0.004)	0.228 (0.005)	0.232 (0.004)	0.220 (0.006)
	Param.	(180)	(150)	(190)	(2, 200)	(2, 200)	(2, 200)
June	Train	0.443 (0.102)	0.586 (0.062)	0.316 (0.027)	0.578 (0.195)	0.736 (0.107)	0.285 (0.040)
	Test	0.224 (0.004)	0.238 (0.004)	0.221 (0.006)	0.220 (0.005)	0.223 (0.004)	0.211 (0.010)
	Param.	(160)	(80)	(160)	(2, 200)	(2, 200)	(2, 200)

shows the average and standard deviation of 20 accuracies on training and test data by month in ${\theta }_{\mathrm{th}}=1.00$. The header of the table indicates the combination of the feature extraction method and the ML model, e.g., GMM-SVC indicates that SVC trained by ${\mathit{f}}_k^{\mathrm{GMM}}$. In addition, the parameters providing the highest average accuracies for each feature extraction method on the test data are indicated in parentheses. These include $\left(n\right)$, representing the number of GMD components for the GMM, and $(\lambda ,n^{\left(l\right)})$ for the IDNN, which denotes the number of hidden layers and the number of neurons in each layer, respectively. These accuracies are obtained via the SVC, RF, and DNN with optimized hyperparameters in space described in Section 5.3. The highest average accuracy for each month is indicated in bold.

The difference in accuracy between the training and test data serves as a critical indicator of the model’s generalization performance. Specifically, when training accuracy is substantially higher than test accuracy, it indicates that the model is overfitting. Conversely, a small difference in accuracy implies that the model exhibits high generalization performance. Therefore, evaluating the accuracy difference between training and test data provides a more detailed assessment of each model’s generalization capabilities.

The results show that the accuracy of the GMM-SVC is the highest in January and February, whereas the GMM-RF achieves the highest accuracy from March to June. GMM-SVC and GMM-RF consistently show higher accuracy than IDNN-SVC and IDNN-RF across all months. The scores of IDNN-SVC fall within the 1–3 SD range of the scores for GMM-SVC and GMM-RF, indicating relatively small differences in accuracy. On the other hand, focusing on the accuracy differences between the training and test data, the range for GMM-SVC and GMM-RF is 0.019–0.392, whereas IDNN-SVC and IDNN-RF exhibit a wider range of 0.358–0.679. This suggests that IDNN-SVC and IDNN-RF are prone to overfitting, whereas feature extraction via the GMM demonstrates better generalization performance. Note that, except for IDNN-SVC in March and May and IDNN-DNN in March and June, all the results have $p>0.050$, confirming that the accuracy distributions on test follow a normal distribution.

This experiment was conducted with a missing ratio of ${\theta }_{\mathrm{th}}=1.00$, resulting in highly sparse trajectories. Consequently, the results indicate that the GMM-based features achieve higher accuracy and generalization performance than the IDNN does for sparse trajectories. The reliance of the IDNN on image-based trajectory transformations likely makes it more susceptible to noise, reducing its generalization performance. In other words, imputed trajectory information, when directly converted into pixel values, can lead to irregular spatial patterns in the image, causing noise accumulation. For sparse trajectories, the lack of structural information on the IDNN may limit its effectiveness. By contrast, GMM-based feature generation identifies RPs on the basis of the distribution of human flows and aggregates trajectories around those RPs. This approach may minimize the accumulation of noise caused by sparsity. As a result, GMM-based features demonstrate superior generalization performance for sparse trajectories.

6.4. Results and Discussions: Performance by the Missing Value Rate of the User Trajectory

On the basis of Section 6.4, we analyze the trend of the ML models using the feature vectors ${\mathit{f}}_k^{\mathrm{GMM}}$ and ${\mathit{f}}_k^{\mathrm{IDNN}}$ performance with respect to trajectory sparsity. First, to consider the effect of trajectory sparsity on classification performance, we visualize the trends of the accuracy values for the test and training data against the missing value ratio ${\theta }_{\mathrm{th}}$ in Figure 11 and Figure 12, respectively. For both figures, the vertical and horizontal axes represent the accuracy $A_{\mathrm{M}}$ and missing value ratio ${\theta }_{\mathrm{th}}$, respectively. In addition, the red, blue and green symbols represent the SVC, RF, and DNN, respectively, and the circle and cross symbols indicate the GMM- and IDNN-based features, respectively. Fewer differences between the accuracy of training and test data indicate higher generalization performance. In Figure 11, the horizontal and vertical dashed lines represent the best accuracy and the corresponding ${\theta }_{\mathrm{th}}$ among the combinations of feature extraction methods and ML models, respectively.

In the SVC and RF methods that use GMM and IDNN features, given that the trends in which a lower ${\theta }_{\mathrm{th}}$ leads to higher accuracy were confirmed between the threshold of the best accuracy, i.e., the dashed line in Figure 11 and ${\theta }_{\mathrm{th}}=1.00$, the sparsity of the trajectory would decrease the model performance. Conversely, if ${\theta }_{\mathrm{th}}$ is too small, low performance can also occur. This is because an extremely small ${\theta }_{\mathrm{th}}$ yields low numbers of samples relative to the feature space. Conversely, for both the GMM and IDNN features, the DNN tends to decrease in accuracy as ${\theta }_{\mathrm{th}}$ decreases in the range of ${\theta }_{\mathrm{th}}=0.95$–$0.98$, which contrasts with the trend observed in SVC and RF, where accuracy increases as ${\theta }_{\mathrm{th}}$ decreases. This suggests that, in this range of missing value ratios, the amount of training data may be insufficient to adequately train the DNN models. This result is supported by the findings shown in Figure 12 for the DNN on the training data. Because the trends in accuracy for the training and test data are similar, it can be inferred that insufficient data availability prevents the DNN from learning effectively, thereby hindering its performance improvement.

Focusing on the best accuracy values in Figure 11, GMM-RF achieves the highest accuracy for all months. In particular, for April to June, the test accuracy of GMM-RF is higher than that of the other combinations for all the missing value ratios. In addition, the test accuracy of GMM-SVC is not as high as that of GMM-RF, but it is confirmed that GMM-SVC has high generalization performance because of the smaller differences between the accuracies of the training and test data. This is implied from Figure 11 and Figure 12. On the other hand, while the test accuracies of the GMM and IDNN features are close, the training accuracy of the IDNN features tends to be greater than that of the GMM features, particularly for ${\theta }_{\mathrm{th}}=0.98$–$1.00$. This implies the IDNN features are prone to overfitting in the instance of this experiments. This trend was confirmed in Section 6.3, and its reasons are the same as those mentioned before; that is, the conversion process of a trajectory to an image cannot eliminate noise. These results also suggest that the GMM-based features are more effective for handling sparse trajectories than IDNN-based features. Overall, the results indicate that GMM-based features provide a more reliable and effective approach for extracting information from sparse trajectories, making them a better fit for scenarios where semantic information is unavailable and data sparsity is a significant challenge. For ${\theta }_{\mathrm{th}}=0.95$–$0.97$, overfitting is confirmed in GMM-SVC and RF and IDNN-SVC and RF. This may be because the number of samples is small in the feature space, and the model is overfitted to the training data.

Here, we analyze the classification results via detailed metrics, such as precision, recall, and F-measure, in addition to accuracy. Note that we focus on only the GMM-RF results, which showed the highest accuracy in all months.

Table 3. Four metric values for the GMM-RF model with an n-dimensional feature vector. This table shows the metric values for each month and each ${\theta }_{\mathrm{th}}$. The highest average accuracy and its standard deviation are indicated in bold. Moreover, the p-value for the Shapiro–Wilk test is used to confirm the distribution normality of 20 accuracy samples.

Month	n	${\mathit{\theta }}_{\mathbf{th}}(>{\mathit{\theta }}_{\mathit{k}})$	${\mathit{A}}_{\mathbf{M}}$	${\mathit{P}}_{\mathbf{M}}$	${\mathit{R}}_{\mathbf{M}}$	$\mathit{F}{1}_{\mathbf{M}}$	Shapiro–Wilk †
January	60	0.95	0.535 (0.087)	0.562 (0.085)	0.535 (0.087)	0.534 (0.086)	n.s.
		0.96	0.545 (0.064)	0.564 (0.062)	0.545 (0.064)	0.547 (0.063)	n.s.
		0.97	0.528 (0.046)	0.543 (0.044)	0.528 (0.046)	0.529 (0.045)	n.s.
		0.98	0.485 (0.036)	0.496 (0.031)	0.485 (0.036)	0.484 (0.034)	n.s.
		0.99	0.421 (0.021)	0.426 (0.022)	0.421 (0.021)	0.418 (0.022)	n.s.
		1.00	0.331 (0.011)	0.330 (0.013)	0.331 (0.011)	0.320 (0.009)	n.s.
February	50	0.95	0.443 (0.040)	0.464 (0.041)	0.443 (0.040)	0.443 (0.040)	n.s.
		0.96	0.476 (0.038)	0.483 (0.042)	0.476 (0.038)	0.473 (0.039)	n.s.
		0.97	0.428 (0.024)	0.434 (0.027)	0.428 (0.024)	0.426 (0.025)	n.s.
		0.98	0.423 (0.018)	0.426 (0.019)	0.423 (0.018)	0.421 (0.018)	n.s.
		0.99	0.371 (0.016)	0.364 (0.018)	0.371 (0.016)	0.363 (0.017)	n.s.
		1.00	0.296 (0.008)	0.297 (0.009)	0.296 (0.008)	0.292 (0.009)	n.s.
March	180	0.95	0.326 (0.024)	0.341 (0.024)	0.326 (0.024)	0.327 (0.023)	n.s.
		0.96	0.342 (0.017)	0.343 (0.020)	0.342 (0.017)	0.338 (0.019)	n.s.
		0.97	0.339 (0.017)	0.339 (0.017)	0.339 (0.017)	0.334 (0.018)	n.s.
		0.98	0.329 (0.012)	0.330 (0.012)	0.329 (0.012)	0.325 (0.012)	n.s.
		0.99	0.292 (0.009)	0.289 (0.008)	0.292 (0.009)	0.287 (0.009)	n.s.
		1.00	0.243 (0.004)	0.245 (0.006)	0.243 (0.004)	0.239 (0.005)	n.s.
April	180	0.95	0.319 (0.023)	0.328 (0.026)	0.319 (0.023)	0.317 (0.024)	n.s.
		0.96	0.335 (0.023)	0.335 (0.025)	0.335 (0.023)	0.332 (0.024)	n.s.
		0.97	0.325 (0.014)	0.325 (0.017)	0.325 (0.014)	0.321 (0.017)	n.s.
		0.98	0.319 (0.012)	0.318 (0.017)	0.319 (0.012)	0.315 (0.014)	n.s.
		0.99	0.292 (0.008)	0.289 (0.011)	0.292 (0.008)	0.283 (0.010)	n.s.
		1.00	0.244 (0.005)	0.241 (0.006)	0.244 (0.005)	0.232 (0.006)	n.s.
May	150	0.95	0.315 (0.013)	0.317 (0.015)	0.315 (0.013)	0.312 (0.014)	**
		0.96	0.331 (0.015)	0.328 (0.016)	0.331 (0.015)	0.327 (0.015)	n.s.
		0.97	0.325 (0.012)	0.323 (0.014)	0.325 (0.012)	0.321 (0.013)	n.s.
		0.98	0.325 (0.014)	0.325 (0.017)	0.325 (0.014)	0.320 (0.016)	*
		0.99	0.298 (0.009)	0.292 (0.010)	0.298 (0.009)	0.288 (0.011)	n.s.
		1.00	0.250 (0.004)	0.245 (0.005)	0.250 (0.004)	0.235 (0.006)	n.s.
June	80	0.95	0.281 (0.021)	0.288 (0.025)	0.281 (0.021)	0.280 (0.022)	**
		0.96	0.310 (0.011)	0.309 (0.013)	0.310 (0.011)	0.307 (0.012)	n.s.
		0.97	0.316 (0.011)	0.311 (0.012)	0.316 (0.011)	0.311 (0.011)	n.s.
		0.98	0.308 (0.010)	0.306 (0.010)	0.308 (0.010)	0.304 (0.010)	n.s.
		0.99	0.278 (0.007)	0.274 (0.009)	0.278 (0.007)	0.270 (0.007)	n.s.
		1.00	0.238 (0.004)	0.232 (0.005)	0.238 (0.004)	0.227 (0.005)	n.s.

^† n.s.: $p\ge 0.050$, *: $p<0.050$, **: $p<0.010$.

shows each metric for the classification results of the test data by threshold ${\theta }_{\mathrm{th}}$ and month. The averages and standard deviations corresponding to datasets of 20 seeds for each metric are indicated in parentheses. Moreover, with respect to the confidence interval of classification accuracy, as in Section 6.2, Table 3 shows the p-values for the Shapiro–Wilk test for the distribution of accuracies for 20 seeds. Because the p-values are $p>0.050$ except for ${\theta }_{\mathrm{th}}=0.95$ and $0.98$ in May and ${\theta }_{\mathrm{th}}=0.95$ in June, we assume that any distribution of accuracies follows a normal distribution.

Focusing on the precision, recall, and F-measure in any parameter setting, we can conclude that the levels of classification performance for the classes are similar because these metric values are almost equivalent to the accuracy. This is a natural result because we handle a balanced dataset. For the missing value rate ${\theta }_{\mathrm{th}}=1.00$, because the average accuracy is at its maximum value of 0.331 (0.011) in January and at its minimum value of 0.238 (0.004) in June, it is confirmed that the classification performance is significantly higher than random classification in any month. Furthermore, through adjustments in the threshold ${\theta }_{\mathrm{th}}$, the classification performance is increased such that the average accuracy reaches its maximum of 0.545 (0.064) in January and minimum of 0.316 (0.011) in June. The normality of the accuracy distribution for the 20 seeds is confirmed, except for ${\theta }_{\mathrm{th}}=0.95$ and $0.98$ in May and ${\theta }_{\mathrm{th}}=0.95$ in June. These results show that with appropriate parameter settings, the constructed GMM-RF models perform with significantly higher classification accuracy than those produced by random classification. Note that, in GMM-SVC, the classification performance is not as high as that of GMM-RF, but it is confirmed that GMM-SVC has the same tendency and high generalization performance, as shown in Figure 11 and Figure 12. As a result, we infer that sparse spatial–temporal information, i.e., a sparse trajectory, would have sufficient information to classify the corresponding user age group.

Finally, the classification performance with respect to the data acquisition period is discussed. Focusing on the maximum accuracy of GMM-RF, the accuracy values from January to June are 0.545 (0.064), 0.476 (0.038), 0.342 (0.017), 0.335 (0.023), 0.331 (0.015), and 0.316 (0.011), respectively. In January and February, the values are particularly high. Conversely, the values from March to June are lower than those from January and February. This occurred because declarations of a state of emergency and a semistate of emergency were enforced from 8 January to 21 March and 20 April to 1 August 2021, respectively, in the target area. The state of emergency involves strong behavioral regulation, including orders to close restaurants and commercial facilities. In contrast, the semistate of an emergency is a behavioral regulation that includes only the required orders regarding hours of operation [45]. Thus, in January and February, in which states of emergency were enforced, the restriction of irregular behavior (e.g., amusement trips to commercial and recreational facilities) among people’s activities clearly differentiated the behavior patterns of people in different user age groups. However, in March, when a state of emergency was also enforced, the classification performance was significantly lower than that in January and February. This could be because the states of emergency migrated to semistates of emergency in the middle of March, and people’s behavior changed remarkably. The classification performance also declined from April to June; this may be because people’s behavior became more diverse than that during the state of emergency. The monthly sample sizes shown in Figure 5 and Figure 6 indicate that the sample size from March to June is larger than that from January to February. This suggests that people’s behavior patterns gradually returned to normal during the period from March to June. Nonetheless, the average accuracy values are significantly higher than those of random classification in the range of 3SD for many parameter settings. These results also reveal that the sparse trajectory contains sufficient information for classifying user age groups without the use of any semantic information other than spatial–temporal information.

7. Conclusions

In this study, we conducted a performance evaluation and analysis of GMM-based features for sparse trajectories without semantic information. The experiments utilized three machine learning models, namely, SVC, RF, and DNN, to perform a detailed examination of the features. We also compared them with the IDNN, an existing feature extraction method for trajectories. The results showed that GMM-based features achieved higher classification accuracy and better generalization performance than IDNN-based features. In particular, the RF algorithm attained high accuracy, whereas the SVC algorithm maintained stable generalization performance. Moreover, a sensitivity analysis of the missing rate ${\theta }_{\mathrm{th}}$, which indicates trajectory sparsity, revealed that the IDNN becomes more prone to overfitting as the sparsity increases, whereas the GMM-based approach can sustain accuracy and generalizability performance. These contributions demonstrate that GMM-based features offer clear advantages over existing methods in terms of both accuracy and generalizability when dealing with sparse trajectories lacking semantic information. Consequently, this study shows that sparse trajectory data can be effectively utilized even without semantic information. This study expands the potential for data utilization in various fields, including urban planning, transportation analysis, and public policy.

The limitations of this study are as follows:

• Although simple linear interpolation was adopted for trajectory imputation, its impact on classification performance was not verified.
• The proposed method was verified only for a single region. In addition, owing to constraints in the data acquisition period, the influence of behavioral changes from prepandemic norms during the COVID-19 era could not be taken into account.
• The features did not contain the connection relationships of RPs.

In light of these points, our future research will verify the generality of the proposed model in terms of the area and period of location data and further consider the utilization of sparse spatial–temporal information.