# AutoSoccerPose: Automated 3D posture Analysis of Soccer Shot Movements

Calvin Yeung  Kenjiro Ide
Nagoya University, Japan
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>Keisuke Fujii
Nagoya University, Japan
RIKEN/JST PRESTO, Japan
<EMAIL_ADDRESS>
###### Abstract

Image understanding is a foundational task in computer vision, with recent
applications emerging in soccer posture analysis. However, existing publicly
available datasets lack comprehensive information, notably in the form of
posture sequences and 2D pose annotations. Moreover, current analysis models
often rely on interpretable linear models (e.g., PCA and regression), limiting
their capacity to capture non-linear spatiotemporal relationships in complex
and diverse scenarios. To address these gaps, we introduce the 3D Shot Posture
(3DSP) dataset in soccer broadcast videos, which represents the most extensive
sports image dataset with 2D pose annotations to our knowledge. Additionally,
we present the 3DSP-GRAE (Graph Recurrent AutoEncoder) model, a non-linear
approach for embedding pose sequences. Furthermore, we propose AutoSoccerPose,
a pipeline aimed at semi-automating 2D and 3D pose estimation and posture
analysis. While achieving full automation proved challenging, we provide a
foundational baseline, extending its utility beyond the scope of annotated
data. We validate AutoSoccerPose on SoccerNet and 3DSP datasets, and present
posture analysis results based on 3DSP. The dataset, code, and models are
available at: https://github.com/calvinyeungck/3D-Shot-Posture-Dataset.

##  1 Introduction

Refining the interpretation of images is a fundamental aim in computer vision
[26]. It encompasses tasks such as image classification [26], object detection
[38], and human pose estimation [23]. In sports, image interpretation has
significantly contributed to foul detection [43], player identification [42],
and tracking [39]. Recently, such computer vision techniques have been applied
to analyze soccer posture effectively [49].

In soccer posture analysis, a high-quality dataset and the analysis method
hold profound effects on the analysis result. Concerning the pose dataset, for
soccer pose and other human movements, the pose sequence is often captured in
2D given the use of a standard RGB camera [30]. Currently, publicly available
2D pose datasets include MPII [2], OCHuman [55], and COCO-WholeBody [23]. For
soccer and sports image datasets, options include LSP [24], Sports-102 [4],
LearningFromThePros [49], and Sport Image [48]. These datasets have been
effectively utilized for their designated tasks. Regarding the analysis
method, existing studies have analyzed the goalkeeper [49, 36] and shooter
[33, 19] postures as they directly affect the match result. With a simplified
scenario and linear model, studies have provided valuable insights into
posture dynamics.

However, concerning the dataset, the large 2D pose datasets often lack a
sufficient number of soccer images, and the soccer-specific datasets provide
single images of postures instead of sequences, failing to represent entire
movements. Regarding the analysis method, the simplified scenarios may not
fully capture the complexities of real-world situations, limiting the
applicability of results, especially at the highest competitive levels.
Additionally, the utilization of linear models may overlook more intricate
relationships in the data.

Figure 1: Overview of AutoSoccerPose.

Therefore, in this study, we propose the 3D Shot Posture (3DSP) dataset,
consisting of annotated 2D pose sequences of shooting instances from
professional soccer matches. We introduce the 3DSP-GRAE (Graph Recurrent
AutoEncoder) model, a non-linear model designed to encode pose sequences into
vectors, and AutoSoccerPose, a semi-automated pipeline for extracting 2D and
3D pose sequences directly from professional soccer broadcast videos for
analysis, as illustrated in Fig. 1. The contributions of this paper are as
follows:

  * •

3DSP (3D Shot Posture) dataset, which comprises 2D pose annotations, 3D pose
estimations, and tracklet information extracted from professional broadcast
videos.

  * •

3DSP-GRAE model, a novel nonlinear spatiotemporal model designed to encode
pose sequences into vector representations effectively.

  * •

AutoSoccerPose, a pipeline tailored for semi-automated 3D posture analysis in
soccer, specifically focusing on player shot movements.

  * •

Extensive evaluations at each stage of AutoSoccerPose, accompanied by an
analysis of shooting postures utilizing the 3DSP dataset.

##  2 Related work

###  2.1 2D pose image datasets

For images consisting of humans, 2D human pose estimation has emerged as a
well-established research topic. This task aims to accurately define human
body part segments or joint keypoints in the image. Commonly used benchmarks
for this task include OCHuman [55], COCO-WholeBody [23], both derived from
COCO [26], and MPII [2] dataset, covering various scenarios such as daily
activities, human interactions, and sports events.

In the realm of sports, image frequently revolves around the posture of
athletes. While traditional computer vision tasks such as image classification
[4], 2D/3D pose estimation [24, 20], and its extension, human parsing [48],
can be undertaken, the intricate anatomy involved, coupled with challenging
image conditions like motion blur [24, 4], amplify the difficulty of these
tasks. In addition, advancements in computer vision techniques and soccer
analytics, particularly those concerning player action value [11, 53, 54],
have opened avenues for leveraging soccer player posture in computer vision-
based player evaluation [49].

Nonetheless, while multiple datasets consist of pose images of professional
soccer athletes, there are a few limitations when applying them for posture
analysis:

  1. 1.

General-purpose dataset: Despite the potential abundance of images in
conventional datasets, the proportion of professional soccer images is low due
to the multitude of classes they encompass. Furthermore, these images do not
specifically cater to distinct soccer actions, posing challenges for
comparison.

  2. 2.

Partial information: Often, a complete soccer action, such as shooting, is
captured through multiple sequential images (poses) [19, 33]. However, current
datasets only offer a single pose image for each discrete action. Moreover, 2D
pose data relies on computer vision-based model estimation, introducing errors
given the unique characteristics of sports images.

To mitigate these issues, we have collected the 3DSP dataset (see Sec. 4),
currently the largest repository of professional soccer human pose data with
2D posture annotations. This dataset focuses on the shooting action, providing
sequential images of the shooting movement. In Tab. 1, we summarize the
existing datasets containing professional soccer pose images along with the
3DPS dataset.

Dataset |  | 2D Pose
---
Annotation
| Sequential
---
Images
| Images in
---
Soccer
| Number of
---
Images
| Soccer
---
Action
Purpose
LSP [24] | ✓ | ✗ |  | 2k |  | 2D Pose Estimation
Sport Image [48] | ✓ | ✗ |  | 1.3k |  | Human Parsing
MPII [2] | ✓ | ✗ | 0.2k | 41k |  | 2D Pose Estimation
OCHuman [55] | ✓ | ✗ |  | 4.7K |  | 2D pose estimation
COCO-WholeBody [23] | ✓ | ✗ |  | 200K |  | 2D pose estimation
LearningFromThePros [49] | ✗ | ✗ | 1k | 1k | GoalKeeping | Posture Analysis
Sports-102 [4] | ✗ | ✗ | 0.2k | 14k |  | Image Classification
3DSP (Ours) | ✓ | ✓ | 4k | 4k | Shooting | Posture Analysis
Table 1: Overview of human pose in sports (professional soccer included)
dataset. The datasets are arranged in order of their publication years. Due to
the challenges associated with collecting 3D pose annotations during
professional football matches, there currently exists no publicly available
dataset providing 3D pose annotations for professional football.

###  2.2 Posture analysis in soccer

Posture analysis holds profound importance in soccer, given its direct impact
on player performance and strategic optimization. Previous studies have
predominantly focused on shooting and goalkeeping postures, recognizing their
pivotal roles in influencing match outcomes.

In shooting analysis, researchers have delved into unraveling the intricate
differences between various kicking techniques. For instance, [33] examined
disparities between non-rotational shots and instep kicks, while [19]
investigated distinctions among straight shots, knuckling shots, and curve
shots. Both studies utilized advanced camera setups and reflective markers to
capture 2D poses of college soccer players’ shooting postures. The
differentiation was based on 2D pose data, kinematic parameters, and principal
component analysis techniques.

On the other hand, studies on goalkeeping posture, such as [49], have focused
on analyzing goalkeeper saving techniques in penalty and 1-vs-1 situations.
Utilizing techniques like estimated 3D pose extracted with PoseHG3D [58],
engineered features, and K-means clustering, researchers examined the various
postures used by goalkeepers. Similarly, [36, 35] examined the strategies
employed by goalkeepers in penalty kicks. They employed logistic regression to
infer the relationship between orientation and strategies, where the
orientation was engineered from the 2D pose of goalkeepers and penalty takers,
utilizing OpenPose [6] for pose estimation.

While previous studies provide valuable insight into the posture of soccer
players, several drawbacks could limit the development of posture analysis:

  1. 1.

Simplified scenario: Previous datasets are mainly collected under simplified
scenarios, such as competitive settings with few players or controlled
experimental environments. However, soccer actions often occur in occluded and
crowded situations captured in broadcast videos.

  2. 2.

Labor-intensive: Collecting data for posture analysis involves labor-intensive
processes including image preprocessing, pose estimation or annotation, and
quality control. Even for the datasets used in this study, the collection
process had taken months to complete.

  3. 3.

Linear model: Previous methods like principal component analysis and logistic
regression offer valuable insights but may fail to capture non-linear
spatiotemporal relationships present in posture data.

Given these limitations, this study proposes the AutoSoccerPose for semi-
automated 3D posture estimation from professional soccer broadcast video (see
Sec. 3) and a non-linear graph-based sequential model 3DSP-GRAE for posture
analysis (see Sec. 3.5).

##  3 Proposed method

AutoSoccerPose111For detailed specifications of the AutoSoccerPose, please
refer to https://github.com/calvinyeungck/3D-Shot-Posture-Dataset. aims to
extract 3D poses from broadcast videos and analyze the extracted posture. Each
stage of AutoSoccerPose corresponds to an existing computer vision task in
soccer. Fig. 1 provides an overview of the stages involved in AutoSoccerPose,
and each stage is elaborated upon in the subsequent sections.

###  3.1 Broadcast video

Similar to all posture analyses, Clips that capture the player’s posture are
fundamental to AutoSoccerPose. For broadcast video, there exist publicly
available professional matches videos that could be collected, as in [36], or
existing publicly available datasets like SoccerNet [12, 17] that provide a
large number of videos (nearly 800 hours of professional soccer matches).
Nonetheless, while the posture of professional players in broadcast videos is
the most informative and commonly available data, other less complex video
footage like those in previous studies (see Sec. 2.2) could be utilized in
AutoSoccerPose.

Besides, accurate timestamp annotation of the action, shooting in this case,
is critical to cutting clips from the broadcast video and the performance of
AutoSoccerPose (see Sec. 5.2). The timestamp of important events like shooting
could easily be retrieved in the post-match report, manually annotated, or
publicly available datasets [12, 17, 34] as it was a fundamental type of
soccer data annotation. Furthermore, accurate timestamp annotation has been a
growing soccer computer vision task in recent times, in SoccerNet Challenge
2023 [10] to facilitate broadcast video understanding, the task action
spotting (locate when and what type of event happened in the broadcast video)
has been proposed. We reserve the integration of action spotting for future
research.

To be more specific, the AutoSoccerPose was currently developed with the
videos and timestamp annotation in SoccerNet [12, 17], with frames resolution
of $1280\times 720$ and 25 fps. The retrieved clips used for AutoSoccerPose
cover 0.5 seconds before and after the annotated timestamps. However via
observation, the first 20 frames cover the entire shot movement for most of
the shots, therefore only the first 20 were utilized for each shot.

###  3.2 Tracking

Tracking in AutoSoccerPose aims to identify soccer players in video clips and
create tracklets: pixel coordinates of each respective player within each
frame of the clips. Conventionally, two main tracking methods are utilized:
tracking-by-detection [50, 5, 13], where objects are initially identified with
bounding boxes in each frame and then associated (tracked) based on various
conditions; and end-to-end tracking [57, 32, 51], where tracklets (bounding
boxes of an object) are generated directly from the model.

State-of-the-art soccer player tracking models, achieving the best performance
in the SoccerNet Challenge 2023, predominantly follow the tracking-by-
detection approach. However, these tracking models are often not open-sourced.
Here, we adopt the tracking by detection method using widely used models:
fintuned YOLOv8222YOLOv8: https://docs.ultralytics.com/. as the detector and
BoTSort [1] as the tracking algorithm. For performance evaluation and
discussion of YOLOv8-BoTSort, refer to Sec. 5.1.

###  3.3 Tracklet selection

Tracklet selection aims to pinpoint the tracklet associated with the shooter.
This task can be viewed as a facet of image comprehension, where the objective
is to classify whether a sequence of images pertains to the shooter. Analogous
image comprehension tasks in soccer encompass sport type classification [4],
foul detection [18, 14], and action spotting [10]. Nonetheless, many existing
studies propose intricate deep-learning architectures that necessitate
substantial and diverse datasets for training or finetuning.

To avoid the need for an extensive dataset, we opt for a simpler approach
employing a Convolutional Neural Network (CNN) model1. Our model processes the
sequence of images through three standard convolutional blocks, each
comprising a convolutional layer, ReLU activation function, max pooling layer,
batch normalization, and dropout layer. Subsequently, the feature maps of the
images within the sequence are concatenated. Finally, the concatenated feature
maps are passed through a three-layer Multi-Layer Perceptron (MLP) with ReLU
activation to predict whether the tracklet contains the shooter.

For each tracklet, we extract a sequence of images spanning frames 10 to 15.
From each frame, we crop a tracklet image of size 96 $\times$ 96 pixels, with
the center of the bounding box of the detected player (as described in Sec.
3.2) serving as the center of the tracklet image. The shooter tracklet for
each clip is determined by selecting the tracklet with the highest estimated
probability from the CNN model. Definitions of clips and frames can be found
in Sec. 3.1, while details regarding the selection of necessary frames are
outlined in Sec. 5.2.

###  3.4 2D & 3D pose estimation

The objective of 2D and 3D pose estimation is to derive posture details from
video clips, with a specific emphasis on capturing the movements of the
shooter. 3D pose estimation methods can be categorized into direct 3D pose
estimation [8, 21, 37] and 2D-3D lifting [30, 52, 41], which involves
converting 2D pose coordinates into 3D space. The MPI-INF-3DHP dataset [31] is
a widely used benchmark for in-the-wild 3D human pose estimation, resembling
scenarios found in soccer broadcasts. Among the best-performing models, those
employing the 2D-3D lifting approach, such as MotionAGFormer [30], are
prevalent and thus adopted for 3D pose estimation in AutoSoccerPose.

For 2D pose estimation, methodologies are typically classified into bottom-up
approaches [6, 7, 16], which detect points in the image and then associate
keypoints to form 2D poses, and top-down approaches [22, 41, 52], which first
detect individuals using object detection models and then perform pose
estimation within cropped bounding boxes. Top-down approaches, often more
efficient, are favored, with RTMPose [22] selected for AutoSoccerPose due to
its superior performance on datasets like COCO-whole body [23] and our
dataset, 3DSP.

In implementation, for shooter tracklets, we extract sequences of 20 images,
each sized 100 $\times$ 100 pixels, covering the detected players’ bounding
boxes and entire bodies, akin to Sec. 3.3. However, cropping images larger
than the bounding box might include other players besides the shooter. To
mitigate this, as the shooter is always centered in the cropped image, we
select the shooter’s 2D pose as the one with the center of body (torso) joint
closest to the image center. The 2D pose from RTMPose [22] is then lifted to
3D using MotionAGFormer [30]. The performance of both 2D and 3D pose
estimation is further discussed in Sec. 5.3.

###  3.5 Posture analysis

Posture analysis seeks to distinguish various shooting styles among soccer
players. Given that shot movements are typically unlabeled, unsupervised
learning represents the most suitable approach for this task. Furthermore, as
outlined in Sec. 2.2, prior research [19, 36] predominantly employs linear
models such as Principal Component Analysis (PCA) and logistic regression to
infer 3D posture, potentially limiting their ability to capture non-linear
spatiotemporal relationships. Hence, to conduct unsupervised learning on
posture data with a model capable of capturing the non-linear spatiotemporal
properties, we propose the 3DSP-GRAE (Graph Recurrent AutoEncoder) model1 to
encode the latent (representation) vector, illustrated in Fig. 2. This model
draws inspiration from the Graph Convolutional Network (GCN) [25] and LSTM
AutoEncoder (LSTM-AE) [40]. While several graph and sequential-based models
exist [30, 47, 9], none are explicitly designed for 3D pose encoding. Thus,
leveraging the relatively simple structure of 3DSP-GRAE offers the advantage
of requiring less data for training compared to fine-tuning existing models.

Figure 2: Overview of the 3DSP-GRAE model. The variables $X$, $H$, $C$, and
$t$ in the LSTM layer denote the input features, hidden state, cell state, and
time step, respectively.

In general, the 3DSP-GRAE model utilizes an encoder-decoder architecture. For
the encoder, 3D poses pass through a GCN to extract spatial relationships.
Subsequently, the processed 3D poses are sequentially fed into an LSTM to
extract nonlinear and temporal relationships. The output of the LSTM at the
last layer and the last hidden state serve as the latent variables (vector
representation) of the respective shot movement. In the decoder, another LSTM
takes the latent variables and the processed 3D poses with GCN as input, but
in reverse order. The LSTM output is then converted back to coordinates using
another GCN. Consequently, the decoder output comprises the predicted 3D
poses. Lastly, by employing the mean square error between the predicted 3D
pose and the original input, the 3DSP-GRAE model can be trained.

In detail, each pose comprises 17 joints, as outlined in [30], with each joint
represented by XYZ coordinates. The GCN consists of 2 simplified graph
convolutional layers with ReLU activation, and the adjacency matrix represents
connections between the joints with self-connections. Compared to the original
graph convolutional layer [25], the parameter for the adjacency matrix is
ignored, considering it to be the identity matrix. Furthermore, the encoder
LSTM takes the pose input from pose 1 to 20 in temporal order. Meanwhile, in
the decoder LSTM, the first input is zero-padded, followed by poses 20 to 2.

##  4 3DSP dataset

In this section, we introduce the 3D Shot Posture (3DSP) dataset3333DSP
dataset available at https://github.com/calvinyeungck/3D-Shot-Posture-
Dataset., a crucial resource for developing AutoSoccerPose, serving as a
benchmark for 2D pose estimation in professional soccer match broadcast
videos. Furthermore, as evidenced by the comparative study outlined in Tab. 1,
3DSP comprises sequential images capturing soccer movements previously
unavailable, with a specific focus on distinct soccer actions such as
shooting. 3DSP also represents the largest collection of soccer images
annotated with 2D pose information. Additionally, the 3DSP dataset facilitates
advancements in soccer posture analysis and offers direct applicability for
future studies in this domain. Moreover, the 3DSP dataset could be extended
easily with AutoSoccerPose.

Data collection process: 3DSP mirrors the process of AutoSoccerPose (with more
details and code provided on GitHub 1), albeit with manual refinement and
annotation to ensure the highest quality of 2D poses. Initially, broadcast
videos were sourced from SoccerNet [12, 17], and utilizing SoccerNet
annotations on actions, the videos were trimmed (0.5 seconds before and after
the annotated timestamp, totaling 25 frames). Subsequently, tracklets for
these clips were generated employing a fine-tuned YOLO v82 in conjunction with
BoT-Sort [1]. Furthermore, the tracklet corresponding to the shooter was
manually selected and refined, and a cropped image was produced using the
bounding box. The 2D poses of the first 20 frames (determined empirically)
were then manually annotated and lifted to 3D utilizing the MotionAGFormer
[30].

To summarize the dataset structure, we collected data on shot movements from
22 unique matches in the 2015-2016 English Premier League. The training set
comprises 20 cropped images (100 $\times$ 100 pixels) per shot movement, with
a total of 200 shot movements, resulting in 4000 images. Annotations encompass
tracklets, 2D poses, and estimated 3D poses. In the test set, there are also
20 cropped images per shot movement, with a total of 10 shot movements;
however, only tracklet information is available. For both sets of data,
references to SoccerNet [12, 17] are provided, along with the code utilized
for extending 3DSP with SoccerNet [12, 17].

##  5 Experiments

In this section, we verify each step of the AutoSoccerPose (see Sec. 3). More
details is available at GitHub1.

###  5.1 Tracking performance on SoccerNet [12]

In this subsection, we assess the performance of AutoSoccerPose tracking using
the test set provided in the SoccerNet Player Tracking Challenge 2023 [10].
This dataset comprises ground truth tracking annotations derived from
broadcast videos in SoccerNet [12, 17]. We conduct evaluations for both player
detection and tracking.

For player detection, we employ evaluation metrics including precision,
recall, and average precision (AP) (see Tab. 2 caption). The evaluated models
encompass state-of-the-art object detection models introduced subsequent to
the SoccerNet Player Tracking Challenge 2023 [10], including YOLOv8 2 and
RT_DETR [28], alongside their finetuned variants trained on the SoccerNet
Player Tracking Challenge 2023 data [10], utilizing the Python package
Ultralytics. The outcomes are delineated in Tab. 2. The result signified that
YOLOv8 with finetuning yielded the most favorable player detection results.
The relatively diminished recall rate before finetuning was attributed to
player occlusion. However, with finetuning, YOLOv8 achieved a 97% AP, in
contrast to RT_DETR [28], which exhibited an increase in AP but a decrease in
precision.

Model | Precision | Recall | $\text{AP}_{0.5}$
---|---|---|---
YOLOv82 | 95.24% | 46.02% | 80.33%
RT_DETR [28] | 92.01% | 50.83% | 80.41%
RT_DETR [28] (finetuned) | 85.00% | 94.30% | 85.10%
YOLOv82 (finetuned) | 95.70% | 93.80% | 97.30%
Table 2: Soccer player detection models performance. Ranked by AP, with the
top-performing result highlighted in bold. Precision quantifies the proportion
of correctly detected bounding boxes out of all detected ones, recall measures
the proportion of correctly detected bounding boxes out of all ground truth
boxes, and AP calculates the area under the precision-recall curve.

Concerning tracking evaluation, the compared methodologies encompass the
baseline and the top-performing approach (Kalisteo [29]) from the SoccerNet
Player Tracking Challenge 2023 [10], in addition to two prevalent tracking
techniques integrated with the finetuned YOLOv8. The outcomes are depicted in
Tab. 3. It was evident that the Kalisteo method [29] continued to outperform
popular methods. However, the absence of publicly available code for the
Kalisteo method [29] posed a limitation. Furthermore, the test set of the
SoccerNet Player Tracking Challenge 2023 [10] comprised numerous clips
featuring substantial player occlusion, whereas in most shooting scenarios,
player occlusion was less pronounced. Hence, we ascertain that BoTSort
suffices for AutoSoccerPose.

Method | DetA | AssA | HOTA
---|---|---|---
SoccerNet Baseline [10] | 38.38% | 49.81% | 43.67%
YOLOv82-ByteTrack [56] | 63.75% | 42.31% | 51.85%
YOLOv82-BoTSort [1] | 68.80% | 53.40% | 60.54%
YOLOX [15]-Kalisteo [29] |  73.64% | 73.76% | 73.65%
Table 3: Tracking performance on SoccerNet player tracking challenge 2023 [10]
test set. The evluation metrics include DetA (Detection Accuracy), AssA
(Association Accuracy), and HOTA (Higher Order Tracking Accuracy) [27]. Ranked
by HOTA, with the top-performing result highlighted in bold.

###  5.2 Tracklet selection required features

| Required
---
Frames
Precision | Recall | ACC |  | CLIP
---
ACC
Frame 12-13 | 20.81% | 77.50% | 76.52% | 52.50%
Frame 11-14 | 38.33% | 57.50% | 90.02% | 65.00%
Frame 10-15 | 57.78% | 65.00% | 93.90% | 75.00%
Frame 9-16 | 41.82% | 57.50% | 90.94% | 65.00%
Table 4: Tracklet selection model performance utilizing different numbers of frames. Ranked by the number of required frames, with the top-performing result highlighted in bold. Model | Head | Sho | El | Wri | Body | Hip | Knee | Ank | PDJ | AUC
---|---|---|---|---|---|---|---|---|---|---
HRNet [41] | 94.56% | 63.38% | 21.14% | 17.29% | 96.35% | 78.68% | 47.46% | 34.94% | 56.08% | 47.38%
DWPose [52] | 86.69% | 85.80% | 72.75% | 58.88% | 89.88% | 71.41% | 45.41% | 32.98% | 67.94% | 55.80%
RTMPose [22] | 96.81% | 95.08% | 86.22% | 76.01% | 96.98% | 95.62% | 88.30% | 79.00% | 89.51% | 73.56%
Table 5: 2D pose estimation performance on 3DSP. Ranked by PDJ, with the top-
performing result highlighted in bold. Sho, El, Wri, and Ank represent
Shoulder, Elbow, Wrist, and Ankle respectively. Each column for body parts
indicates the average PDJ for the specified body part, while the PDJ column
denotes the mean PDJ across all 17 keypoints. Figure 3: AutoSoccerPose
qualitative results. Each row depicts the AutoSoccerPose 2D and 3D pose
estimation from a broadcast video, and the columns denote which frame of the
video. The top photo represents the broadcast video frame for each cell, while
the bottom-left and bottom-right images correspond to the 2D and 3D pose
estimations, respectively.

After verifying the tracking, the next crucial step is tracklet selection,
performed using a CNN model in AutoSoccerPose (see Sec. 3.3). Given the
absence of existing models tailored for this specific task, we focus on
validating essential features, namely the frames containing pivotal
information. Our evaluation criteria encompass precision and recall metrics,
alongside accuracy (ACC) and clip accuracy (CLIP ACC). The ACC metric assesses
the CNN model’s accuracy in classifying a tracklet as the shooter’s tracklet.
Additionally, CLIP ACC evaluates the accuracy of tracklet selection,
specifically in identifying the shooter’s tracklet from all tracklets within a
clip (see Sec. 3.3). 200 clips are sourced from SoccerNet [12, 17],
complemented by tracklet information from the 3DSP train set (see Sec. 4).

The feature set is based on tracklet images spanning frames 12 and 13,
bounding the timestep of the shot movement (see Sec. 3). We iteratively enrich
the feature set by progressively incorporating the frames preceding and
succeeding the aforementioned frames. For each feature set, the CNN model is
trained using PyTorch with an 80/20 train-validation split. Grid searching is
employed to optimize the number of convolutional layers, MLP layers, and the
hidden size of the MLP layer. The results are summarized in Table Tab. 4.
Notably, employing tracklet images from frames 10 to 15 yields the best
performance. While the tracklet selection model demonstrates reasonable
performance, there remains potential for further enhancement.

###  5.3 2D pose estimation performance on 3DSP

This subsection aims to evaluate the performance of 2D pose estimation models
under zero-shot conditions. We utilized the image and 2D annotations from the
3DSP dataset train set (see Sec. 4). The dataset provides annotations for 17
keypoints representing various body joints. To assess the performance, we
employed the Percent of Detected Joints (PDJ) evaluation metric [45]. In this
metric, a keypoint is considered detected if the normalized distance between
the predicted and ground truth keypoints is under 0.5. Additionally, the area
under the PDJ curve (AUC) was calculated to report the performance across
different thresholds.

Given the dynamic nature of the shooting movements, especially significant
rotation of the shoulders and hips during shots, the distance between the
contralateral shoulder and hip keypoints may not accurately represent the
torso’s largest distance consistently. Therefore, we normalized the distance
between detected keypoints and ground truth keypoints by the distance between
the center of the shoulder and the center of the hip, as suggested in [44].

The evaluated models include DWPose [52] and RTMPose [22], which are the top-
performing models in the COCO-WholeBody dataset [23] for 2D pose estimation
tasks. Additionally, we included HRNet [41], which is the 2D pose estimation
model integrated into MotionAGFormer [30]. The results were summarized in Tab.
5, indicating that RTMPose [22] exhibited the best performance across all
joint detections. The detection errors predominantly occurred in limb
detection. Through observation, two primary reasons were identified: when the
shooter faced backward, the left and right limbs could be misidentified, and
when the shooter faced sideways, one side of the limbs was covered by the
torso.

###  5.4 3D pose estimation on non-annotated data

Here, we validate the results obtained by AutoSoccerPose for estimating 3D
poses from broadcast videos, instead of relying on 2D annotations from 3DSP.
This approach directly showcases what AutoSoccerPose users can retrieve
without annotated data. Qualitative results are depicted in Fig. 3, comprising
the outcomes of three clips: a) a shot towards the goal on the right from
outside the box, b) a shot towards the goal on the left from outside the box,
and c) a shot towards the goal on the left from inside the box. While the
majority of 2D and 3D pose estimations were satisfactory, two limitations of
AutoSoccerPose were identified: motion blur, as in Clip (a) Frame 15, which
impeded the accurate detection of limb keypoints, and occlusion between the
shooter and other players, as in Clip (b) Frame 20, hindered the detection of
limb positions.

###  5.5 Soccer shots posture analysis

This subsection is dedicated to analyzing the 3D posture of professional
soccer players. We leverage the 3D pose data available in the 3DSP dataset
(see Sec. 4) to ensure high data fidelity for insightful analysis. Our goal is
to identify various shooting techniques. To achieve this, we employ the 3DSP-
GRAE model (see Sec. 3.5) to extract representative vectors from sequential
shot movements (comprising 3D postures from 20 frames). Subsequently, we apply
K-means clustering [3] to categorize these vectors. The resulting clusters are
visualized in a 2-dimensional space using t-distributed stochastic neighbor
embedding (t-SNE) [46] and illustrated in Fig. 4. Furthermore, we showcase the
shot movements closest to each cluster centroid by displaying frames 11 to 15
in Fig. 5. Based on observation, Cluster 1 and Cluster 2 represent inside
shots, whereas Cluster 3 denotes instep shots. The definitions are as follows:

Figure 4: K-means clustering visualization. Figure 5: Shot movement closest to
the cluster centroid. The row and column denote the cluster and frame,
respectively.

  * •

Inside shot: Precise technique where the ball is struck with the inside of the
foot, allowing for accurate placement into goal corners. Involves swinging the
foot inward toward the opposite foot for finesse shots.

  * •

Instep shot: Powerful technique involving striking the ball with the area of
the foot containing the laces. Commonly used for long-range or forceful shots.
Players approach the ball straight on and use the laces to generate
significant power.

Comparing Cluster 1 and 2 with statistics, Cluster 1 exhibited a 16% greater
shooting feet ankle average travel distance. Moreover, the maximum vertical
coordinate in Cluster 1 surpassed that of Cluster 2 by 33%, while the minimum
knee angle was 17% lower. Statistics revealed that shots in Cluster 1
generated a larger swing motion compared to those in Cluster 2, likely
attributable to deeper knee bending and a more pronounced vertical trajectory.
This enabled enhanced momentum transfer and force generation during the shot.
The identified shooting style holds promise for further utilization in player
evaluation for optimizing posture and style, as demonstrated in [49].
Meanwhile, when the aforementioned analysis was performed with linear model
PCA and K-means [3], the clustering result could only reflect shots to the
left, center, and right, demonstrating the need for a non-linear
spatiotemporal model (see GitHub1).

##  6 Conclusion

In this paper, we introduce AutoSoccerPose, a framework designed for
estimating both 2D and 3D shooting postures of soccer players and conducting
comprehensive posture analysis. Additionally, we present the 3DSP dataset,
which stands as the largest annotated 2D pose dataset for professional soccer
images. Through experimental validation, we showcase the performance of
AutoSoccerPose across each stage of the framework, culminating in an
assessment of its overall efficacy. Given that each component of
AutoSoccerPose corresponds to a well-established task in computer vision, the
seamless integration of future state-of-the-art models into AutoSoccerPose
would readily enhance its capabilities. We envision that our research will
inspire automated posture analysis in soccer and various other sports.

##  7 Acknowledgement

This work was financially supported by JSTSPRING GrantNumber JPMJSP2125 and
JSPS KAKENHI 23H03282.

## References

  * Aharon et al. [2022] Nir Aharon, Roy Orfaig, and Ben-Zion Bobrovsky.  Bot-sort: Robust associations multi-pedestrian tracking.  _arXiv preprint arXiv:2206.14651_ , 2022.
  * Andriluka et al. [2014] Mykhaylo Andriluka, Leonid Pishchulin, Peter Gehler, and Bernt Schiele.  2d human pose estimation: New benchmark and state of the art analysis.  In _Proceedings of the IEEE Conference on computer Vision and Pattern Recognition_ , pages 3686–3693, 2014.
  * Arthur et al. [2007] David Arthur, Sergei Vassilvitskii, et al.  k-means++: The advantages of careful seeding.  In _Soda_ , pages 1027–1035, 2007.
  * Bera et al. [2023] Asish Bera, Mita Nasipuri, Ondrej Krejcar, and Debotosh Bhattacharjee.  Fine-grained sports, yoga, and dance postures recognition: A benchmark analysis.  _IEEE Transactions on Instrumentation and Measurement_ , 2023.
  * Cao et al. [2023] Jinkun Cao, Jiangmiao Pang, Xinshuo Weng, Rawal Khirodkar, and Kris Kitani.  Observation-centric sort: Rethinking sort for robust multi-object tracking.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 9686–9696, 2023.
  * Cao et al. [2017] Zhe Cao, Tomas Simon, Shih-En Wei, and Yaser Sheikh.  Realtime multi-person 2d pose estimation using part affinity fields.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 7291–7299, 2017.
  * Cheng et al. [2020] Bowen Cheng, Bin Xiao, Jingdong Wang, Honghui Shi, Thomas S Huang, and Lei Zhang.  Higherhrnet: Scale-aware representation learning for bottom-up human pose estimation.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 5386–5395, 2020.
  * Chun et al. [2023] Sungho Chun, Sungbum Park, and Ju Yong Chang.  Learnable human mesh triangulation for 3d human pose and shape estimation.  In _Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision_ , pages 2850–2859, 2023.
  * Ci et al. [2019] Hai Ci, Chunyu Wang, Xiaoxuan Ma, and Yizhou Wang.  Optimizing network structure for 3d human pose estimation.  In _Proceedings of the IEEE/CVF international conference on computer vision_ , pages 2262–2271, 2019.
  * Cioppa et al. [2023] Anthony Cioppa, Silvio Giancola, Vladimir Somers, Floriane Magera, Xin Zhou, Hassan Mkhallati, Adrien Deliège, Jan Held, Carlos Hinojosa, Amir M Mansourian, et al.  Soccernet 2023 challenges results.  _arXiv preprint arXiv:2309.06006_ , 2023.
  * Decroos et al. [2020] Tom Decroos, Lotte Bransen, Jan Van Haaren, and Jesse Davis.  Vaep: an objective approach to valuing on-the-ball actions in soccer.  In _Proceedings of the twenty-ninth international joint conference on artificial intelligence, IJCAI-20_ , pages 4696–4700. International Joint Conferences on Artificial Intelligence Organization, 2020.
  * Deliege et al. [2021] Adrien Deliege, Anthony Cioppa, Silvio Giancola, Meisam J Seikavandi, Jacob V Dueholm, Kamal Nasrollahi, Bernard Ghanem, Thomas B Moeslund, and Marc Van Droogenbroeck.  Soccernet-v2: A dataset and benchmarks for holistic understanding of broadcast soccer videos.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 4508–4519, 2021.
  * Du et al. [2023] Yunhao Du, Zhicheng Zhao, Yang Song, Yanyun Zhao, Fei Su, Tao Gong, and Hongying Meng.  Strongsort: Make deepsort great again.  _IEEE Transactions on Multimedia_ , 2023.
  * Fang et al. [2024] Jiale Fang, Calvin Yeung, and Keisuke Fujii.  Foul prediction with estimated poses from soccer broadcast video.  _arXiv preprint arXiv:2402.09650_ , 2024.
  * Ge et al. [2021] Zheng Ge, Songtao Liu, Feng Wang, Zeming Li, and Jian Sun.  Yolox: Exceeding yolo series in 2021.  _arXiv preprint arXiv:2107.08430_ , 2021.
  * Geng et al. [2021] Zigang Geng, Ke Sun, Bin Xiao, Zhaoxiang Zhang, and Jingdong Wang.  Bottom-up human pose estimation via disentangled keypoint regression.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 14676–14686, 2021.
  * Giancola et al. [2018] Silvio Giancola, Mohieddine Amine, Tarek Dghaily, and Bernard Ghanem.  Soccernet: A scalable dataset for action spotting in soccer videos.  In _Proceedings of the IEEE conference on computer vision and pattern recognition workshops_ , pages 1711–1721, 2018.
  * Held et al. [2023] Jan Held, Anthony Cioppa, Silvio Giancola, Abdullah Hamdi, Bernard Ghanem, and Marc Van Droogenbroeck.  Vars: Video assistant referee system for automated soccer decision making from multiple views.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 5085–5096, 2023.
  * Hong et al. [2011] Sungchan Hong, Chulsoo Chung, Keiko Sakamoto, and Takeshi Asai.  Analysis of the swing motion on knuckling shot in soccer.  _Procedia Engineering_ , 13:176–181, 2011.
  * Ingwersen et al. [2023] Christian Keilstrup Ingwersen, Christian Møller Mikkelstrup, Janus Nørtoft Jensen, Morten Rieger Hannemose, and Anders Bjorholm Dahl.  Sportspose-a dynamic 3d sports pose dataset.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 5218–5227, 2023.
  * Iskakov et al. [2019] Karim Iskakov, Egor Burkov, Victor Lempitsky, and Yury Malkov.  Learnable triangulation of human pose.  In _Proceedings of the IEEE/CVF international conference on computer vision_ , pages 7718–7727, 2019.
  * Jiang et al. [2023] Tao Jiang, Peng Lu, Li Zhang, Ningsheng Ma, Rui Han, Chengqi Lyu, Yining Li, and Kai Chen.  Rtmpose: Real-time multi-person pose estimation based on mmpose.  _arXiv preprint arXiv:2303.07399_ , 2023.
  * Jin et al. [2020] Sheng Jin, Lumin Xu, Jin Xu, Can Wang, Wentao Liu, Chen Qian, Wanli Ouyang, and Ping Luo.  Whole-body human pose estimation in the wild.  In _Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part IX 16_ , pages 196–214. Springer, 2020.
  * Johnson and Everingham [2010] Sam Johnson and Mark Everingham.  Clustered pose and nonlinear appearance models for human pose estimation.  In _bmvc_ , page 5. Aberystwyth, UK, 2010.
  * Kipf and Welling [2016] Thomas N Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  _arXiv preprint arXiv:1609.02907_ , 2016.
  * Lin et al. [2014] Tsung-Yi Lin, Michael Maire, Serge Belongie, James Hays, Pietro Perona, Deva Ramanan, Piotr Dollár, and C Lawrence Zitnick.  Microsoft coco: Common objects in context.  In _Computer Vision–ECCV 2014: 13th European Conference, Zurich, Switzerland, September 6-12, 2014, Proceedings, Part V 13_ , pages 740–755. Springer, 2014.
  * Luiten et al. [2021] Jonathon Luiten, Aljosa Osep, Patrick Dendorfer, Philip Torr, Andreas Geiger, Laura Leal-Taixé, and Bastian Leibe.  Hota: A higher order metric for evaluating multi-object tracking.  _International journal of computer vision_ , 129:548–578, 2021.
  * Lv et al. [2023] Wenyu Lv, Shangliang Xu, Yian Zhao, Guanzhong Wang, Jinman Wei, Cheng Cui, Yuning Du, Qingqing Dang, and Yi Liu.  Detrs beat yolos on real-time object detection.  _arXiv preprint arXiv:2304.08069_ , 2023.
  * Maglo et al. [2023] Adrien Maglo, Astrid Orcesi, Julien Denize, and Quoc Cuong Pham.  Individual locating of soccer players from a single moving view.  _Sensors_ , 23(18):7938, 2023.
  * Mehraban et al. [2024] Soroush Mehraban, Vida Adeli, and Babak Taati.  Motionagformer: Enhancing 3d human pose estimation with a transformer-gcnformer network.  In _Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision_ , pages 6920–6930, 2024.
  * Mehta et al. [2017] Dushyant Mehta, Helge Rhodin, Dan Casas, Pascal Fua, Oleksandr Sotnychenko, Weipeng Xu, and Christian Theobalt.  Monocular 3d human pose estimation in the wild using improved cnn supervision.  In _2017 international conference on 3D vision (3DV)_ , pages 506–516. IEEE, 2017.
  * Meinhardt et al. [2022] Tim Meinhardt, Alexander Kirillov, Laura Leal-Taixe, and Christoph Feichtenhofer.  Trackformer: Multi-object tracking with transformers.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 8844–8854, 2022.
  * Nakamura et al. [2018] Taishu Nakamura, Tasuku Miyoshi, Shota Sato, Motoki Takagi, Yasuhisa Kamada, and Yoshiyuki Kobayashi.  Differences in soccer kicking type identified using principal component analysis.  _Sports Engineering_ , 21:149–159, 2018.
  * Pappalardo et al. [2019] Luca Pappalardo, Paolo Cintia, Alessio Rossi, Emanuele Massucco, Paolo Ferragina, Dino Pedreschi, and Fosca Giannotti.  A public data set of spatio-temporal match events in soccer competitions.  _Scientific data_ , 6(1):236, 2019.
  * Pinheiro et al. [2021] Guilherme de Sousa Pinheiro, Vitor Bertoli Nascimento, Matt Dicks, Varley Teoldo Costa, and Martin Lames.  Design and validation of an observational system for penalty kick analysis in football (ospaf).  _Frontiers in Psychology_ , 12:661179, 2021.
  * Pinheiro et al. [2022] Guilherme de Sousa Pinheiro, Xing Jin, Varley Teoldo Da Costa, and Martin Lames.  Body pose estimation integrated with notational analysis: A new approach to analyze penalty kicks strategy in elite football.  _Frontiers in sports and active living_ , 4:818556, 2022.
  * Reddy et al. [2021] N Dinesh Reddy, Laurent Guigues, Leonid Pishchulin, Jayan Eledath, and Srinivasa G Narasimhan.  Tessetrack: End-to-end learnable multi-person articulated 3d pose tracking.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 15190–15200, 2021.
  * Redmon et al. [2016] Joseph Redmon, Santosh Divvala, Ross Girshick, and Ali Farhadi.  You only look once: Unified, real-time object detection.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 779–788, 2016.
  * Scott et al. [2022] Atom Scott, Ikuma Uchida, Masaki Onishi, Yoshinari Kameda, Kazuhiro Fukui, and Keisuke Fujii.  Soccertrack: A dataset and tracking algorithm for soccer with fish-eye and drone videos.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 3569–3579, 2022.
  * Srivastava et al. [2015] Nitish Srivastava, Elman Mansimov, and Ruslan Salakhudinov.  Unsupervised learning of video representations using lstms.  In _International conference on machine learning_ , pages 843–852. PMLR, 2015.
  * Sun et al. [2019] Ke Sun, Bin Xiao, Dong Liu, and Jingdong Wang.  Deep high-resolution representation learning for human pose estimation.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 5693–5703, 2019.
  * Suzuki et al. [2023] Tomohiro Suzuki, Kazushi Tsutsui, Kazuya Takeda, and Keisuke Fujii.  Runner re-identification from single-view video in the open-world setting.  _arXiv preprint arXiv:2310.11700_ , 2023.
  * Tanaka et al. [2023] Ryota Tanaka, Tomohiro Suzuki, Kazuya Takeda, and Keisuke Fujii.  Automatic edge error judgment in figure skating using 3d pose estimation from a monocular camera and imus.  In _Proceedings of the 6th International Workshop on Multimedia Content Analysis in Sports_ , pages 41–48, 2023.
  * Tompson et al. [2014] Jonathan J Tompson, Arjun Jain, Yann LeCun, and Christoph Bregler.  Joint training of a convolutional network and a graphical model for human pose estimation.  _Advances in neural information processing systems_ , 27, 2014.
  * Toshev and Szegedy [2014] Alexander Toshev and Christian Szegedy.  Deeppose: Human pose estimation via deep neural networks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 1653–1660, 2014.
  * Van der Maaten and Hinton [2008] Laurens Van der Maaten and Geoffrey Hinton.  Visualizing data using t-sne.  _Journal of machine learning research_ , 9(11), 2008.
  * Wang et al. [2020] Jingbo Wang, Sijie Yan, Yuanjun Xiong, and Dahua Lin.  Motion guided 3d pose estimation from videos.  In _European Conference on Computer Vision_ , pages 764–780. Springer, 2020.
  * Wang et al. [2011] Yang Wang, Duan Tran, and Zicheng Liao.  Learning hierarchical poselets for human parsing.  In _CVPR 2011_ , pages 1705–1712. IEEE, 2011.
  * Wear et al. [2022] Matthew Wear, Ryan Beal, Tim Matthews, Tim Norman, and Sarvapali Ramchurn.  Learning from the pros: Extracting professional goalkeeper technique from broadcast footage.  _arXiv preprint arXiv:2202.12259_ , 2022.
  * Wojke et al. [2017] Nicolai Wojke, Alex Bewley, and Dietrich Paulus.  Simple online and realtime tracking with a deep association metric.  In _2017 IEEE International Conference on Image Processing (ICIP)_ , pages 3645–3649. IEEE, 2017.
  * Xu et al. [2020] Yihong Xu, Aljosa Osep, Yutong Ban, Radu Horaud, Laura Leal-Taixé, and Xavier Alameda-Pineda.  How to train your deep multi-object tracker.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 6787–6796, 2020.
  * Yang et al. [2023] Zhendong Yang, Ailing Zeng, Chun Yuan, and Yu Li.  Effective whole-body pose estimation with two-stages distillation.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 4210–4220, 2023.
  * Yeung and Fujii [2023] Calvin Yeung and Keisuke Fujii.  A strategic framework for optimal decisions in football 1-vs-1 shot-taking situations: An integrated approach of machine learning, theory-based modeling, and game theory.  _arXiv preprint arXiv:2307.14732_ , 2023.
  * Yeung et al. [2023] Calvin Yeung, Tony Sit, and Keisuke Fujii.  Transformer-based neural marked spatio temporal point process model for football match events analysis.  _arXiv preprint arXiv:2302.09276_ , 2023.
  * Zhang et al. [2019] Song-Hai Zhang, Ruilong Li, Xin Dong, Paul Rosin, Zixi Cai, Xi Han, Dingcheng Yang, Haozhi Huang, and Shi-Min Hu.  Pose2seg: Detection free human instance segmentation.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 889–898, 2019.
  * Zhang et al. [2022] Yifu Zhang, Peize Sun, Yi Jiang, Dongdong Yu, Fucheng Weng, Zehuan Yuan, Ping Luo, Wenyu Liu, and Xinggang Wang.  Bytetrack: Multi-object tracking by associating every detection box.  In _European Conference on Computer Vision_ , pages 1–21. Springer, 2022.
  * Zhang et al. [2023] Yuang Zhang, Tiancai Wang, and Xiangyu Zhang.  Motrv2: Bootstrapping end-to-end multi-object tracking by pretrained object detectors.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 22056–22065, 2023.
  * Zhou et al. [2017] Xingyi Zhou, Qixing Huang, Xiao Sun, Xiangyang Xue, and Yichen Wei.  Towards 3d human pose estimation in the wild: a weakly-supervised approach.  In _Proceedings of the IEEE international conference on computer vision_ , pages 398–407, 2017.


n\epsilon_{n}\int_{-\infty}^{0}G_{n}(t)dt\leq 4e\pi cn^{-\alpha+\frac{1}{2}},$
(41)

where (a) follows by using Lemma 10.b. This proves the probability bound
presented in Eq. (13). Next, we provide a bound on the tail probability. We
also have,

$\displaystyle\int_{-\infty}^{0}$ $\displaystyle
G_{n}(t)dt\allowdisplaybreaks$
$\displaystyle=\int_{-\infty}^{0}\exp\Big{(}-n\epsilon_{n}t-n\rho_{n}(e^{-t}+t-1)\Big{)}dt\allowdisplaybreaks$
$\displaystyle=\int_{0}^{\infty}\exp\Big{(}n\epsilon_{n}t-n\rho_{n}(e^{t}-t-1)\Big{)}dt\allowdisplaybreaks$
$\displaystyle\geq\exp\big{(}\zeta_{n}\big{)}\int_{1/\sqrt{n\rho_{n}}}^{\infty}\exp\Big{(}-n\rho_{n}(e^{t}-t-1)\Big{)}dt\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\geq}}\exp\big{(}\zeta_{n}\big{)}\int_{1/\sqrt{n\rho_{n}}}^{\infty}\exp\Big{(}-\frac{1}{2}n\rho_{n}t^{2}-n\rho_{n}t^{3}e^{t}\Big{)}dt\allowdisplaybreaks$
$\displaystyle\geq\exp\big{(}\zeta_{n}\big{)}\int_{1/\sqrt{n\rho_{n}}}^{1}\exp\Big{(}-\frac{1}{2}n\rho_{n}t^{2}-en\rho_{n}t^{2}\Big{)}dt\allowdisplaybreaks$
$\displaystyle\geq\exp\big{(}\zeta_{n}\big{)}\int_{1/\sqrt{n\rho_{n}}}^{1}\exp\big{(}-4n\rho_{n}t^{2}\big{)}dt\allowdisplaybreaks$
$\displaystyle=\frac{1}{\sqrt{n\rho_{n}}}\exp\big{(}\zeta_{n}\big{)}\int_{1}^{\sqrt{n\rho_{n}}}\exp\big{(}-4t^{2}\big{)}dt,$

where (a) follows by using $e^{t}-t-1\leq\frac{1}{2}t^{2}+t^{3}e^{t}$ for all
$t\geq 0$. This gives us that

$\displaystyle\int_{-\infty}^{0}$ $\displaystyle
G_{n}(t)dt\geq\frac{1}{\tilde{\kappa}_{1}}\times\frac{1}{\sqrt{n\rho_{n}}}\exp\big{(}\zeta_{n}\big{)},$
(42)

where
$\tilde{\kappa}_{1}:=\big{(}\int_{1}^{2}\exp\big{(}-4t^{2}\big{)}dt\big{)}^{-1}$.
Further, we have that for any $\theta\in(0,\sqrt{n\rho_{n}})$,

$\displaystyle\int_{0}^{\theta/\sqrt{n\rho_{n}}}$ $\displaystyle
G_{n}(t)dt\allowdisplaybreaks$
$\displaystyle=\int_{0}^{\theta/\sqrt{n\rho_{n}}}\exp\Big{(}-n\epsilon_{n}t-n\rho_{n}(e^{-t}+t-1)\Big{)}dt\allowdisplaybreaks$
$\displaystyle\leq\int_{0}^{\theta/\sqrt{n\rho_{n}}}\exp\Big{(}-\frac{1}{2}n\rho_{n}t^{2}-n\rho_{n}\big{(}e^{-t}-\frac{1}{2}t^{2}+t-1\big{)}\Big{)}dt\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\int_{0}^{\theta/\sqrt{n\rho_{n}}}\exp\Big{(}-\frac{1}{2}n\rho_{n}t^{2}+\frac{1}{6}n\rho_{n}t^{3}\Big{)}dt\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}\frac{1}{\sqrt{n\rho_{n}}}\int_{0}^{\theta}\exp\Big{(}-\frac{1}{2}t^{2}+\frac{1}{6\sqrt{n\rho_{n}}}t^{3}\Big{)}dt\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}\frac{1}{\sqrt{n\rho_{n}}}\int_{0}^{\theta}\exp\Big{(}-\frac{1}{3}t^{2}\Big{)}dt\allowdisplaybreaks$
$\displaystyle\leq\frac{1}{\sqrt{n\rho_{n}}}\int_{0}^{\infty}\exp\Big{(}-\frac{1}{3}t^{2}\Big{)}dt\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(d)}}{{\leq}}\frac{\sqrt{3\pi}}{2\sqrt{n\rho_{n}}},$
(43)

where (a) follows by using $e^{-t}-\frac{1}{2}t^{2}+t-1\geq-\frac{1}{6}t^{3}$;
(b) follows by substituting $t\rightarrow\frac{1}{\sqrt{n\rho_{n}}}t$; (c)
follows for $\theta<\sqrt{n\rho_{n}}$; and finally, (d) follows by using Eq.
(D.1). Now, by using Lemma 10.a and substituting
$\tilde{r}_{n}=r_{n}-n\epsilon_{n}$, for $\theta\in(0,\sqrt{n\rho_{n}})$, we
have

$\displaystyle\mathbb{E}$
$\displaystyle\Big{[}\exp\Big{(}\frac{1}{\sqrt{n\rho_{n}}}\theta\tilde{r}_{n}\Big{)}\Big{|}r_{n}>0\Big{]}\allowdisplaybreaks$
$\displaystyle=\exp\Big{(}n\rho_{n}(e^{-\frac{\theta}{\sqrt{n\rho_{n}}}}+\frac{\theta}{\sqrt{n\rho_{n}}}-1)\Big{)}\Big{[}1+\left(\int_{-\infty}^{0}G_{n}(t)dt\right)^{-1}\int_{0}^{\frac{\theta}{\sqrt{n\rho_{n}}}}G_{n}(t)dt\Big{]}\allowdisplaybreaks$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\exp\Big{(}\frac{\theta^{2}}{2}\Big{)}\Big{[}1+\tilde{\kappa}_{1}\frac{\sqrt{3\pi}}{2}\exp\big{(}-\zeta_{n}\big{)}\Big{]},$
(44)

where (a) follows by using $e^{-x}+x-1\leq\frac{x^{2}}{2}$ for any $x\geq 0$
and Eq. (42) and (D.1). Note that the result in Eq. (D.1) is valid for all
three regimes.

Next, by using Markov’s inequality, for any $\theta>0$ and $x\geq 0$,

$\displaystyle\mathbb{P}\Big{(}\frac{1}{\sqrt{n\rho_{n}}}r_{n}>x+\zeta_{n}\Big{|}r_{n}>0\Big{)}\leq\Big{[}1+\tilde{\kappa}_{1}\frac{\sqrt{3\pi}}{2}\Big{]}\exp\Big{(}-\theta
x+\frac{\theta^{2}}{2}\Big{)}.$

As $r_{n}\leq n$, we only consider the range of $x$ for which
$x+\zeta_{n}\leq\frac{1}{\sqrt{n\rho_{n}}}n$, which gives
$x<\sqrt{n\rho_{n}}$. Then, by substituting $\theta=x<\sqrt{n\rho_{n}}$ in the
previous equation, we have

$\displaystyle\mathbb{P}\Big{(}\frac{1}{\sqrt{n\rho_{n}}}r_{n}>x+\zeta_{n}\Big{|}r_{n}>0\Big{)}\leq\Big{[}1+\tilde{\kappa}_{1}\frac{\sqrt{3\pi}}{2}\Big{]}e^{-\frac{x^{2}}{2}}.$
(45)

Now, result in Theorem 8.a follows by using Eq. (41) and picking

$\kappa_{1}:=4e\pi\Big{[}1+\tilde{\kappa}_{1}\frac{\sqrt{3\pi}}{2}\Big{]}.$

This completes the proof of Theorem 8.a.

Next, we consider the HW regime. In this case, we have
$\epsilon_{n}=cn^{-\frac{1}{2}}$, and

$\displaystyle\zeta_{n}=\frac{c}{\sqrt{1-cn^{-\frac{1}{2}}}}.$

Further,

$\displaystyle\int_{-\infty}^{0}G_{n}(t)dt-\sqrt{\frac{2\pi}{n\rho_{n}}}\exp\Big{(}\frac{1}{2}\zeta^{2}\Big{)}\Phi\big{(}\zeta_{n}\big{)}\allowdisplaybreaks$
$\displaystyle=\int_{0}^{\infty}\exp\Big{(}n\epsilon_{n}t-n\rho_{n}(e^{t}-t-1)\Big{)}dt-\int_{0}^{\infty}\exp\Big{(}n\epsilon_{n}t-\frac{1}{2}n\rho_{n}t^{2}\Big{)}dt\allowdisplaybreaks$
$\displaystyle=\int_{0}^{\infty}\exp\Big{(}n\epsilon_{n}t-\frac{1}{2}n\rho_{n}t^{2}\Big{)}\Big{(}\exp\Big{(}-n\rho_{n}\big{(}e^{t}-\frac{1}{2}t^{2}-t-1\big{)}\Big{)}-1\Big{)}dt$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\geq}}-n\rho_{n}\int_{0}^{\infty}\big{(}e^{t}-\frac{1}{2}t^{2}-t-1\big{)}\exp\Big{(}n\epsilon_{n}t-\frac{1}{2}n\rho_{n}t^{2}\Big{)}dt$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\geq}}-n\rho_{n}\int_{0}^{\infty}t^{3}e^{t}\exp\Big{(}n\epsilon_{n}t-\frac{1}{2}n\rho_{n}t^{2}\Big{)}dt$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\geq}}-\frac{1}{n\rho_{n}}\int_{0}^{\infty}t^{3}\exp\Big{(}(1+2c)t-\frac{1}{2}t^{2}\Big{)}dt$
$\displaystyle\stackrel{{\scriptstyle(d)}}{{\geq}}-\frac{\tilde{\kappa}_{2}}{n\rho_{n}},$

where (a) follows by using $e^{x}-1\geq x$ for all $x\in\mathbb{R}$; (b)
follows by using $e^{t}-\frac{1}{2}t^{2}-t-1\leq t^{3}e^{t}$ for all $t\geq
0$; (c) follows by substituting $t\rightarrow\frac{1}{\sqrt{n\rho_{n}}}t$ and
using $\frac{1}{\sqrt{n\rho_{n}}}+\zeta_{n}\leq(1+2c)$ when the system is in
HW regime and $n\geq 4c^{2}$; and (d) follows by taking
$\tilde{\kappa}_{2}:=\int_{0}^{\infty}t^{3}\exp\Big{(}(1+2c)t-\frac{1}{2}t^{2}\Big{)}dt$.
This gives us that

$\displaystyle
n\epsilon_{n}\int_{-\infty}^{0}G_{n}(t)dt\geq\sqrt{2\pi}\zeta_{n}\exp\Big{(}\frac{1}{2}\zeta^{2}\Big{)}\Phi\left(\zeta_{n}\right)-\tilde{\kappa}_{2}\frac{\epsilon_{n}}{\rho_{n}}.$

Combining the above inequality with Eq. (40), we have

$\displaystyle-\tilde{\kappa}_{2}\frac{\epsilon_{n}}{\rho_{n}}\leq
n\epsilon_{n}\int_{-\infty}^{0}G_{n}(t)dt-\sqrt{2\pi}\zeta_{n}\exp\Big{(}\frac{1}{2}\zeta^{2}\Big{)}\Phi\big{(}\zeta_{n}\big{)}\leq
0.$

Thus, we have

$\displaystyle\lim_{n\rightarrow\infty}n\epsilon_{n}\int_{-\infty}^{0}$
$\displaystyle
G_{n}(t)dt=\lim_{n\rightarrow\infty}\sqrt{2\pi}\zeta_{n}\exp\Big{(}\frac{1}{2}\zeta_{n}^{2}\Big{)}\Phi\left(\zeta_{n}\right)=\sqrt{2\pi}c\exp\left(\frac{c^{2}}{2}\right)\Phi\left(c\right),$

and so,

$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}(r_{n}>0)=\lim_{n\rightarrow\infty}\frac{n\epsilon_{n}\int_{-\infty}^{0}G_{n}(t)dt}{1+n\epsilon_{n}\int_{-\infty}^{0}G_{n}(t)dt}=\frac{\sqrt{2\pi}c\exp\left(\frac{c^{2}}{2}\right)\Phi\left(c\right)}{1+\sqrt{2\pi}c\exp\left(\frac{c^{2}}{2}\right)\Phi\left(c\right)}.$

Next, by using Eq. (D.1) and Markov’s inequality, we have

$\displaystyle\mathbb{P}\Big{(}\frac{1}{\sqrt{n\rho_{n}}}r_{n}>x+\zeta_{n}\Big{)}$
$\displaystyle\leq\mathbb{P}\Big{(}\frac{1}{\sqrt{n\rho_{n}}}r_{n}>x+\zeta_{n}\Big{|}r_{n}>0\Big{)}$
$\displaystyle\leq\Big{[}1+\tilde{\kappa}_{1}\frac{\sqrt{3\pi}}{2}\Big{]}e^{-\frac{x^{2}}{2}},$

where the last inequality follows by same calculations as in Eq. (45). Now,
result in Theorem 8.b follows by picking

$\kappa_{2}:=\Big{[}1+\tilde{\kappa}_{1}\frac{\sqrt{3\pi}}{2}\Big{]}.$

This completes the proof of Theorem 8.b.

Finally, we consider the Sub-HW regime, i.e., $\epsilon_{n}=cn^{-\alpha}$
where $\alpha\in\Big{(}0,\frac{1}{2}\Big{)}$. In this case,

$\displaystyle\zeta_{n}=\Big{(}\frac{c^{2}n^{1-2\alpha}}{1-cn^{-\alpha}}\Big{)}^{\frac{1}{2}}\geq
cn^{\frac{1}{2}-\alpha}.$

As such $\zeta\rightarrow\infty$ as $n\rightarrow\infty$. From Eq. (42), we
have that

$\displaystyle n\epsilon_{n}\int_{-\infty}^{0}$ $\displaystyle
G_{n}(t)dt\geq\frac{1}{\tilde{\kappa}_{1}}\zeta_{n}\exp\big{(}\zeta_{n}\big{)}.$

Then,

$\displaystyle
1-\mathbb{P}(r_{n}>0)=\frac{1}{1+n\epsilon_{n}\int_{-\infty}^{0}G_{n}(t)dt}\leq\tilde{\kappa}_{1}\zeta_{n}^{-1}\exp\big{(}-\zeta_{n}\big{)}\leq\frac{\tilde{\kappa}_{1}}{c}n^{\alpha-\frac{1}{2}}e^{-cn^{\frac{1}{2}-\alpha}}.$

Thus, for $\alpha\in\big{(}0,\frac{1}{2}\big{)}$, by choosing
$\kappa_{3}=\tilde{\kappa}_{1}$, we have

$\displaystyle\mathbb{P}(r_{n}>0)\geq
1-\frac{\kappa_{3}}{c}n^{\alpha-\frac{1}{2}}e^{-cn^{\frac{1}{2}-\alpha}},$
$\displaystyle\lim_{n\rightarrow\infty}\mathbb{P}(r_{n}>0)=1.$

Next, by using Eq. (D.1) and Markov’s inequality, we have

$\displaystyle\mathbb{P}\Big{(}\frac{1}{\sqrt{n\rho_{n}}}\tilde{r}_{n}>x\Big{)}$
$\displaystyle\leq\mathbb{P}\Big{(}\frac{1}{\sqrt{n\rho_{n}}}r_{n}>x\Big{|}r_{n}>0\Big{)}$
$\displaystyle\leq\Big{[}1+\tilde{\kappa}_{1}\frac{\sqrt{3\pi}}{2}e^{-cn^{\frac{1}{2}-\alpha}}\Big{]}e^{-\frac{x^{2}}{2}}.$

Further, for $\theta>0$, we have
$\left(\int_{-\infty}^{0}G_{n}(t)dt\right)^{-1}\int_{\infty}^{-\theta}G_{n}(t)dt\leq
1$. Then,

$\displaystyle\mathbb{E}\Big{[}\exp\Big{(}-\frac{1}{\sqrt{n\rho_{n}}}\theta\tilde{r}_{n}\Big{)}\Big{|}r_{n}>0\Big{]}$
$\displaystyle\leq\exp\Big{(}n\rho_{n}\big{(}e^{\frac{\theta}{\sqrt{n\rho_{n}}}}-\frac{\theta}{\sqrt{n\rho_{n}}}-1\big{)}\Big{)}$
$\displaystyle\leq\exp\Big{(}\frac{1}{2}\theta^{2}+\frac{\theta^{3}}{\sqrt{n\rho_{n}}}e^{\frac{\theta}{\sqrt{n\rho_{n}}}}\Big{)},$

This gives us that

$\displaystyle\mathbb{P}\Big{(}-\frac{1}{\sqrt{n\rho_{n}}}\tilde{r}_{n}>x\Big{|}r_{n}>0\Big{)}\leq
e^{-\theta
x}\exp\Big{(}\frac{1}{2}\theta^{2}+\frac{\theta^{3}}{\sqrt{n\rho_{n}}}e^{\frac{\theta}{\sqrt{n\rho_{n}}}}\Big{)}.$

By substituting $\theta=x$ in the above equation, we get that for any $x\geq
0$,

$\displaystyle\mathbb{P}\Big{(}-\frac{1}{\sqrt{n\rho_{n}}}\tilde{r}_{n}>x\Big{|}r_{n}>0\Big{)}\leq\exp\Big{(}-\frac{1}{2}x^{2}+\frac{x^{3}}{\sqrt{n\rho_{n}}}e^{\frac{x}{\sqrt{n\rho_{n}}}}\Big{)}.$

Note that $-\tilde{r}_{n}\leq n\epsilon_{n}$ as $r_{n}\geq 0$. Hence, we
consider the range of $x$ such that $\sqrt{n\rho_{n}}x<n\epsilon_{n}$, or
$x<\zeta_{n}$. For any $x\geq\zeta_{n}$, the above inequality holds trivially.
Then, for $x<\zeta_{n}$, we have
$\frac{x}{\sqrt{n\rho_{n}}}\leq\frac{\epsilon_{n}}{\rho_{n}}\leq
2\epsilon_{n}\leq 1$, assuming $n^{\alpha}>2c$ (which holds when $n\geq
4c^{2}$). Then,

$\displaystyle\mathbb{P}\Big{(}-\frac{1}{\sqrt{n\rho_{n}}}\tilde{r}_{n}>x\Big{|}r_{n}>0\Big{)}\leq\exp\Big{(}-\frac{1}{2}x^{2}+2e\epsilon_{n}x^{2}\Big{)}.$

This completes the proof.

∎


#  Curvature dominance DE-model in f(R)-gravity

G. K. Goswami<EMAIL_ADDRESS>Department of Mathematics, Netaji Subhas
University of Technology, New Delhi-110 078, India  Rita Rani
<EMAIL_ADDRESS>Department of Mathematics, Netaji Subhas University of
Technology, New Delhi-110 078, India  Harshna Balhara<EMAIL_ADDRESS>Department of Mathematics, Netaji Subhas University of Technology, New
Delhi-110 078, India  J. K. Singh<EMAIL_ADDRESS>Department of
Mathematics, Netaji Subhas University of Technology, New Delhi-110 078, India

###### Abstract

We have probed a cosmological model in $f(R)$-gravity, which is a cubic
equation in scalar curvature $R$. The terms arise due to nonlinear $f(R)$
function are treated as energy due to curvature inspired geometry. As a
result, we find accelerating expansion in the universe, which creates an anti-
gravitating negative pressure in it. Some of the physical parameters are
solved using numerical methods. The evolution of the model are examined by the
latest observational Hubble data (46-data points) and Pantheon data (the
latest compilation of SNIa with 40 binned in the redshift range
$0.014\leqslant z\leqslant 1.62$). Some important features of the model have
been discussed by analyzing the plots of various dynamical parameters. The
plots of deceleration parameter $q$ and the Hubble parameter $H$ describe the
accelerating expansion in the evolution of the Universe at the present epoch.
The transition from deceleration to acceleration for our model is obtained at
redshift $z_{tr}\simeq 0.694069$, which is in good agreement with
$\Lambda$CDM. We have also carried out state finder analysis for our model.
The analysis of specific features of the model confirms that our model is
consistent with $\Lambda$CDM in late times.

PACS number: 98.80 cq
Keywords: $f(R)$-gravity, FLRW metric, Accelerating universe, State finder
diagnostics.

##  I Introduction

A study of origin, evolution and the large scale structure of the universe is
the subject of cosmology. It is a well established fact that our universe is
spatially homogeneous and isotropic (SHI) at large scale in the unit of
$Megaparsecs$ and $Giga\,years$. A veteran Friedman-Lema$\hat{i}$tre-
Robertson-Walker $(FLRW)$ space-time represents the SHI universe. When
Einstein Field Equations are solved for FLRW space-time, they describe an
expanding and decelerating universe which means that the various galaxies
present in it are moving away from us at a slowly decreasing pace with time.
The present day surveys and observation ref1 ; ref2 ; ref3 ; ref4 ; ref5
suggest a different type of scenario which tell us that the various galaxies
are moving away from us at a relatively increasing pace with time. The
literatures tell us that cold dark matters of gravitating origin and dark
energy of repulsive anti gravitating origin are present in abundance in our
universe which are responsible for the change in the mode of the universe.

The cosmological constant once introduced by Einstein in his field equations
to describe a physical static universe was resurrected in order to explain
acceleration in it known as $\Lambda$CDM model, which fits best on
observational ground ref6 ; ref7 ; ref8 ; ref9 ; ref10 ; ref11 . Despite, it
fails in explaining fine-tuning that vacuum energy is very small compared to
typical particle physics scales ref12 . Cosmological tracking solutions ref13
; ref14  were proposed to overcome the issues against the $\Lambda$CDM. By
that time, a different way of thinking was developed. It was thought that Dark
energy must arise from gravitational origin and some non-linear term of Ricci
scalar $R$ in the form of an arbitrary function $f(R)$ was introduced in the
Einstein Hilbert Action to arrive at an alternative to Einstein theory. The
name was given $f(R)$ gravity. In this theory acceleration is obtained
spontaneously from the gravitational sector.

In 2003, the $f(R)=R+\frac{a}{R}+{b}{R^{2}}$ formalism was first enunciated by
Nojiri and Odintsov ref15  to explain the early inflation and the acceleration
in late times. Starobinsky ref16  used a different form of $f(R)$ and
developed a cosmology without a cosmological constant. Sotiriou and Liberati
ref17 ; ref18  presented a Palatini Metric-affine formalism of $f(R)$ gravity.
It is noteworthy to mention the work of Srivastava ref19  who has developed a
model in $f(R)$ gravity, which represents both early (inflation) and late
acceleration in the Universe. Some remarkable works on $f(R)$-gravity have
also been carried out by many authors, which are stated in ref20 ; ref21 ;
ref22 ; ref23 ; ref24 ; ref25 ; ref26 .

In this paper, we have developed a model filled with perfect fluid in $f(R)$
gravity by taking $R+aR^{2}+{b}{R^{3}}$ as a particular form of $f(R)$. In the
background FLRW space-time, we have considered two field equations which
resemble Einstein field equations (EFE). The additional terms which arise due
to non-linear $f(R)$ function are shifted to the right-hand side of field
equations that are treated as energy expressions due to curvature inspired
geometry. As a result, they bring acceleration in the universe and put anti-
gravitating negative pressure in it. The paper is presented in the following
sequences. In section II, we have developed $f(R)$ field equations for FLRW
space-time, then taking a particular case $f(R)=R+\alpha R^{2}+\beta R^{3}$,
we have found numerical solutions for Hubble, deceleration and jerk
parameters. In section III, we have used a data set of $46$ Hubble parameters
to estimate present value of matter energy parameter $(\Omega_{m})$, density
$(\rho_{m})$ and equation of state (EoS) parameter for curvature
$(\omega_{k})$. On the basis of these estimated values, we present Error plots
for Hubble parameters and distant modulus. The plots show that our theoretical
plots fit well with the observed Hubble and Pantheon data. Moreover, the plots
match at par with the $\Lambda$CDM model. The deceleration parameter $q$
versus redshift $z$ plot describes the accelerating universe at the present
epoch. The red shift transition $z_{t}\simeq 0.694069$ in our model causes
transient acceleration in early and late times, which is in good agreement
with $\Lambda$CDM. In section III, we have also carried out the state finder
analysis for our model. The analysis confirms that our model lies near
$\Lambda$CDM, and it fits well on observational ground. In the last section,
we interpret our work for our obtained curvature dominated model.

##  II  Einstein Field Equations in $f(R)$-formalism

The $f(R)$-gravity general action with matter Lagrangian $S_{m}$ yields ref15

$S=\int\Big{(}\frac{1}{16\pi G}f(R)+S_{m}\Big{)}\sqrt{-g}dx^{4},$ (1)

where the function $f(R)=R+\alpha R^{2}+\beta R^{3}$ is used for inflation,
which was first suggested by Starobinsky ref16 . The function $f(R)$ can be
treated as a quantum correction in General Relativity (GR), where $R^{2}$
produces an inflationary scenario in the early times. The Starobinsky model is
most suitable model according as the latest observations pla  and may be
studied to describe inflation in addition to scalar field models jdb  .
Clearly, the function $f(R)$ containing negative powers of $R$ is capable to
explain the current accelerating expansion of the Universe.

The Einstein Field Equation (EFE) in $f(R)$-gravity are given as:

${}f^{\prime}(R)R_{ij}-\frac{1}{2}f(R)g_{ij}=T_{ij}-g_{ij}\square
f^{\prime}(R)+\nabla_{i}\nabla_{j}f^{\prime}(R),$ (2)

where we have taken $8\pi G$ and velocity of light c as unity. In this
equation $f^{\prime}(R)$ is the differentiation of $f(R)$ with respect to $R$
and other mathematical symbols have their usual meanings.

We consider $FLRW$ flat space-time, which is a SHI filled with perfect fluid
is given by

${}ds^{2}=dt^{2}-a^{2}(t)(dx^{2}+dy^{2}+dz^{2}).$ (3)

The energy momentum tensor of perfect fluid is as follows

$T_{ij}=(\rho_{m}+p_{m})u_{i}u_{j}-pg_{ij},$ (4)

where $\rho_{m}$ and $p_{m}$ are the matter density and the pressure of
perfect fluid. We use co moving coordinates so that the four velocity $u^{i}$
satisfies $u^{i}u_{i}=1$.

The field equations given by Eq.(2) are solved for metric Eq. (3) and energy
momentum tensor given by Eq.(4) are obtained as

${}3\frac{\dot{a^{2}}}{a^{2}}=\frac{\rho_{m}}{f^{\prime}}+\frac{1}{f^{\prime}}\bigg{[}\frac{1}{2}(f-f^{\prime}R)-3\frac{\dot{a}}{a}\dot{R}f^{\prime\prime}\bigg{]}$
(5)

and

${}2\frac{\ddot{a}}{a}+\frac{\dot{a^{2}}}{a^{2}}=-\frac{1}{f^{\prime}}\bigg{[}3\frac{\dot{a}}{a}f^{\prime\prime}\dot{R}+f^{\prime\prime\prime}\dot{R}^{2}+f^{\prime\prime}\ddot{R}+\frac{Rf^{\prime}}{2}-\frac{f}{2}\bigg{]},$
(6)

where prime and dot stand for derivative w.r.t. $R$ and proper time $t$
respectively. At present, we take pressure-less fluid as $p_{m}\simeq 0$ .

The Field Eqns. (5) and (6) becomes Einstein Field Eqns. for FLRW space-time
when $f(R)=R$. The additional expressions appearing in these Eqns. on the RHS
are present due to non-linear function f(R). We associate them to curvature
pressure $p_{k}$ and curvature energy $\rho_{k}$ as follows:

${}\rho_{k}=\frac{1}{2}(f-Rf^{\prime})-3\frac{\dot{a}}{a}\dot{R}f^{\prime\prime}$
(7)

and

${}p_{k}=3\frac{\dot{a}}{a}f^{\prime\prime}\dot{R}+f^{\prime\prime\prime}\dot{R}^{2}+f^{\prime\prime}\ddot{R}+\frac{Rf^{\prime}}{2}-\frac{f}{2}.$
(8)

More specifically, we can express Field Eqns. (5) and (6) in terms of the
Hubble parameter $H=\frac{\dot{a}}{a}$ and its derivative $\dot{H}$ as
follows:

${}3H^{2}=\frac{\rho_{m}+\rho_{k}}{f^{\prime}},$ (9)

and

${}2\dot{H}+3H^{2}=-\frac{p_{k}}{f^{\prime}}.$ (10)

The purpose of this work is to form a cosmological model of the universe which
fits best on the observational ground. At present, we have two types of data
sets which are described as follows:

  * •

Hubble Observational Data Set (OHD): A data set of $46$ observational points
of Hubble parameters in the redshift range $0\leq z\leq 2.36$, which are
procured from cosmic chronometric technique (CCT) ref27 .

  * •

Pantheon data Set: A data set of $1048$ SNIa apparent magnitude measurements
known as the Pantheon compilation ref28  includes $276$ SNIa data set of the
Pan-STARRS1 Medium Deep Survey in the redshift range $0.03<z<0.65$ along with
SNIa data set from SDSS, SNLS and low-$z$ HST samples.

Now, we express the field Equations (9) and (10) in terms of derivatives of
the $H$ w.e.t. redshift $z$. For this we use the following transformation
formula for redshift $1+z=\frac{a_{0}}{a}$, where $a_{0}$ is the present value
of $a$. The following field equations are obtained as per our requirement:

$3f^{\prime}H^{2}+0.5\bigg{(}f+6f^{\prime}(2H^{2}-(1+z)H\frac{dH}{dz}\bigg{)}+18(1+z)f^{\prime\prime}H^{2}\bigg{(}3H\frac{dH}{dz}-(1+z)\bigg{(}\frac{dH}{dz}\bigg{)}^{2}-(1+z)H\frac{d^{2}H}{dz^{2}}\bigg{)}=\rho_{m},$
(11)

$\frac{d^{3}H}{dz^{3}}=\frac{3Hf^{\prime}-0.5f-6f^{\prime}H^{2}}{6(1+z)^{3}f^{\prime\prime}H^{3}}+\frac{dH}{dz}\bigg{(}\frac{f^{\prime}+36f^{\prime\prime}H^{2}-6(1+z)f^{\prime\prime}H^{2}}{6(1+z)^{2}f^{\prime\prime}H^{2}}\bigg{)}+\left(\frac{dH}{dz}\right)^{2}\bigg{(}\frac{1}{(1+z)H}+\frac{54f^{\prime\prime\prime}H}{(1+z)f^{\prime\prime}}\bigg{)}+\\\
\left(\frac{dH}{dz}\right)^{3}\bigg{(}\frac{-36f^{\prime\prime\prime}}{f^{\prime\prime}}+\frac{1}{H^{2}}\bigg{)}+\bigg{(}\frac{dH}{dz}\bigg{)}^{4}\bigg{(}\frac{6(1+z)f^{\prime\prime\prime}}{f^{\prime\prime}H}\bigg{)}+\frac{d^{2}H}{dz^{2}}\bigg{(}\frac{-3-(1+z)}{(1+z)^{2}}\bigg{)}\\\
+\bigg{(}\frac{d^{2}H}{dz^{2}}\bigg{)}^{2}\bigg{(}\frac{6(1+z)f^{\prime\prime\prime}H}{f^{\prime\prime}}\bigg{)}+\frac{dH}{dz}\frac{d^{2}H}{dz^{2}}\bigg{(}\frac{-36f^{\prime\prime\prime}H^{2}+4f^{\prime\prime}}{f^{\prime\prime}H}\bigg{)}+\bigg{(}\frac{dH}{dz}\bigg{)}^{2}\frac{d^{2}H}{dz^{2}}\bigg{(}\frac{6(1+z)f^{\prime\prime\prime}}{f^{\prime\prime}}\bigg{)},$
(12)

In obtaining Eqns. (11) and (12), the following expressions in terms of $R$,
$\dot{R}$, $H$ and $\frac{dH}{dz}$ are used:

$\displaystyle R$ $\displaystyle=-6(\dot{H}+2H^{2})$ (13)
$\displaystyle=6H\left(\frac{dH}{dz}(1+z)-2H\right),$ (14)
$\displaystyle\dot{R}$ $\displaystyle=-6(\ddot{H}+4H\dot{H})$ (15)
$\displaystyle=-6\left(-3H^{2}(1+z)\frac{dH}{dz}+(1+z)^{2}H\left(\frac{dH}{dz}\right)^{2}+H^{2}(1+z)^{2}\frac{d^{2}H}{dz^{2}}\right),$
(16)

and

$\ddot{R}=-24(1+z)^{2}H^{2}\bigg{(}\frac{dH}{dz}\bigg{)}^{2}\\\
-18(1+z)H^{3}\frac{dH}{dz}-6(1+z)^{2}H^{3}\frac{d^{2}H}{dz^{2}}-24(1+z)^{3}H^{2}\bigg{(}\frac{dH}{dz}\bigg{)}\bigg{(}\frac{d^{2}H}{dz^{2}}\bigg{)}+6(1+z)^{3}H\bigg{(}\frac{dH}{dz}\bigg{)}^{3}+6(1+z)^{3}H^{3}\frac{d^{3}H}{dz^{3}}.$
(17)

##  III  Numerical solutions and the Dynamics of the Model

As stated earlier, we want to develop a model which provides an acceleration
in the universe at present. For this purpose, we consider a particular form
$f(R)=R+\alpha R^{2}+\beta R^{3}$. Earlier, some authors ref15 ; ref16  have
used second order non-linear Ricci scalar to develope an accelerating universe
model. So, we got motivation from the earlier work ref19 ; ref20  who have
used $f(R)$-gravity in developing their model. So, $f^{\prime}(R)=1+2\alpha
R+3\beta R^{2}$, $f^{\prime\prime}(R)=2\alpha+6\beta R$ and
$f^{\prime\prime\prime}(R)=6\beta$.

###  III.0.1  The deceleration and jerk parameters

The deceleration parameter $q$ and jerk parameter $j$ in terms of $a$ and $H$
are defined as:

$\displaystyle q$ $\displaystyle=-\frac{\ddot{a}}{aH^{2}}$
$\displaystyle=-1+(1+z)\frac{1}{H}\frac{dH}{dz},$

and

$\displaystyle j$ $\displaystyle=\frac{\dddot{a}}{aH^{3}}$
$\displaystyle=1+3\frac{\dot{H}}{H^{2}}+\frac{\ddot{H}}{H^{3}}$
$\displaystyle=1-2\frac{(1+z)}{H}\frac{dH}{dz}+\frac{(1+z)^{2}}{H^{2}}\left(\frac{dH}{dz}\right)^{2}+\frac{(1+z)^{2}}{H}\frac{d^{2}H}{dz^{2}}.$

The present values of the deceleration parameter $q_{0}=-0.55$ and jerk
parameter $j_{0}=1$ of the obtained model is consistent with $\Lambda$CDM
model. Planck’s latest survey ref27  provides us the present value of Hubble
constant $H_{0}\simeq 68\,Mpc/Km/sec\simeq 0.07\,Gyr^{-1}$. From these values,
we can estimate the present values of first and second derivatives of Hubble
parameter as $(\frac{dH}{dz})_{0}=0.0315$ and
$(\frac{d^{2}H}{dz^{2}})_{0}=0.48825$.

The energy parameter $\Omega_{m}$ is defined as
$\Omega_{m}=\frac{\rho_{m}}{\rho_{c}}$, where $\rho_{c}=\frac{8\pi G}{3H^{2}}$
is the critical density. Its present value is given as $\rho_{c0}=1.88\times
10^{-29}h_{0}^{2}\,gm/cm^{3}$, where $0.5<h_{0}<1$. EoS-parameter for
curvature inspired pressure and density may be defined as
$\omega_{k}=\frac{p_{k}}{{\rho}_{k}}$. Since $\rho_{c0}$ is negligibly small,
so we shall make certain corrections as $\rho_{m0}=\Omega_{m0}$ for further
calculations. If we put the present values of $H$ and its derivatives in Eqns.
(7),(8) and (11), we get the following relations amongst the parameters
$\alpha,\beta,\Omega_{m0}$ and $\omega_{k0}$ as:

${}\beta=261.437+1.55647~{}\alpha-17784.8~{}\Omega_{m0}$ (18)
${}\alpha=\frac{22.5937+32.2766~{}\omega_{k0}-952.293~{}\Omega_{m0}-2195.69~{}\omega_{k}~{}\Omega_{m0}}{0.7+\omega_{k0}}$
(19)

$\begin{array}[]{c@{\hspace{0.1in}}cc}\includegraphics[width=216.81pt,height=180.67499pt]{qz}\hfil\hskip
7.22743pt&\includegraphics[width=216.81pt,height=180.67499pt]{jz}\\\
\mbox{(}a)\hfil\hskip 7.22743pt&\mbox{(}b)\end{array}$

Figure 1: The plots of $q-z$ and $j-z$.

In Fig. 2, we see that the deceleration parameter $q$ is negative at present.
This shows that we are living in an accelerating universe. The red shift
transition is shifted from positive to negative values at $z_{tr}\simeq
0.694069$, which shows that the universe transforms from decelerating mode in
the past to accelerating mode in late times. The jerk parameter approaches to
$\Lambda$CDM because $j\rightarrow 1$ as $z\rightarrow 0$.

###  III.0.2  Statistical Analysis of the parameter

As stated earlier in Section(II), we have Hubble parameters data set of $46$
points indicated as $H_{obs}$. Since Eqn. (12) is a non-linear third order
differential equation in derivatives of $H$ and it is not possible to find its
analytical solution, therefore, we go for getting its numerical solution by
putting initial values of first and second derivatives of $H$. We call the
corresponding values of Hubble parameter obtained on the basis of Eqn. (12) as
$H_{th}$.

Now, we use Chi-square test to get better approximations of the parameters
$\Omega_{m0}$ and $\omega_{k0}$. It is defined as

$\chi^{2}=\sum\limits_{i=1}^{46}\frac{[H_{obs}(z_{i})-H_{th}(z_{i})]^{2}}{err{(z_{i})}^{2}},$
(20)

where $err{(z_{i})}$ stands for the standard error in OHD.

The minimum value of Chi-square is found to be $\chi^{2}_{min}=19.348$ for
$\omega_{k0}=-0.44$ and $\Omega_{m0}=0.25$. From the eqns. (18) and (19), we
get $\alpha=1.0027$ and $\beta=-5074.44$. The process is shown in the Fig. 2.
in which blue dot is the estimated point $(-0.44,0.25,19.348)$ for minimum
Chi-square.

$\begin{array}[]{c@{\hspace{0.1in}}c}\includegraphics[width=216.81pt,height=180.67499pt]{chi}\hfil\hskip
7.22743pt&\includegraphics[width=216.81pt,height=180.67499pt]{rgn}\\\
\mbox{(}a)\hfil\hskip 7.22743pt&\mbox{(}b)\end{array}$

Figure 2: Blue and Red dots indicate the present values of $\Omega_{m0}$ and
$\omega_{k0}$ of our model in the Chi-square and Region plots respectively.

###  III.0.3  SNIa Supernova distance modulus and apparent magnitude:

The luminosity distance of any object is defined as:

$D_{l}(z)=(1+z)H_{0}\int_{0}^{z}{\frac{1}{H(z*)}dz*},$

and the distance modulus $\mu(z)$ of a luminous object is correlated to the
luminosity distance $D_{l}(z)$ through the following equation as:

$\mu(z)=m-M=5LogD_{l}(z)+\mu_{0},$ (21)

where $m$ and $M$ are the apparent and absolute magnitude of the object and
$\mu_{0}=25+5Log\left\\{\frac{c}{H_{0}}\right\\}$.

By carrying numerical integration of tabulated values of the Hubble parameter,
we have been able to find a plot $\mu(z)-z$ for our model. We have an error
bar plot for $620$ supernova data sets which include $580$ SNIa union
compilation plus $40$ binned data form ref29 ; ref21 . It is interesting to
see that our theoretical plot coincides with the corresponding $\Lambda$CDM
plot, which is displayed in Fig 2b.

We have presented the Hubble Error and the distance modulus plots versus red
shift in Fig.2. Our plots are matching with graph prepared on the basis of
$\Lambda$CDM model, which shows that these are in a good agreement with our
theoretical and observed results.

$\begin{array}[]{c@{\hspace{0.1in}}c}\includegraphics[width=216.81pt,height=180.67499pt]{er-
hz}\hfil\hskip
7.22743pt&\includegraphics[width=216.81pt,height=180.67499pt]{er-muz}\\\
\mbox{(}a)\hfil\hskip 7.22743pt&\mbox{(}b)\end{array}$

Figure 3: The blue color error bars of observational data $H(z)$ and $SNIa$ in
the plots indicate the deviations of the model from $\Lambda CDM$. Red and
black lines represent our obtained model and the $\Lambda CDM$ respectively.

###  III.0.4  EoS-parameter and Matter energy density

Equation of state parameter $\omega_{k}$ for curvature dominated energy is
defined as $\omega_{k}=\frac{p_{k}}{{\rho}_{k}}$, where $\rho_{k}$ and $p_{k}$
are given by Eqns. (7) and (8). In $\omega_{k}\sim z$ plot, $\omega_{k}$ is
negative and is equal to $-0.44$ at present (see Fig. 3). The curvature
dominated pressure $p_{k}$ becomes negative from positive at the redshift
transition $z\simeq 1.65$. Therefore, we say that the accelerated dark energy
model is developed due to negative pressure.

$\begin{array}[]{c@{\hspace{0.1in}}cc}\includegraphics[width=216.81pt,height=180.67499pt]{rho-m}\hfil\hskip
7.22743pt&\includegraphics[width=216.81pt,height=180.67499pt]{wz}\\\
\mbox{(}a)\hfil\hskip 7.22743pt&\mbox{(}b)\end{array}$

Figure 4: The plots of $\rho_{m}-z$ and $\omega_{k}-z$ .

The numerical values of $H$, $\frac{dH}{dz}$ and $\frac{d^{2}H}{dz^{2}}$ for
redshifts $z$ are calculated from (12) and then using (11), we evaluate matter
density $\rho_{m}$. The related plot of $\rho_{m}\sim z$ can be seen in Fig.
4. As per observations, the value of matter energy density should be
$\rho_{m0}=\Omega_{m0}\times\rho_{c0}$ at present. So putting
$\Omega_{m0}\simeq 0.3$, $\rho_{c0}=1.88\times 10^{-29}\times
h_{0}^{2}\,gm/cm^{3}$ and $h_{0}\simeq 0.68$ we get $\rho_{m0}\simeq
0.27744\times 10^{-29}$, $gm/cm^{3}$, which is consistent with recent
observations.

###  III.0.5 State finder analysis

In recent years, various observations have certified that our universe is
accelerating at the present epoch, which caused a huge list of various DE-
models have been proposed so far. Based on EoS-parameter $\omega_{de}$ for
dark energy, the most eminent models $\Lambda$CDM ($\omega_{de}$=1), $QCDM$
($\omega_{de}\geq-1/3\leq-1$), $KCDM$ ($\omega_{de}$ is negative and
variable), Chaplygin gas energy model and Brain world model etc. are proposed
so far. The references related to these models can be seen in ref30 ; ref31 .
In the year $2003$, Sahni et al. ref25 ; ref26  proposed an analysis to
examine the different category of established dark energy models. This
analysis is based on deceleration parameter $q$ and jerk parameter $r$ as
defined earlier. They have introduced a new parameter $s$, which is defined as
$s=\frac{r-1}{3(q-1/2)},\,q\neq 1/2$, and focused on two-dimensional $s\sim r$
and $q\sim r$ plots. We have presented the plots for our model (see Fig. 5).
It is interesting to see that our model approaches to $\Lambda$CDM from both
ends i.e. from Chaplygin gas to $\Lambda$CDM and from quintessence to
$\Lambda$CDM in $s\sim r$ plot. In $q\sim r$ plot, our model passes nearby
$\Lambda$CDM. So we can say that our model lies near $\Lambda$CDM, and it fits
well on observational ground.

$\begin{array}[]{c@{\hspace{0.1in}}c}\includegraphics[width=216.81pt,height=180.67499pt]{s-r}\hfil\hskip
7.22743pt&\includegraphics[width=216.81pt,height=180.67499pt]{q-r}\\\
\mbox{(}a)\hfil\hskip 7.22743pt&\mbox{(}b)\end{array}$

Figure 5: The plots of $s-r$ and $q-r$ for state finder analysis.

##  IV conclusion

In this work, we have probed a cosmological model with the perfect fluid
filled universe in $f(R)$-gravity by taking $R+aR^{2}+{b}{R^{3}}$ as a
particular form of $f(R)$ function. The field equations (EFE) in
$f(R)$-gravity are solved for a flat FLRW spatially homogeneous and isotropic
space-time. The terms which arise due to non-linear $f(R)$ function are
shifted to the RHS of the field equations that are treated as energy due to
curvature inspired geometry. As a result, it produces acceleration and anti-
gravitating negative pressure in the universe. The plots show that our
theoretical plots fit well with the observational Hubble and Pantheon data.
Moreover, the plots match at par with the $\Lambda$CDM model. The $q-z$ plot
describes an accelerating universe at the present epoch. The transition
redshift $z$ for our model is obtained as $z_{t}\simeq 0.694069$, which is in
good agreement with $\Lambda$CDM. We have also carried out the state finder
analysis for our model, which confirms that our model lies near $\Lambda$CDM,
and it fits well on the observational ground.

## References

  * (1) S. Perlmutter et al., Nature 391, 51 (1998)
  * (2) S. Perlmutter et al., Astrophys. J., 517, 5 (1999)
  * (3) A. G. Riess et al., Astron. J., 116, 1009 (1998)
  * (4) J. P. Ostriker and P. J. Steinhardt, Nature 377, 600 (1995), and references therein.
  * (5) M. S. Turner, G. Steigman, and L. Krauss, Phys. Rev. Lett., 52, 2090 (1984)
  * (6) D. N. Spergel et al., [WMAP Collaboration], Astrophys. J. Suppl. Ser., 170, (2007) 377
  * (7) M. Tegmark et al., [SDSS Collaboration], Phys. Rev. D, 69 , 103501 (2004)
  * (8) E. J. Copeland, M. Sami and S. Tsujikawa, Int. J. Mod. Phys. D, 15, 1753 (2006)
  * (9) O. Gron and S. Hervik, Einstien’s General Theory of Relativity With Modern Application in Cosmology, (New York: Springer) (2007)
  * (10) K. Abazajian et al., [SDSS Collaboration], Astron. J., 128, 502 (2004)
  * (11) V. Sahni and A. A. Starobinsky, Int. J. Mod. Phys. D, 9, 373 (2000)
  * (12) S. Weinberg, Rev. Mod. Phys. 61, 1 ( 1989)
  * (13) P. Steinhardt, L. Wang and I. Zlatev, Phys. Rev. D, 59, 123504 (1999)
  * (14) V. B. Johri, Phys. Rev. D, 63 103504 (2001)
  * (15) S. Nojiri and S. D. Odintsov, Phys. Rev. D, 68 123512 (2003)
  * (16) A. A. Starobinsky, JETP Lett., 86 157 (2007)
  * (17) T. P. Sotiriou and S. Liberati, J. Phys. Conf. Ser.,  68 012022 (2007)
  * (18) T. P. Sotiriou and S. Liberati, Annals Phys., 322 935 (2007)
  * (19) S. K. Srivastava and Phys. Lett, B, 648 (2007) 119 (2007)
  * (20) A. Mukherjee and N. Banerjee, Astrophys. Space Sci., 352 893 (2014)
  * (21) G. C. Samantha and N. Godani, Ind. J. Phys., 94 1303 (2020)
  * (22) Ciprian A. Sporea, arXiv:1403.3852v2 [gr-qc] (2014)
  * (23) M. Mali and S. Shankaranarayanan, Nuclear Phys. B, 937, 422 (2018)
  * (24) S. Capozziello, Int. J. Mod. Phys. D, 11, 483 (2002)
  * (25) Jonathan D. Evans, Lisa M. H. Hall, and Philippe Caillol, Phys. Rev. D, 77, 083514 (2008)
  * (26) J. D. Barrow and S. Cotsakis, Phys. Lett. B, 214 515518, (1988)
  * (27) P. A. R. Ade et al., (Planck 2015 Collaboration), A. A., 594 A20 (2016)
  * (28) J. D. Barrow and S. Cotsakis, Phys. Lett. B, 214 515 (1988)
  * (29) P. Biswas, P. Roy and R. Biswas, arXiv: 1908.00408 [gr-qc] (2019)
  * (30) D. M. Scolnic et al., Astrophys. J, 859 101 (2018)
  * (31) P. A. R. Ade, et al., Planck $\Lambda$ Collaboration, Astron. Astrophys., 594 A13 (2016)
  * (32) N. Suzuki et al., Astrophys. J., 746 (2012) 85
  * (33) V. Sahni  et al., JETP Lett., 77 201 (2003)
  * (34) U. Alam et al., Mon. Not. Roy. Astron. Soc., 344 1057 (2003)


that are not included in the main text. The figure captions contain further
descriptions of the results.

Synthetic Benchmarks. In the synthetic datasets, where all features are
independent, FDT, RF, BART, BNN, GAM, BRR and BL perform better and more
stable than they do in the real datasets where there’s feature correlation.
The better performance of FDT compared to RF-impurity and RF-knockoff
illustrates the advantage of the proposed integrated partial derivative metric
for variable selection. For the synthetic-continuous and synthetic-mixture
cases, FDT has higher AUROC scores across most scenarios, especially when data
are generated having high complexity with quickly-varying local fluctuations
(rbf, matern32). Moreover, all 11 methods perform only moderately well in
complex data settings. The two tree-based methods, RF and BART also have high
AUROC scores across scenarios, since the tree-based methods naturally rank by
how well the features improve the purity of the node. Note that under low
dimension case ($d=25$), BKMR is comparable to FDT when $f_{0}\in\Hsc$
(linear, rbf, matern32). However, when it comes to medium- or relatively high-
dimension settings ($d=50,100$), BKMR produces low AUROC scores due to
suffering from the issue of curse of dimensionality [71]. RFNN, also a kernel-
based method, has similar trend as BKMR. Finally, BAKR performs consistently
poorly and has lowest AUROC scores in relatively low-dimension setting
($d=25,50$). Linear models (GAM, BRR and BL) achieve comparable or superior
performance under the linear data setting. However, for more complicated data
generation functions, BRR and BL consistently perform poorly with low AUROC
scores.

Socio-economic and Healthcare Datasets In the adult, heart and mi cases, where
the features are correlated, the performances of all 11 methods are worse than
in the synthetic-mixture and synthetic-continuous cases (where the features
are independent). Their performance tends to saturate earlier and are less
stable with respect to the sample size. In relatively low-dimension settings
($d=25,50$), the standard methods such as BART and BNN have higher AUROC
scores than FDT. However, when the dimension is higher ($d=100$), FDT
consistently performs better.

Figure 6: AUROC scores for synthetic-mixture data. FDT generally outperforms
other methods in most of the data settings in relatively higher dimension
($d=50,100$). Knockoff with random forest statistics produce lower AUROC
scores than in synthetic-continuous, even in linear data settings. Additive
models BRR, BL and GAM have mediocre scores under the nonlinear settings. Some
model (e.g., GAM) reports missing result in $n>p$ setting due to the
restriction of their implementations.

Figure 7: AUROC scores for synthetic-continuous data. FDT generally
outperforms other methods in most of the data settings, with BKMR as the
comparable one when $d=25$. However, BKMR performs poorly in higher dimension.
Tree-based methods RF, BART and Knockoff with random forest statistics have
high AUROC scores. The performance of BNN is between tree-based methods and
additive models while BAKR performs poorly consistently. Additive models BRR,
BL and GAM have mediocre scores under the nonlinear settings. Some model
(e.g., GAM) reports missing result in $n>p$ setting due to the restriction in
their implementations.

Figure 8: AUROC scores for adult data. In relatively low-dimension settings
($d=25,50$), the standard methods such as BART and BNN have higher AUROC
scores than FDT. However, when the dimension is higher ($d=100$), FDT performs
better consistently. Some model (e.g., GAM) reports missing result in $n>p$
setting due to the restriction in their implementations.

Figure 9: AUROC scores for heart data. In relatively low-dimension settings
($d=25,50$), the standard methods such as BART and BNN have higher AUROC
scores than FDT. However, when the dimension is higher ($d=100$), FDT performs
better consistently. Some model (e.g., GAM) reports missing result in $n>p$
setting due to the restriction in their implementations.

Figure 10: AUROC scores for mi data. In low-dimension setting ($d=25$), the
standard methods such as BART and BNN have higher AUROC scores than FDT.
However, when the dimension is higher ($d=50,100$), FDT and GAM perform better
consistently. Some model (e.g., GAM) reports missing result in $n>p$ setting
due to the restriction in their implementations.

Figure 11: Testing MSE for synthetic-mixture data. FDT generally performs
better or competitively with baselines, except in the linear case where BL
unsurprisingly does best. BKMR consistently performs worse than other methods,
except in the low data, high dimension setting when BAKR performs worst. Some
model (e.g., GAM) reports missing result in $n>p$ setting due to the
restriction in their implementations. Notice that this dataset contains a
setting $n=p$, which can lead to the double descent phenonmenon for some
random-feature-based models [17].

Figure 12: Testing MSE for synthetic-continuous data. A method will not be
shown if they share the model fit with another method (RF-impurity and RF-
knockoff), or if the method does not produce valid result due to small sample
size (GAM).

Figure 13: Testing MSE for adult data. A method will not be shown if they
share the model fit with another method (RF-impurity and RF-knockoff), or if
the method does not produce valid result due to small sample size (GAM).

Figure 14: Testing MSE for heart data. A method will not be shown if they
share the model fit with another method (RF-impurity and RF-knockoff), or if
the method does not produce valid result due to small sample size (GAM).
Notice that this dataset contains a setting $n=p$, which can lead to the
double descent phenonmenon for some random-feature-based models [17].

Figure 15: Testing MSE for mi data. A method will not be shown if they share
the model fit with another method (RF-impurity and RF-knockoff), or if the
method does not produce valid result due to small sample size (GAM).

##  Appendix I Additional Experiments: Regularization Path for Bangladesh
birth cohort study

We propose a way to visualize the selection path that incorporates the
uncertainty of variable importance scores. Specifically, we consider the
posterior survival function $S(s)=P(\psi_{j}>s),j=1,\ldots,d$ for increasing
$s$ starting from $0$. Larger value of $S(s)$ indicates larger probability of
that certain feature being relevant. This is analogous to the regularization
path under the LASSO method. However, our approach incorporates posterior
uncertainty, and does not require repeated model fitting at different levels
of regularization strength [49].

We apply this to Bangladesh birth cohort study [41] (a well-established
dataset in the environmental health literature), where we fit models to learn
the association between infant’s neural development scores and key
environmental factors such as hospital location (clinic), sex (sex), levels of
macro nutrient intake (prot, fat, carb, fib, ash) and levels of measured
concentration of environmental toxins in body fluids (as$\\_$ln, mn$\\_$ln,
pb$\\_$ln), while controlling for other socio-economic and biological factors
(family income, parent education levels, etc). In general, the level of macro-
nutrient intake (in particular fiber and protein) indicates a child’s general
nutrition status (i.e., whether he/she is eating well), and is known to be
positively associated with neural development. On the other hand, the existing
studies in the Bangladesh population have established a neurotoxic effect
between arsenic exposure (i.e., as$\\_$ln), through drink water) on the early-
stage cognitive development [30], as well as weak but significant effect of
the joint mixture of other environmental toxins (manganese (mn$\\_$ln) and
lead (pb$\\_$ln))) [28, 70]. Furthermore, due the fact that the model has
already controlled for biological and socio-economic confounding factors, non-
nutrient-related factors such as hospital location and sex should not have a
significant effect on the children’s neural development status.

The variable selection result is shown in Figure 16, where we plot the
posterior survival function $P(\psi_{j}>s)$ for $s\in(0,1)$, and compare it to
the survival function under Bayesian Approximate Kernel Regression (BAKR),
Bayesian Ridge Regression (BRR), Bayesian LASSO (BL), and also the frequentist
LASSO regularization path under the GAM model. We normalized all variable
importance scores within the range $(0,1)$. As a result, the variable
selection performance is indicated by the relative magnitudes of the area
under the curve for each variable (and not by the absolute magnitude due to
the normalization).

As shown in Figure 16, the top variables selected by our method (FDT)
correspond well with existing conclusions in the literature: it correctly
picked up the larger impact of macro-nutrients (in particular, fibre, fat and
protein) and smaller but still significant effects of environmental toxins
(arsenic, manganese and lead), also notice that it ranked known non-causal
factors such as hospital location and sex to be the lowest. In comparison, the
linear methods (GAM, BRR and BL) all incorrectly reported high effect from
hospital location on children’s neural developement outcome (likely due to
their restrict model form), while the nonlinear model (BAKR, based on RBF
kernel) did not properly pick up the effect of environmental toxins.

Figure 16: Regularization path for Bangladesh birth cohort study. The top
variables selected by our method (FDT) correspond well with established
toxicology pathways in the literature.