# GrASPE: Graph based Multimodal Fusion for Robot Navigation in Unstructured
Outdoor Environments

Kasun Weerakoon, Adarsh Jagan Sathyamoorthy, Jing Liang, Tianrui Guan, Utsav
Patel, and Dinesh Manocha
Supplemental version including Tech Report, and Video at
http://gamma.umd.edu/graspe/

###### Abstract

We present a novel trajectory traversability estimation and planning algorithm
for robot navigation in complex outdoor environments. We incorporate
multimodal sensory inputs from an RGB camera, 3D LiDAR, and robot’s odometry
sensor to train a prediction model to estimate candidate trajectories’ success
probabilities based on partially reliable multi-modal sensor observations. We
encode high-dimensional multi-modal sensory inputs to low-dimensional feature
vectors using encoder networks and represent them as a connected graph to
train an attention-based Graph Neural Network (GNN) model to predict
trajectory success probabilities. We further analyze the image and point cloud
data separately to quantify sensor reliability to augment the weights of the
feature graph representation used in our GNN. During runtime, our model
utilizes multi-sensor inputs to predict the success probabilities of the
trajectories generated by a local planner to avoid potential collisions and
failures. Our algorithm demonstrates robust predictions when one or more
sensor modalities are unreliable or unavailable in complex outdoor
environments. We evaluate our algorithm’s navigation performance using a Spot
robot in real-world outdoor environments. We observe an increase of x% in
terms of success rate and a y% decrease in false positive estimations compared
to the state-of-the-art navigation methods.

##  I Introduction

Mobile robots have increasingly been utilized in numerous outdoor applications
such as delivery[1], agriculture[2], surveillance[3], exploration[4], rescue
missions[5], etc. These applications need the ability for the robots to
navigate challenging outdoor environmental conditions such as inclement
weather, low lighting, cluttered vegetation, etc. In this work, we consider
such environments as unstructured outdoor environments.

The robots’ perception could encounter noise, occlusions, or other mode of
error or failures when navigating in such unstructured outdoor environments.
Especially, cameras undergo motion blur, low lighting and occlusions [6, 7]
while LiDAR point cloud experience heavy distortions in cluttered vegetation.
A key issue is to develop methods that can perform reliable perception and
planning computations by taking into account such sensor uncertainties.

Figure 1: Spot robot navigation in unstructured outdoor terrains using GrASPE.

In terms of robot perception, deep learning (DL) techniques have been proposed
for multi-modal sensor fusion. These methods combine distinct feature
representations from different sensor modalities such as camera, IMU, and
lidar [8, 9, 10, 11] or perform only camera-lidar fusion [12, 13]. Some of
these methods use Graph Neural Networks (GNN) [14], which have been used for
multi-modal sensor fusion in indoor scenes [15, 16]. A key challenge is to
develop reliable methods for outdoor scenes. Moreover, the existing DL and GNN
methods do not consider the variable reliability of the cameras and lidars in
highly unstructured environments.

In terms of planning in unstructured outdoor environments, the robot must
compute trajectories that avoid collisions and prevent the robot from getting
stuck (e.g., in vegetation such as bushes, vines, etc). To this end, several
recent works have incorporated anomaly detection algorithms to identify
collisions or failures during navigation as anomalies based on multi-sensory
observations [17, 18, 19]. Such algorithms can be trained using simple
positively and negatively labeled samples which can be created trivially as
opposed to the extensive labeling required for training supervised learning
methods [20, 21, 22]. These positive and negative samples are typically
corresponding to traversable and non-trversabale observations (e.g. due to
possible collisions or other mode of failures). However, since these
observations do not account for potential future navigation failures and
collisions, such methods could lead to catastrophic accidents during
navigation.

To deal with possible future navigation failures due to sensor uncertainties,
proactive anomaly detection methods have been utilized to estimate the
probability of future navigation failure using predicted trajectories, and
multi-sensor observations [7, 23]. However, robots navigating in unstructured
outdoor environments frequently experience varying lighting conditions, and
partial or complete sensor occlusions that lead to distorted sensory
observations.

Main Contributions: We present GrASPE(Graph Atenntion based Sensor fusion for
Path Evaluation), a novel trajectory traversability estimation and planning
algorithm for legged robot navigation in unstructured outdoor environments. We
incorporate multi-modal observations from an RGB camera, 3D LiDAR, and robot
odometry to train a prediction model to estimate candidate trajectories’
success probability for navigation. Our model learns the correlation between
multi-sensor data under unstructured outdoor conditions where the camera
undergoes occlusions, motion-blur, and low-lights while the LiDAR experience
heavy point cloud distortions. We further evaluate the reliability of camera
and LiDAR inputs to achieve reliability-aware multi-modal fusion in our
prediction model via a graph attention formulation. The key contributions of
our approach include:

  * •

A novel graph-based prediction model that can estimate a given trajectory’s
probability of avoiding navigation failures (i.e. success probabilities) such
as collisions and getting stuck based on partially reliable multi-modal sensor
observations in unstructured outdoor environments. The prediction model is
capable of estimating the success probabilities of a trajectory even when the
sensory observations are distorted due to varying environmental conditions.

  * •

Novel reliability estimation methods to quantify the usefulness of image and
3D point cloud data for sensor fusion. Our model incorporates feature encoding
networks to project high-dimensional sensor inputs into feature vectors. The
encoded feature vectors are represented as a connected undirected graph to
train an attention-based GNN for sensor fusion. The proposed reliability
estimations are employed to augment the weights of this graph representation.
In other words, the graph representation is modified with improved
reliability-awareness.

  * •

A Dynamic Window Approach (DWA) based planner to perform navigation on
unstructured outdoor terrains while avoiding potential failures and collisions
using our prediction model. The overall pipeline generates dynamically
feasible, collision free and traversable actions to navigate a legged robot.
We implement our algorithm on a Boston Dynamics Spot robot to evaluate its
performance in real-world unstructured outdoor scenarios. Our method results
in an x% improvement in terms of success rate compared to state-of-the-art
methods.

##  II Related Work

In this section, we discuss the existing work on multi-modal sensor fusion and
anomaly detection in the context of robot navigation

###  II-A Multi-modal Sensor Fusion for Navigation

Robot perception in real world environments can be challenging due to varying
conditions such as lighting, noise, occlusions, motion blur etc. To this end,
prior algorithms incorporate uncertainty modeling and adaptation techniques
for robot navigation using a single sensor modality [24, 25]. However, such
methods perform well only under controlled outdoor settings, where the given
sensor inputs experience limited perturbation in terms of noise from the
environment. Therefore, multi-sensor fusion methods are used to deal with real
environmental conditions and to mitigate unreliable perception from individual
sensors [26, 27, 28]. For instance, different variants of the Kalman Filter
[29] are widely used to fuse Odometry (from wheels or visual/LiDAR
estimations), IMU, and GPS sensor data to compute improved localization for
navigation on slippery and complex terrains [30]. However, these techniques
require the different sensors to provide the same type of feature observations
(e.g., odometry) which limits their applicability to combine sensors that can
provide distinct environmental features (e.g. camera and LiDAR).

Deep learning methods have also been widely used for sensor fusion [26, 28,
31, 32]. Particularly, the camera-lidar fusion is employed in navigation and
simultaneous localization and mapping (SLAM)[33] literature to obtain combined
perception from visual and geometric features [9, 12]. However, these methods
typically assume that both camera and LiDAR perception are reliable in terms
of feature availability. In contrast, our proposed algorithm attempts to deal
with real unstructured outdoor environments where the above assumption is
inapplicable.

###  II-B Anomaly Detection based Prediction Models

Anomaly detection has gained increased attention in recent literature due to
the advancements in Deep Learning[34]. Some methods [20] utilize RGB and Depth
images as multi-modal observations to train a predictive model that can
generate an anomaly mask that reflects ”known” and ”unknown(anomaly)” regions
in the scene. However, such methods assume constant illumination and feature-
rich image observations for accurate predictions. A supervised variational
autoencoder (SVAE) is used [35] for failure identification in unstructured
environments using 2D LiDAR data. This SVAE is used [7] as a LiDAR feature
encoder to identify navigation anomalies proactively in crop fields using
multi-modal fusion. However, the 2D lidar observations are heavily distorted
in cluttered outdoor environments such as tall grass and 3D point-cloud can
provide richer information, particularly in outdoor scenarios. Other
navigation models [36, 21] execute the actions by avoiding undesirable
maneuvers in cluttered environments based on the current visual sensory
observations. However, such vision-based systems lead to erroneous predictions
during camera occlusions and low light conditions. To this end, we propose a
camera, 3D LiDAR, and odometry sensor fusion strategy to achieve better
perception in unstructured outdoor environments for robot navigation.

##  III Background and Problem Formulation

In this section, we state the problem formulation and provide background to
the Graph Neural Networks (GNNs) used in our approach.

###  III-A Problem Formulation

In this paper, we address the problem of legged robot navigation in
unstructured outdoor environments using multimodal sensor fusion. We solve the
problem in two steps: 1.) Trajectory success probability prediction from
multi-modal sensor data 2.) Planning while avoiding unsuccessful trajectories
identified using step 1. We incorporate RGB image, 3D Pointcloud, and robot’s
velocity history in step 1 to train a prediction model to estimate a given
trajectory’s success probability at each time step $t$ as a vector. In this
context, we define success probabilities as the navigation probability of not
encountering collisions, stumbling, or similar modes of failure.

During runtime, a dynamic window approach(DWA) based planner is utilized to
generate velocities to reach a goal. However, instead of executing those
velocities directly to the robot, we calculate the corresponding future
trajectories and feed them to the prediction model along with the multi-sensor
data. The prediction model estimates the candidate trajectory’s success
probability at each time step $t$ as a vector. Only the velocities
corresponding to the trajectories with higher successful probabilities will be
executed for planning. Hence, the prediction model trained in step 1 acts as a
traversability evaluator for the actions generated by the planner in step 2.
The overall system architecture of the prediction model is presented in Fig. 2
and details of our planner are discussed in Section IV-E.

###  III-B Graph Neural Networks (GNNs)

Graph Neural Networks (GNNs) are a class of Neural Networks that can be
operated directly on data represented with graphs. In our work, we represent
the dimension-reduced multi-sensor features and their reliability measures as
a connected graph to learn the correlations between the sensor features. These
connections in the graph are mathematically described by a weighted graph
$\mathcal{G}:=(\mathcal{V},\mathcal{W})$, where $\mathcal{V}={1,2,…,N}$ is a
set of $N$ nodes/vertices associated with elements of the sensors’ feature
vectors. $\mathcal{W}$ is a weighted adjacency matrix with entries
$\mathcal{W}_{i,j}$ representing the strength of the connection (edge) between
feature elements $i$ and $j$.

Figure 2: GrASPE System Architecture for candidate trajectory’s success
probability estimation. We utilize RGB image, 3D point cloud, robot’s velocity
history, and the predicted trajectory of the robot as an image to train a
novel graph based prediction model to estimate the given trajectory’s probably
of success in terms of navigation.

We incorporate Graph Convolution Network (GCN) [37] and Graph Attention
Network [38] operators in our GNN to encode and pay attention to the local
graph structure and node embeddings. GCN is a generalized form of the
convolution in CNNs, where the convolution is operated on irregular or non-
euclidean data structures (such as graphs) to learn features from a node and
its neighbours. We use Graph Attention Network(GAT) [39] to identify the
hidden representations of the graph nodes by attending over its neighbors
using a self-attention strategy.

##  IV Proposed Method

In this section, we present the details of our proposed navigation algorithm.
We fist discuss the components of the prediction model trained for
trajectory’s success probability estimation: 1. Multi-modal data acquisition
and prepossessing, 2. Reduced-dimensional feature representation for multi-
modal data, 3. Camera and LiDAR reliability estimation, 4. Graph-based
prediction model. Finally, we explain the details of our planning algorithm.

###  IV-A Multi-modal Sensors

We design our prediction model (i.e. success probability estimation model)
using the following multi-sensor inputs as observations: RGB image
$I_{t}^{rgb}\in\mathbb{R}^{w\times h\times 3}$, 3D pointcloud
$P_{t}^{lidar}\in\mathbb{R}^{3\times N_{p}}$, robot’s velocity history
$V_{t}=\\{(v,\omega)_{t-T+1},...,(v,\omega)_{t}\\}$ and predicted trajectory
from the robot’s current velocity as an image
$I_{t}^{traj}\in\mathbb{R}^{w\times h\times 1}$. To obtain the predicted
trajectory image $I_{t}^{traj}$, we first extrapolate the robot’s trajectory
for the next $T$ times steps based on the current velocity $(v,\omega)_{t}$.
Then, the resulting trajectory is projected to a blank image using the
homography projection to generate $I_{t}^{traj}$. The RGB image and 3D point
cloud captures visual and geometric features of the environment respectively.
The robot’s velocity history is used to account for the capture the robot’s
recent behavior. Finally, the prediction model incorporates these observations
$O_{t}=\big{[}I_{t}^{rgb},P_{t}^{lidar},V_{t},I_{t}^{traj}\big{]}$ and
estimates the navigability of the predicted trajectory
$\hat{E}_{t}=[p_{t},p_{t+1},..,p_{t+T-1}]$ (i.e. success probability vector).
Details of the observation data collection and labeling are described in
Section V-B.

Figure 3: Image Reliability Estimation: We estimate input image reliability
based on two factors: 1. Overall lighting condition 2. Feature richness. (a)
and (b) demonstrate images of the same scene under two different lighting
conditions. (c) and (d) presents the FAST features extracted from the images
(a) and (b). We observe that the images obtained under low light conditions
have a significantly less number of features. Similarly, the motion blur and
camera occlusions result in images with low reliability measures. (see cite
arxiv version for detailed results).  Figure 4: Point cloud Reliability
Estimation:We estimate point cloud reliability based on the edge and planar
feature availability. (a) Top view of a feature-rich point cloud observation
(b) Top view of a distorted point cloud observed in tall grass environments.
We observe that the 3D point cloud is feature-rich in open and structured
outdoor environments. However, the cluttered vegetation leads to heavy
distortions in the point cloud.  Figure 5: Navigability predictions from our
GrASPE algorithm, GrASPE w/o reliability, GA-Nav[40], and PSPNet[41] under
different environmental conditions. (a) Hanging leaves; (b) Pliable tall
grass; (c) Camera occlusion; (d)Low-light condition. We observe that the
GrASPE algorithm predictions outperform all the other comparison methods under
varying environmental conditions that are critical for robot navigation in
unstructured outdoor environments.

###  IV-B Multimodal Feature Vector Generation

Different sensor modalities capture different environmental features (images
capture visual features, point clouds capture edges and surfaces), and fusing
them is not straightforward. Moreover, raw image and point cloud data
processing is computationally expensive due to the high dimensionality.
Therefore, we process each sensor input using separate feature encoding
networks to obtain dimension-reduced feature vector representations.

####  IV-B1 Visual Feature Extraction

We utilize a ResNet[42] based pipeline to extract image features $f_{img}$
from the camera RGB image $I_{t}^{rgb}$. We incorporate CBAM[43] module
between the network layers to perform spatial and channel attention toward the
important visual features.

####  IV-B2 LiDAR Feature Extraction

Pointnet [44] based network is incorporated to extract pointcloud features
$f_{point}$ from the lidar data $P_{t}^{lidar}$.

####  IV-B3 Velocity Feature Extraction

Velocity vector is fed into a set of linear convolutional layers in VelNet to
obtain the velocity features $f_{vel}$.

####  IV-B4 Trajectory Feature Extraction

The predicted trajectory image $I_{t}^{traj}$ passes through a set of pooling
layers and CNN layers to obtain the trajectory feature vector $f_{traj}$.

###  IV-C Sensor Reliability Estimation

We observe that the RGB images and the point cloud data are unreliable in
certain unstructured outdoor environments due to the inability of visual and
geometric feature detection. Even though the feature encoding networks
proposed in the literature and in Section IV-B1 are capable of feature
extraction from images and point clouds, they do not account for the
reliability of the input sensory data. To this end, we perform reliability
estimation on image and point cloud observations to quantify the reliability
as a scalar value.

We estimate the reliability of the image and point cloud inputs at each time
step quantitatively using classical image and point processing methods that
execute in real-time [45, 46]. We further assume that the instantaneous robot
velocities obtained from its odometry are not significantly affected by the
environment. Therefore, it is considered as reliable during operation.

####  IV-C1 Image Reliability Estimation

We consider that the images captured by the camera have a high reliability if
they have: 1. High image brightness, and 2. Availability of visual features
such as corners, edges etc [47]. See Fig. 3 for a sample reliability
comparison.

To estimate the image brightness, we first convert the input RGB image
$I_{t}^{rgb}$ to a gray-scale image $I_{t}^{gray}$. Then we calculate the Root
Mean Square(RMS) value of the histogram distribution of $I_{t}^{gray}$ as
follows to obtain brightness estimation,

$r_{bright}=RMS(hist(I_{t}^{gray}))$ (1)

Here, $hist$ is the histogram operator and $RMS$ is the root mean square
operator where $RMS(x)=\sqrt{\frac{1}{n}\sum_{i=1}^{n}x_{i}^{2}}$.

To estimate the availability of useful features in the image, we calculate the
number of corner features $n_{c}$ in the input image $I_{t}^{rgb}$ using
FAST(Features from Accelerated Segment Test) algorithm proposed in [48, 45].
We consider the image is feature rich if the $n_{c}$ is higher than a
threshold. Hence, $r_{corners}=\frac{n_{c}}{(w\times h)}$. Here, $w$ and $h$
are the height and width of the input RGB image. The final image reliability
measure is obtained as a normalized scalar value $r_{img}\in[0,1]$ using the
weighted sum of $r_{bright}$ and $r_{corners}$ :

$r_{img}=\alpha_{b}r_{bright}+\alpha_{c}r_{corners}$ (2)

####  IV-C2 Pointcloud Reliability

We observe that the lidar point cloud is heavily distorted in the presence of
unstructured vegetation leading to poor estimation of the surrounding objects’
geometry. Therefore, to evaluate the point cloud reliability, we perform edge
and planar feature extraction to calculate the number of 3D features available
in the point cloud $P_{t}^{lidar}$ at a given time $t$.

Let $l$ be a 3D point in the input point cloud $P_{t}^{lidar}$ and
$\mathcal{M}$ be the set of coordinates of point $l$ acquired from a single
LiDAR scan. Then, we can obtain the point feature evaluation factor $c_{l}$ of
the point $l$ as,

$c_{l}=\frac{1}{||X_{l}||.|\mathcal{M}|}\Big{\lVert}\sum_{k\in\mathcal{M},k\neq
l}\big{(}X_{l}-X_{k}\big{)}\Big{\rVert},$ (3)

where, $X_{k}$ with $k=1,2,..,|\mathcal{M}|$ are the coordinates of the 3D
points in the set $\mathcal{M}$. Zhang et al. [46] demonstrate that the points
with higher and lower $c_{l}$ values obtained from the equation 3 belong to
edge and planar features respectively. Hence, we define two threshold values
$c_{max}$ and $c_{min}$ to extract the edge and planar points using $c_{l}$.
The resulting edge and planar feature point sets can be denoted as:
$\mathcal{S}_{edge}=\\{l|c_{l}\geq c_{max},l\in P_{t}^{lidar}\\}$ and
$\mathcal{S}_{planar}=\\{l|c_{l}\leq c_{min},l\in P_{t}^{lidar}\\}$
respectively.

Then, the point cloud reliability measures, $r_{edge}$ and $r_{planar}$, are
calculated using the ratio between the cardinality of each feature point set
and the input point cloud.

$r_{edge}=\frac{|\mathcal{S}_{edge}|}{|P_{t}^{lidar}|},\quad
r_{planar}=\frac{|\mathcal{S}_{planar}|}{|P_{t}^{lidar}|}.$ (4)

Finally, the point cloud reliability measure $r_{point}$ is derived as a
weighted combination of $r_{edge}$ and $r_{planar}$ using two tunable
parameters $\beta_{e}$ and $\beta_{p}$.

$r_{point}=\beta_{e}r_{edge}+\beta_{p}r_{planar}.$ (5)

Fig. 4 presents an example comparison of $r_{edge}$ and $r_{planar}$ in real
outdoor scenarios.

###  IV-D Graph Attention based Prediction Model

Once the sensor inputs are encoded into feature vectors (Section IV-B), we
concatenate them to obtain a combined vector
$f_{vec}=[f_{img},f_{point},f_{traj},f_{vel}]$. This $f_{vec}\in[0,1]^{N}$ is
utilized to generate the graph representation for fusion.

Let $\mathcal{G}_{t}:=(\mathcal{V},\mathcal{W})_{t}$ be a graph defined in
III-B at time $t$. We consider the graph nodes $\mathcal{V}={1,2,…,N}$ as a
set of $N$ nodes/vertices associated with elements of the feature vector
$f_{vec}$. Then, we define the feature vector $f_{vec}$ as the graph signal
residing on the graph nodes $\mathcal{V}$. i.e. the $i^{th}$ component of the
vector $f_{vec}$ denotes the value at the $i^{th}$ node in $\mathcal{V}$.

The weighted adjacency matrix $\mathcal{W}$ should represent the strength
between the graph nodes (i.e. edge weights) i.e, the connectivity between the
encoded feature elements in $f_{vec}$. Typically, the edge weights are
calculated using a distance measure (absolute difference or L2 norm). However,
we modify this edge weights (especially between the nodes corresponding to the
image and point cloud features) based on the reliability measures (Section
IV-C) as follows.

Let $\mathcal{W}_{i,j}$ be the edge weight between the feature elements $i$
and $j$ in $f_{vec}$.

$\mathcal{W}_{i,j}=\begin{cases}exp^{-\\{\lambda|f_{vec}(i)-f_{vec}(j)|\times
r_{i,j}\\}},&\text{if }i\neq j\\\ 0,&\text{otherwise},\end{cases}$ (6)

where, $\lambda$ is a tunable scalar and $r_{i,j}\in[0,1]$ represents the
reliability measure between the node $i$ and $j$ as follow,

$r_{i,j}=\begin{cases}r_{img},&\text{if }\\{i,j|i\in f_{img},j\notin
f_{point}\\}\\\ r_{point},&\text{if }\\{i,j|i\notin f_{img},j\in
f_{point}\\}\\\ \frac{1}{2}\big{(}r_{point}+r_{img}\big{)},&\text{if
}\\{i,j|i\in f_{img},j\in f_{point}\\}\\\ 1,&\text{otherwise}.\\\ \end{cases}$
(7)

Hence, we refer the derived graph representation $\mathcal{G}_{t}$ of the
feature vector $f_{vec}$ as reliability-aware feature graph. We design a
light-weight Graph Neural Network with a set of Graph Convolution and Graph
Attention layers defined in III-B to predict the navigability estimation
vector $\hat{E}_{t}$ for a given trajectory using the reliability-aware
feature graph $\mathcal{G}_{t}$. The overall prediction pipeline can be
treated as a mapping function $\psi:O_{t}\rightarrow\hat{E}_{t}$.

###  IV-E Outdoor Navigation with DWA

We adapt a modified version of the well-known Dynamic Window Approach (DWA)
[49] to perform navigation while evaluating the navigability of the generated
trajectories using GrASPE.

We represent the robot’s actions as a linear and angular velocity pairs
$(v,\omega)$. Let $V_{s}$ be the space of all the possible robot velocities.
The DWA algorithm considers two velocity space constraints to obtain
dynamically feasible and collision-free velocities: (1). The dynamic window
$V_{d}$ contains the dynamically feasible velocities during the next $\Delta
t$ time window; (2). The admissible velocity space $V_{a}$ includes the
collision-free velocities. The restricted velocity space $V_{r}=V_{s}\cap
V_{d}\cap V_{a}$.

These collisions for $V_{a}$ are generally identified by evaluating the
potential future trajectory collisions using 2D LiDAR scan data. However, the
LiDAR scan is highly distorted and unreliable in cluttered outdoor terrains
such as thin grass and vegetation. Hence, the traditional DWA will identify
traversable cluttered regions as obstacles. To avoid such misclassifications,
we omit the admissible velocity space constraint when the point cloud is
highly unreliable.

$V_{r}=\begin{cases}V_{s}\cap V_{d},&\text{if }r_{point}\leq r_{th}\\\
V_{s}\cap V_{d}\cap V_{a}&\text{otherwise},\end{cases}$ (8)

where $r_{th}$ is a reliability threshold that indicates the point cloud data
is highly distorted. The resulting velocity space $V_{r}$ is utilized to
calculate the optimal velocity pair $(v^{*},\omega^{*})$ by maximizing the
objective function below,

$Q(v,\omega)=\sigma\big{(}\gamma_{1}.heading(v,\omega)+\gamma_{2}.dist(v,\omega)+\gamma_{3}.vel(v,\omega)\big{)},$
(9)

where $heading(.),dist(.)$ and $vel(.)$ are the cost functions [49] to
quantify the robot’s heading towards the goal, distance to the closest
obstacle in the trajectory, and the forward velocity of the robot. $\sigma$ is
a smoothing function and $\gamma_{i},i=1,2,3$ are adjustable weights. The
resulting optimal velocity pair can be denoted as
$(v^{*},\omega^{*})=argmax(Q(v,\omega))$.

We calculate the predicted trajectory for the next $T$ time steps considering
the optimal $(v^{*},\omega^{*})$ as the robot’s immediate next velocity. Then,
the obtained trajectory is fed into our GrASPE algorithm to estimate the
trajectories’ success probabilities $\hat{E}_{t}$. If $min(\hat{E}_{t})\geq
e_{th}$ and $mean(\hat{E}_{t})\geq e_{th}$, the velocity pair
$(v^{*},\omega^{*})$ corresponding to that trajectory is considered as
navigable. Otherwise, the $(v^{*},\omega^{*})$ is removed from velocity space
$V_{r}$ and a new optimal velocity pair is calculated using the objective
function defined in equation 9. Here, $e_{th}$ is an adjustable threshold.

##  V Results and Analysis

We explain our method’s implementation and experimental details on a Spot
robot. Then, we conduct ablation studies and comparisons to highlight benefits
of our algorithm.

###  V-A Implementation

GrASPE algorithm is implemented using PyTorch and PyTorch Geometric(PyG). The
prediction model is trained in a workstation with an Intel Xeon 3.6 GHz
processor and an Nvidia Titan GPU using the data collected from a Spot robot.
The robot is equipped with an Intel NUC 11 (a mini PC with Intel i7 CPU and
NVIDIA RTX 2060 GPU), a VLP16 Velodyne LiDAR and a L515 Intel RealSense LiDAR
Camera.

###  V-B Dataset

RGB image $I_{t}^{rgb}$ and the point cloud $P_{t}^{lidar}$ are obtained from
the aforementioned camera and the LiDAR respectively. The robot’s velocity
history $V_{t}$ is collected from the odometry. The ground truth labels
$E_{t}$ are generated manually using element-wise binary labeling for
simplicity. i.e. the trversable portion of the predicted trajectory is labeled
as 1’s and the non-traversable or possible failure regions are labeled with
0’s. We incorporate a randomized planner during the data collection to
minimize the human intervention. Hence the observations corresponding to the
sucessful navigation scenarios are automatically labeled using $E_{t}$ vectors
with all 1’s. However, the non-traversable or failure cases were manually
labeled by evaluating the observations after the data collection. The final
dataset includes a set of observation and label pairs $(O_{t},E_{t})$ for each
time step $t$. The test and training data sets are collected in University of
Maryland outdoor fields. The training set includes 28721 positive samples and
8463 negative samples whereas the test set contains 4537 positive samples and
1096 negative samples.

###  V-C Evaluations

We use following evaluation metrics to compare our method’s performance
against three state of the art navigation methods: GA-Nav[40], PSPNet[41] and
DWA[49]. We further compare GrASPE model without reliability estimations for
ablation studies. GA-Nav is a multi-modal navigation framework that combines
vision-based terrain segmentation with point cloud based elevation estimation
for outdoor navigation. PSPNet is a vision based segmentation algorithm
trained for outdoor terrain segmentation. DWA is a local planner that utilizes
2D LiDAR scan for obstacle avoidance.

Success Rate \- The number of successful goal reaching attempts (while
avoiding non-pliable regions and collisions) over the total number of trials.

Normalized Trajectory length \- The ratio between the robot’s trajectory
length and the straight-line distance to the goal in both successful and
unsuccessful trajectories.

##  VI Conclusions, Limitations and Future Work

## References

  * [1] R. Limosani, R. Esposito, A. Manzi, G. Teti, F. Cavallo, and P. Dario, “Robotic delivery service in combined outdoor–indoor environments: technical analysis and user evaluation,” _Robotics and autonomous systems_ , vol. 103, pp. 56–67, 2018.
  * [2] J. J. Roldán, J. del Cerro, D. Garzón-Ramos, P. Garcia-Aunon, M. Garzón, J. De León, and A. Barrientos, “Robots in agriculture: State of art and practical experiences,” _Service robots_ , pp. 67–90, 2018\.
  * [3] M. Z. Zaheer, A. Mahmood, M. H. Khan, M. Astrid, and S.-I. Lee, “An anomaly detection system via moving surveillance robots with human collaboration,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2021, pp. 2595–2601.
  * [4] Y. Gu, J. Strader, N. Ohi, S. Harper, K. Lassak, C. Yang, L. Kogan, B. Hu, M. Gramlich, R. Kavi, _et al._ , “Robot foraging: Autonomous sample return in a large outdoor environment,” _IEEE Robotics & Automation Magazine_, vol. 25, no. 3, pp. 93–101, 2018.
  * [5] B. Choi, W. Lee, G. Park, Y. Lee, J. Min, and S. Hong, “Development and control of a military rescue robot for casualty extraction task,” _Journal of Field Robotics_ , vol. 36, no. 4, pp. 656–676, 2019.
  * [6] M. Aladem, S. Baek, and S. A. Rawashdeh, “Evaluation of image enhancement techniques for vision-based navigation under low illumination,” _Journal of Robotics_ , vol. 2019, 2019.
  * [7] T. Ji, A. N. Sivakumar, G. Chowdhary, and K. Driggs-Campbell, “Proactive anomaly detection for robot navigation with multi-sensor fusion,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 4975–4982, 2022.
  * [8] A. J. Sathyamoorthy, K. Weerakoon, T. Guan, J. Liang, and D. Manocha, “Terrapn: Unstructured terrain navigation through online self-supervised learning,” _arXiv preprint arXiv:2202.12873_ , 2022.
  * [9] S. S. Mansouri, C. Kanellakis, D. Kominiak, and G. Nikolakopoulos, “Deploying mavs for autonomous navigation in dark underground mine environments,” _Robotics and Autonomous Systems_ , vol. 126, p. 103472, 2020.
  * [10] K. Weerakoon, A. J. Sathyamoorthy, U. Patel, and D. Manocha, “Terp: Reliable planning in uneven outdoor environments using deep reinforcement learning,” in _2022 International Conference on Robotics and Automation (ICRA)_ , 2022, pp. 9447–9453.
  * [11] K. Weerakoon, S. Chakraborty, N. Karapetyan, A. J. Sathyamoorthy, A. S. Bedi, and D. Manocha, “Htron: Efficient outdoor navigation with sparse rewards via heavy tailed adaptive reinforce algorithm,” _arXiv preprint arXiv:2207.03694_ , 2022.
  * [12] P. Cai, S. Wang, Y. Sun, and M. Liu, “Probabilistic end-to-end vehicle navigation in complex dynamic environments with multimodal sensor fusion,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 3, pp. 4218–4224, 2020\.
  * [13] C. Debeunne and D. Vivet, “A review of visual-lidar fusion based simultaneous localization and mapping,” _Sensors_ , vol. 20, no. 7, p. 2068, 2020.
  * [14] Z. Wu, S. Pan, F. Chen, G. Long, C. Zhang, and S. Y. Philip, “A comprehensive survey on graph neural networks,” _IEEE transactions on neural networks and learning systems_ , vol. 32, no. 1, pp. 4–24, 2020.
  * [15] S. Casas, C. Gulino, R. Liao, and R. Urtasun, “Spagnn: Spatially-aware graph neural networks for relational behavior forecasting from sensor data,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_ , 2020, pp. 9491–9497.
  * [16] Z. Ravichandran, L. Peng, N. Hughes, J. D. Griffith, and L. Carlone, “Hierarchical representations and explicit memory: Learning effective navigation policies on 3d scene graphs using graph neural networks,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 9272–9279.
  * [17] T. Ji, S. T. Vuppala, G. Chowdhary, and K. Driggs-Campbell, “Multi-modal anomaly detection for unstructured and uncertain environments,” in _Conference on Robot Learning_. PMLR, 2021, pp. 1443–1455.
  * [18] D. Park, H. Kim, and C. C. Kemp, “Multimodal anomaly detection for assistive robots,” _Autonomous Robots_ , vol. 43, no. 3, pp. 611–629, 2019.
  * [19] S. Zhou, J. Xi, M. W. McDaniel, T. Nishihata, P. Salesses, and K. Iagnemma, “Self-supervised learning to visually detect terrain surfaces for autonomous robots operating in forested terrain,” _Journal of Field Robotics_ , vol. 29, no. 2, pp. 277–297, 2012.
  * [20] L. Wellhausen, R. Ranftl, and M. Hutter, “Safe robot navigation via multi-modal anomaly detection,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 1326–1333, 2020.
  * [21] G. Kahn, P. Abbeel, and S. Levine, “Badgr: An autonomous self-supervised learning-based navigation system,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 1312–1319, 2021.
  * [22] K. Jin, H. Wang, C. Liu, Y. Zhai, and L. Tang, “Graph neural network based relation learning for abnormal perception information detection in self-driving scenarios,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 8943–8949.
  * [23] M. Sorokin, J. Tan, C. K. Liu, and S. Ha, “Learning to navigate sidewalks in outdoor environments,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 3906–3913, 2022.
  * [24] K. Katyal, K. Popek, C. Paxton, P. Burlina, and G. D. Hager, “Uncertainty-aware occupancy map prediction using generative networks for robot navigation,” in _2019 International Conference on Robotics and Automation (ICRA)_ , 2019, pp. 5453–5459.
  * [25] L. Rusli, B. Nurhalim, and R. Rusyadi, “Vision-based vanishing point detection of autonomous navigation of mobile robot for outdoor applications,” _Journal of Mechatronics, Electrical Power, and Vehicular Technology_ , vol. 12, no. 2, pp. 117–125, 2021.
  * [26] F. Liu, X. Li, S. Yuan, and W. Lan, “Slip-aware motion estimation for off-road mobile robots via multi-innovation unscented kalman filter,” _IEEE Access_ , vol. 8, pp. 43 482–43 496, 2020.
  * [27] Y. Yan, B. Zhang, J. Zhou, Y. Zhang, X. Liu, _et al._ , “Real-time localization and mapping utilizing multi-sensor fusion and visual–imu–wheel odometry for agricultural robots in unstructured, dynamic and gps-denied greenhouse environments,” _Agronomy_ , vol. 12, no. 8, p. 1740, 2022.
  * [28] Y. Qu, M. Yang, J. Zhang, W. Xie, B. Qiang, and J. Chen, “An outline of multi-sensor fusion methods for mobile agents indoor navigation,” _Sensors_ , vol. 21, no. 5, p. 1605, 2021.
  * [29] G. Welch, G. Bishop, _et al._ , “An introduction to the kalman filter,” 1995\.
  * [30] A. E. B. Velasquez, V. A. H. Higuti, M. V. Gasparino, A. N. Sivakumar, M. Becker, and G. Chowdhary, “Multi-sensor fusion based robust row following for compact agricultural robots,” _arXiv preprint arXiv:2106.15029_ , 2021\.
  * [31] Y. Su, T. Wang, S. Shao, C. Yao, and Z. Wang, “Gr-loam: Lidar-based sensor fusion slam for ground robots on complex terrain,” _Robotics and Autonomous Systems_ , vol. 140, p. 103759, 2021.
  * [32] J. Liang, K. Weerakoon, T. Guan, N. Karapetyan, and D. Manocha, “Adaptiveon: Adaptive outdoor navigation method for stable and reliable motions,” _arXiv preprint arXiv:2205.03517_ , 2022.
  * [33] Y.-S. Shin, Y. S. Park, and A. Kim, “Direct visual slam using sparse depth for camera-lidar system,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 5144–5151.
  * [34] D. Park, H. Kim, and C. C. Kemp, “Multimodal anomaly detection for assistive robots,” _Autonomous Robots_ , vol. 43, no. 3, pp. 611–629, 2019.
  * [35] T. Ji, S. T. Vuppala, G. Chowdhary, and K. Driggs-Campbell, “Multi-modal anomaly detection for unstructured and uncertain environments,” in _Conference on Robot Learning_. PMLR, 2021, pp. 1443–1455.
  * [36] G. Kahn, P. Abbeel, and S. Levine, “Land: Learning to navigate from disengagements,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 2, pp. 1872–1879, 2021.
  * [37] M. Welling and T. N. Kipf, “Semi-supervised classification with graph convolutional networks,” in _J. International Conference on Learning Representations (ICLR 2017)_ , 2016.
  * [38] K. K. Thekumparampil, C. Wang, S. Oh, and L.-J. Li, “Attention-based graph neural network for semi-supervised learning,” _arXiv preprint arXiv:1803.03735_ , 2018.
  * [39] D. Busbridge, D. Sherburn, P. Cavallo, and N. Y. Hammerla, “Relational graph attention networks,” _arXiv preprint arXiv:1904.05811_ , 2019.
  * [40] T. Guan, D. Kothandaraman, R. Chandra, A. J. Sathyamoorthy, K. Weerakoon, and D. Manocha, “Ga-nav: Efficient terrain segmentation for robot navigation in unstructured outdoor environments,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 3, pp. 8138–8145, 2022.
  * [41] H. Zhao, J. Shi, X. Qi, X. Wang, and J. Jia, “Pyramid scene parsing network,” in _2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2017, pp. 6230–6239.
  * [42] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [43] S. Woo, J. Park, J.-Y. Lee, and I. S. Kweon, “Cbam: Convolutional block attention module,” in _Proceedings of the European conference on computer vision (ECCV)_ , 2018, pp. 3–19.
  * [44] C. R. Qi, H. Su, K. Mo, and L. J. Guibas, “Pointnet: Deep learning on point sets for 3d classification and segmentation,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 652–660.
  * [45] E. Rosten, R. Porter, and T. Drummond, “Faster and better: A machine learning approach to corner detection,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 32, no. 1, pp. 105–119, 2008.
  * [46] J. Zhang and S. Singh, “Loam: Lidar odometry and mapping in real-time.” in _Robotics: Science and Systems_ , vol. 2, no. 9. Berkeley, CA, 2014, pp. 1–9.
  * [47] A. Mojsilovic, J. Gomes, and B. E. Rogowitz, “Isee: Perceptual features for image library navigation,” in _Human Vision and Electronic Imaging VII_ , vol. 4662. SPIE, 2002, pp. 266–277.
  * [48] E. Rosten and T. Drummond, “Machine learning for high-speed corner detection,” in _European conference on computer vision_. Springer, 2006, pp. 430–443.
  * [49] D. Fox, W. Burgard, and S. Thrun, “The dynamic window approach to collision avoidance,” _IEEE Robotics & Automation Magazine_, vol. 4, no. 1, pp. 23–33, 1997.


# MAPPING PROPERTIES OF CERTAIN INTEGRAL OPERATORS

Victor Polunin  Chair of Applied Mathematics and Computer Modeling
Belgorod State National Research University
Pobedy street 85, Belgorod 308015, Russia<EMAIL_ADDRESS>,  Vladimir
Vasilyev  Chair of Applied Mathematics and Computer Modeling
Belgorod State National Research University
Pobedy street 85, Belgorod 308015, Russia<EMAIL_ADDRESS>and  Nelly Erygina
Chair of Mathematics and Education in Natural Sciences
Belgorod State National Research University
Studencheskaya street 14, Belgorod 308007, Russia<EMAIL_ADDRESS>
###### Abstract.

We study mapping properties of two-dimensional linear integral operators in
some weighted spaces with special kernels. The considered spaces are certain
variant of Sobolev–Slobodetskii spaces and their generalizations related to
Banach spaces. Sufficient conditions for boundedness for such operators in
these spaces are obtained.

###### Key words and phrases:

integral operator, weighted norm inequality

###### 2010 Mathematics Subject Classification:

Primary: 45P05; Secondary: 47G10

##  1\. Introduction

Some systems of integral equations arise under studying model pseudo-
differential equations in plane corners [1, 2], for example when we consider
simple boundary conditions like Dirichlet or Neumann conditions. According to
this fact it seems interesting to study some mapping properties of such
integral operators. These two-dimensional integral operator act in Fourier
images of Sobolev–Slobodetski spaces which are weighted spaces with prescribed
weight. These mapping properties play important role under studying discrete
boundary value problems [4, 5] since we try to compare such integral operators
with their truncated analogues, these arise under considering discrete
boundary value problems. Let us note that integral operators and corresponding
equations do not arise under studying model boundary value problems in a half-
space, there are systems of linear algebraic equations as a result of
reduction [3].

We will give here a simple example. Let $\Omega\subset\mathbb{R}^{2}$ be the
first quadrant. We consider the following Dirichlet boundary value problem [1,
5] in Sobolev–Slobodetskii space $H^{s}(\Omega)$

(1) $\left\\{\begin{array}[]{rcl}(Au)(x)=0,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ x\in\Omega,\\\
{u}_{|_{x_{1}=0}}=f(x_{2}),\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
{u}_{|_{x_{2}=0}}=g(x_{1})\end{array}\right.$

where $A$ is a pseudo-differential operator with the symbol $A(\xi)$
satisfying the condition

$c_{1}(1+|\xi|)^{\alpha}\leq|A(\xi)|\leq c_{2}(1+|\xi|)^{\alpha}$

and admitting the wave factorization [1] with respect to $\Omega$

$A(\xi)=A_{\neq}(\xi)A_{=}(\xi)$

with the index $\ae$ such that $\ae-s=1+\delta,|\delta|<1/2$ then the
Dirichlet problem (1) can be reduced to the system of linear integral
equations

$\left\\{\begin{array}[]{rcl}&\int\limits_{-\infty}^{\infty}K_{1}(\xi)C(\xi_{1})d\xi_{1}+D(\xi_{2})=F(\xi_{2})\\\
&C(\xi_{1})+\int\limits_{-\infty}^{\infty}K_{2}(\xi)D(\xi_{2})d\xi_{2}=G(\xi_{1})\end{array}\right.$

with respect to two unknown functions $C,D$. These functions should belong to
Fourier image of certain $H^{s}$-space with certain $s<0$.

##  2\. Integral operators

We study the following integral operator

$(Kf)(x)=\int\limits_{-\infty}^{+\infty}K(x,y)f(y)dy$

in functional spaces with the following norm

$||f||_{s}=\left(\int\limits_{\mathbb{R}}|f(x)|^{2}(1+|x|)^{2s}dx\right)^{1/2}.$

We assume that the kernel $K(x,y)$ satisfies the following condition

(2) $|K(x,y)|\sim(1+|x|+|y|)^{-\ae},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \ae\in\mathbb{R}.$

Theorem 1. The operator $K$ is a linear bounded operator,

$K:{H}^{s_{1}}(\mathbb{R})\rightarrow{H}^{s_{2}}(\mathbb{R}),$

where $s_{1}<0,\ae>\max\\{1/2-s_{1},1+s_{2}-s_{1}\\}$.

###### Proof.

Let us verify.

$||Kf||^{2}_{s_{2}}=\int\limits_{-\infty}^{+\infty}(1+|x|)^{2s_{2}}|(Kf)(x)|^{2}dx=$
$=\int\limits_{-\infty}^{+\infty}(1+|x|)^{2s_{2}}\left|\int\limits_{-\infty}^{+\infty}K(x,y)f(y)dy\right|^{2}dx\leq$
$\leq\int\limits_{-\infty}^{+\infty}(1+|x|)^{2s_{2}}\left(\int\limits_{-\infty}^{+\infty}|K(x,y)||f(y)|dy\right)^{2}dx\leq$
$\leq\leavevmode\nobreak\
const\int\limits_{-\infty}^{+\infty}(1+|x|)^{2s_{2}}\left(\int\limits_{-\infty}^{+\infty}(1+|x|+|y|)^{-\ae}||f(y)|dy\right)^{2}dx.$

In the inner integral, we apply the Cauchy–Schwartz inequality introducing the
factors $(1+|y|)^{-s_{1}}$ and $(1+|y|)^{s_{1}}$ form the first and the second
term and taking into account that $(1+|y|)^{-s_{1}}\leq(1+|x|+|y|)^{-s_{1}}$.
Then we obtain

$||Kf||^{2}_{s_{2}}\leq\leavevmode\nobreak\
cons||f||^{2}_{s_{1}}\int\limits_{-\infty}^{+\infty}(1+|x|)^{2s_{2}}\left(\int\limits_{-\infty}^{+\infty}(1+|x|+|y|)^{-2(s_{1}+\ae)}dy\right)dx\leq$
$\leq\leavevmode\nobreak\
const||f||^{2}_{s_{1}}\int\limits_{0}^{+\infty}(1+|x|)^{2s_{2}}(1+|x|)^{-2(s_{1}+\ae-1/2)}dx$

since $-2(s_{1}+\ae)<-1$. Thus, we conclude

$||Kf||^{2}_{s_{2}}\leq\leavevmode\nobreak\
cons||f||^{2}_{s_{1}}\int\limits_{0}^{+\infty}(1+|x|)^{2(s_{2}-s_{1}-\ae+1/2)}dx\leq\leavevmode\nobreak\
cons||f||^{2}_{s_{1}}$

since $2(s_{2}-s_{1}-\ae+1/2)<-1$. ∎

Corollary 1. The operator $K$ is a linear bounded operator,

$K:{H}^{s}(\mathbb{R})\rightarrow{H}^{s}(\mathbb{R}),$

where $s<0,\ae>\max\\{1/2-s,1\\}$.

##  3\. Some generalizations: $H^{s,p}$-spaces

###  3.1. First variant

This section is devoted to a generalization of the above result to Banach
spaces $H^{s,p}(\mathbb{R}),1<p<\infty$. Let us write that for $f\in
H^{s,p}(\mathbb{R})$

$||f||_{s,p}=\left(\int\limits_{\mathbb{R}}|f(x)|^{p}(1+|x|)^{ps}dx\right)^{1/p}.$

Obviously, for $p=2$ we have $H^{s,2}(\mathbb{R})\equiv H^{s}(\mathbb{R})$.
Let us remind also that the pair $(p,q)$ is so that $1/p+1/q=1$.

Theorem 2. Let the kernel $K(x,y)$ satisfies the condition (2). Then the
operator $K$ is a linear bounded operator

$K:{H}^{s_{1},p_{1}}(\mathbb{R})\rightarrow{H}^{s_{2},p_{2}}(\mathbb{R}),$

where $s_{1},s_{2},\ae,p_{1},p_{2}$ such that
$s_{1}<0,\ae>\max\\{1/{q_{1}}-s_{1},1/{p_{2}}+1/{q_{1}}+s_{2}-s_{1}\\}$.

###### Proof.

We estimate

$||Kf||^{p_{2}}_{s_{2},p_{2}}=\int\limits_{-\infty}^{+\infty}(1+|x|)^{p_{2}s_{2}}|(Kf)(x)|^{p_{2}}dx=$
$=\int\limits_{-\infty}^{+\infty}(1+|x|)^{p_{2}s_{2}}\left|\int\limits_{-\infty}^{+\infty}K(x,y)f(y)dy\right|^{p_{2}}dx\leq$
$\leq\int\limits_{-\infty}^{+\infty}(1+|x|)^{p_{2}s_{2}}\left(\int\limits_{-\infty}^{+\infty}|K(x,y)||f(y)|dy\right)^{p_{2}}dx\leq$
$\leq\leavevmode\nobreak\
const\int\limits_{0}^{+\infty}(1+|x|)^{p_{2}s_{2}}\left(\int\limits_{-\infty}^{+\infty}(1+|x|+|y|)^{-\ae}|f(y)|dy\right)^{p_{2}}dx.$

Let us consider the inner integral

$\int\limits_{-\infty}^{+\infty}(1+|x|+|y|)^{-\ae}|f(y)|dy=$
$=\int\limits_{-\infty}^{+\infty}(1+|y|)^{-s_{1}}(1+|x|+|y|)^{-\ae}|f(y)|(1+|y|)^{s_{1}}dy$

Applying the Hölder inequality and taking into account the inequality
$(1+|y|)^{-s_{1}}\leq(1+|x|+|y|)^{-s_{1}}$ we have

$\int\limits_{-\infty}^{+\infty}(1+|x|+|y|)^{-\ae}||f(y)|dy\leq||f||_{s_{1},p_{1}}\left(\int\limits_{0}^{+\infty}(1+|x|+|y|)^{-q_{1}(s_{1}+\ae)}dy\right)^{\frac{1}{q_{1}}}\leq$
$\leq\leavevmode\nobreak\
const||f||_{s_{1},p_{1}}(1+|x|)^{\frac{-q_{1}(s_{1}+\ae)+1}{q_{1}}},$

since $-q_{1}(s_{1}+\ae)<-1\Longleftrightarrow\ae>1/{q_{1}}-s_{1}$.

If so then we have

$||Kf||^{p_{2}}_{s_{2},p_{2}}\leq\leavevmode\nobreak\
cons||f||^{p_{2}}_{s_{1},p_{1}}\int\limits_{0}^{+\infty}(1+|x|)^{p_{2}(1/{q_{1}}-s_{1}-\ae+s_{2})}dx\leq\leavevmode\nobreak\
const||f||^{p_{2}}_{s_{1},p_{1},}$

since
$p_{2}(1/{q_{1}}-s_{1}-\ae+s_{2})<-1\Longleftrightarrow\ae>1/{p_{2}}+1/{q_{1}}+s_{2}-s_{1}$.
∎

Corollary 2. Let the kernel $K(x,y)$ satisfies the condition (2). Then the
operator $K$ is a linear bounded operator

$K:{H}^{s,p}(\mathbb{R})\rightarrow{H}^{s,p}(\mathbb{R}),$

where $s,\ae,p$ such that $s<0,\ae>\max\\{1/{q}-s,1\\}$.

###  3.2. Second variant

This section is devoted to another generalization of the above result to
Banach spaces $H^{s,p}(\mathbb{R}),1<p<\infty$. We say that for $f\in
H^{p,s}(\mathbb{R})$

$||f||_{p,s}=\left(\int\limits_{\mathbb{R}}|f(x)|^{p}(1+|x|)^{2s}dx\right)^{1/p}.$

Obviously, for $p=2$ we have $H^{s,2}(\mathbb{R})\equiv H^{s}(\mathbb{R})$.
Let us remind also that the pair $(p,q)$ is so that $1/p+1/q=1$.

Theorem 3. Let the kernel $K(x,y)$ satisfies the condition (2). Then the
operator $K$ is a linear bounded operator

$K:{H}^{p_{1},s_{1}}(\mathbb{R})\rightarrow{H}^{p_{2},s_{2}}(\mathbb{R}),$

where $s_{1},s_{2},\ae,p_{1},p_{2}$ such that
$s_{1}<0,\ae>\max\\{1/{q_{1}}-2s_{1}/{p_{2}},1/{p_{2}}+1/{q_{1}}+2s_{2}/{p_{2}}-2s_{1}/{p_{1}}\\}$.

###### Proof.

Let us estimate the norm

$||Kf||^{p_{2}}_{p_{2},s_{2}}=\int\limits_{-\infty}^{+\infty}(1+|x|)^{2s_{2}}|(Kf)(x)|^{p_{2}}dx=$
$=\int\limits_{-\infty}^{+\infty}(1+|x|)^{2s_{2}}\left|\int\limits_{-\infty}^{+\infty}K(x,y)f(y)dy\right|^{p_{2}}dx\leq$
$\leq\int\limits_{-\infty}^{+\infty}(1+|x|)^{2s_{2}}\left(\int\limits_{-\infty}^{+\infty}|K(x,y)||f(y)|dy\right)^{p_{2}}dx\leq$
$\leq\leavevmode\nobreak\
const\int\limits_{0}^{+\infty}(1+|x|)^{2s_{2}}\left(\int\limits_{-\infty}^{+\infty}(1+|x|+|y|)^{-\ae}|f(y)|dy\right)^{p_{2}}dx.$

Let us consider the inner integral

$\int\limits_{-\infty}^{+\infty}(1+|x|+|y|)^{-\ae}|f(y)|dy=$
$=\int\limits_{-\infty}^{+\infty}(1+|y|)^{-2s_{1}/{p_{1}}}(1+|x|+|y|)^{-\ae}|f(y)|(1+|y|)^{2s_{1}/{p_{1}}}dy$

Applying the Hölder inequality and taking into account the inequality
$(1+|y|)^{-2s_{1}/{p-1}}\leq(1+|x|+|y|)^{-2s_{1}/{p-1}}$ we have

$\int\limits_{-\infty}^{+\infty}(1+|x|+|y|)^{-\ae}||f(y)|dy\leq||f||_{p_{1},s_{1}}\left(\int\limits_{0}^{+\infty}(1+|x|+|y|)^{-q_{1}(2s_{1}/{p_{1}}+\ae)}dy\right)^{\frac{1}{q_{1}}}\leq$
$\leq\leavevmode\nobreak\
const||f||_{p_{1},s_{1}}(1+|x|)^{\frac{-q_{1}(2s_{1}/{p_{1}}+\ae)+1}{q_{1}}},$

since
$-q_{1}(2s_{1}/{p_{1}}+\ae)<-1\Longleftrightarrow\ae>1/{q_{1}}-2s_{1}/{p_{2}}$.

If so then we have

$||Kf||^{p_{2}}_{p_{2},s_{2}}\leq\leavevmode\nobreak\
cons||f||^{p_{2}}_{p_{1},s_{1}}\int\limits_{0}^{+\infty}(1+|x|)^{p_{2}(1/{q_{1}}-2s_{1}/{p_{1}}-\ae)+2s_{2}}dx\leq\leavevmode\nobreak\
const||f||^{p_{2}}_{p_{1},s_{1},}$

since
$p_{2}(1/{q_{1}}-2s_{1}/{p_{1}}-\ae)+2s_{2}<-1\Longleftrightarrow\ae>1/{p_{2}}+1/{q_{1}}+2s_{2}/{p_{2}}-2s_{1}/{p_{1}}$.
∎

Corollary 3. Let the kernel $K(x,y)$ satisfies the condition (2). Then the
operator $K$ is a linear bounded operator

$K:{H}^{p,s}(\mathbb{R})\rightarrow{H}^{p,s}(\mathbb{R}),$

where $s,\ae,p$ such that $s<0,\ae>\max\\{1/{q}-2s/{p},1\\}$.

## Conclusion

The considered spaces $H^{s,p}(\mathbb{R})$ and $H^{p,s}(\mathbb{R})$ like
almost the analogous spaces considered in [6, 7]. But we work in Fourier
images, and the introduced spaces are usual weighted spaces with power
weights.

## References

  * [1] V. B. Vasil’ev, Wave Factorization of Elliptic Symbols: Theory and Applications. Introduction to the Theory of Boundary Value Problems in Non-Smooth Domains. Dordrecht–Boston–London: Kluwer Academic Publishers, 2000.
  * [2] A. V. Vasilyev, V. B. Vasilyev, Difference equations and boundary value problems. In: Pinelas, S., Dŏslá, Z., Dŏslý, O., Kloeden, P. (eds.) Differential and Difference Equations and Applications, Springer Proc. Math. & Stat. 2016, 164 421–432.
  * [3] G. Eskin, Boundary Value Problems for Elliptic Pseudodifferential Equations. Providence, RI: AMS, 1981.
  * [4] A. A. Mashinets, A. V. Vasilyev, V. B. Vasilyev, On discrete Neumann problem in a quadrant. Lobachevskii J. Math. 2023, 44 (3) 1011–1021.
  * [5] A. V. Vasilyev, V. B. Vasilyev, A.‘A. Khodyrevz, The discrete Dirichlet problem. Solvability and approximation properties. J. Math. Sci. 2023, 270 (5) 654–664.
  * [6] H. Triebel, Theory of Function Spaces. Basel–Boston–Stuttgart: Birkhäuser Verlag, 1983.
  * [7] M.W. Wong, An Introduction to Pseudo-Differential Operators. Singapore–New Jersey–London–Hong Kong: World Scientific, 1991.


#  How big does a population need to be before demographers can ignore
individual-level randomness in demographic events?

John Bryant  Bayesian Demography Limited, Christchurch, New Zealand  Tahu
Kukutai  Te Ngira Institute for Population Research, University of Waikato,
Hamilton, New Zealand  Junni L. Zhang  National School of Development, Peking
University, Beijing, China

###### Abstract

BACKGROUND When studying a national-level population, demographers can safely
ignore the effect of individual-level randomness on age-sex structure. When
studying a single community, or group of communities, however, the potential
importance of individual-level randomness is less clear.

OBJECTIVE We seek to measure the effect of individual-level randomness in
births and deaths on standard summary indicators of age-sex structure, for
populations of different sizes. We focus on demographic conditions typical of
historical populations.

METHODS We conduct a microsimulation experiment where we simulate events and
age-sex structure under a range of settings for demographic rates and
population size.

RESULTS Individual-level randomness strongly affects age-sex structure for
populations of about 100, but has a much smaller effect on populations of
1,000, and a negligible effect on populations of 10,000.

CONCLUSIONS Analyses of age-sex structure in historical populations with sizes
on the order 100 must account for individual-level randomness in demographic
events. Analyses of populations with sizes on the order of 1,000 may need to
make some allowance for individual-level variation, but other issues, such as
measurement error, probably deserve more attention. Analyses of populations of
10,000 can safely ignore individual-level variation.

CONTRIBUTION We provide rules of thumb for when it is safe to ignore
individual-level randomness in demographic events when analysing age-sex
structure in small historical populations.

##  1 Introduction

Demography has traditionally emphasised national populations, which are large
enough that randomness in individual demographic events such as births and
deaths averages out and can be ignored. In such populations, variation in
population size and structure is almost entirely driven by variation in the
underlying birth, death, and migration rates (Lee,, 1998; Preston et al.,,
2001, p. 157).

The study of demographic change can, however, often benefit from narrowing the
focus down to individual communities or localities, and reconstructing events
and processes at the micro level (Caldwell et al.,, 1987; Caldwell,, 1990).
Local-level demographic studies have generated fundamental knowledge about
historical patterns of survivorship and fertility (Henry,, 1968). Historical
demography in France began with multiple studies of villages, using parish
records to reconstitute families and their demographic behaviours (Rosenthal,,
1997). These studies, and later ones in England (Wrigley and Schofield,, 1989)
and continental Europe (Coale and Watkins,, 1986), profoundly changed
understanding about population change and its social and economic drivers.
Local-level demographic studies can contribute to the understanding of social
change elsewhere in the world, such as documenting the demographic effects of
colonisation on Indigenous peoples (Cook and Borah,, 1971; Covey et al.,,
2011; Petersen et al.,, 1975), where large-scale systematic demographic data
are typically lacking (Dobyns,, 1966; Roberts,, 1989; Thornton,, 1980).

If a population is small, ignoring individual-level randomness may be unwise.
In a population of 100 people, for instance, the number of deaths can easily
halve or double from one year to the next with no change in the underlying
risk of dying. (For a demonstration, see Appendix A.) If deaths, births, or
migrations vary randomly from year to year, we would expect population size
and structure to vary randomly as well. A run of higher-than-average birth
counts, for instance, should increase population size and produce a younger
age structure.

The relationship between individual-level randomness in demographic events and
population-level randomness in size and structure is, however, complicated,
because the relations between demographic events and population size and
structure are subject to lags and feedback loops. The number of births,
deaths, and migrations that occur this year, for instance, reflects the size
and structure of the population this year, but the size and structure of the
population this year is itself a product of births, deaths, and migrations
during previous years.

The complexity of the relationship between births, deaths, migration, and
population makes it difficult to know how much randomness to expect in small
populations. Having reliable intuitions about random variation in small
populations is, however, crucial to making inferences about the demography of
these populations. Consider, for instance, a demographer who has data on the
population share of old people in two communities of around 1,000 people. If
5% of people in one community are aged 60+ and 10% in the other are aged 60+,
how confident should the demographer be that the two communities have
different demographic regimes?

One way to build up intuitions about the amount of random variation to expect
in small populations would be to assemble large quantities of accurate data on
birth rates, death rates, and population size and structure from a wide
variety of settings, and see how much variability there is across communities.
However, for many populations,, particularly historical populations, such data
are simply not available, so the direct observational approach is not
feasible.

Another way to build intuitions is to use mathematical theory. Caswell, (2001,
ch. 15) presents a sophisticated mathematical framework for studying the
effects of individual-level variability on population aggregates.

We, however, take a more direct approach that requires only elementary
mathematics, and hence is hopefully more accessible to applied demographers.
Our approach is computer simulation. We randomly generate artificial life
histories which we aggregate up into populations classified by age and sex.
The strategy of using computer power as an alternative to mathematical theory
is common in contemporary applied statistics (Gelman and Vehtari,, 2021),
though there is also a long tradition in demography of using micro-simulation
to study demographic processes (Wachter et al.,, 1978; Van Imhoff and Post,,
1998; Zagheni,, 2015).

Our simulations use three different initial population sizes: 100, 1,000, and
10,000. Having different population sizes allows us to examine the
relationship between population size and random variability. We also simulate
a variety of different mortality and fertility regimes. This gives us a
yardstick for assessing the importance of randomness: we can compare the
variation attributable to randomness in births and deaths with the variation
attributable to differences in underlying demographic rates. For simplicity,
we assume zero migration, and we assume that the underlying age-specific birth
rates and death rates are constant throughout the simulation period.

The study is designed to mimic the demographic conditions that are studied by
historical demographers, with high mortality and high fertility. We also
assume that, due to limitations in available data, demographers must work with
broad indicators of population structure, such as dependency ratios.

Our simulations suggest that, in populations of around 100 people, variation
due to the randomness of demographic events is large relative to the variation
due to differences in underlying rates. When population size reaches 1,000,
however, variation due to the randomness of demographic events is already
relatively modest. By the time population size reaches 10,000, variation due
to the randomness of demographic events is negligible.

All data and code for the simulations are available at [names of github
repositories].

##  2 Design of the simulations

###  2.1 Inputs

The simulations require values for population size, mortality, the sex ratio
at birth, and fertility. The inputs are summarized in Table 1.

Table 1: Inputs to the simulation Input |  Description
---|---
$N$ |  Initial population size, set by us
$m_{as}$, $L_{as}$ |  Mortality rates and expected person-years lived, from West model life tables
$p_{s}$ |  Proportion of births of sex $s$, based on a sex ratio at birth of 105
$g_{a}$ |  Age-specific fertility rates, scaled to sum to 1, from the ‘Booth standard’
$F$ |  The total fertility rate, supplied by us.

Population size $N$ is the total number of people in the population at the
start of the simulation. Values for $N$ are chosen by us to be 100, 1,000, and
10,000.

The $m_{as}$ in the second row of Table 1 is the mortality rate for people of
sex $s$ in age group $a$. All age groups have widths of 5 years, apart from
the final age group, which includes everyone aged 95 years and over. The
quantity $L_{as}$ is the number of years that a newborn baby of sex $s$ could
expect to live in age group $a$, under prevailing mortality rates.

We obtain values for $m_{as}$ and $L_{as}$ from West model life tables (Coale
et al.,, 1983), as reported in the R package demogR (Jones,, 2007). The West
model life tables are a set of mortality rates, and derived quantities, that
demographers use to represent typical mortality patterns across a wide range
of settings. West model life tables are indexed by level, with lower lowers
indicating higher mortality rates, and hence lower life expectancies. Table 2
shows the West model life table levels that we use in the simulations, along
with the associated life expectancies. These life tables represent mortality
conditions in pre-industrial societies.

Table 2: Life expectancy for each level of the West model life tables used in the simulations | Level 1 | Level 5 | Level 9
---|---|---|---
Female life expectancy | 20.0 | 30.0 | 40.0
Male life expectancy | 18.0 | 27.7 | 37.3

We set $p_{\text{F}}$, the proportion of births that are female, to 0.488.
This is equivalent to assuming that there are 105 male births per 100 female
births, which is a standard assumption in demography.

The ‘Booth standard’ age pattern for fertility describes a typical
distribution of fertility across women’s reproductive years. It is the
starting point for some widely-used methods for fertility estimation (Booth,,
1984; Moultrie et al.,, 2013). Values for the Booth standard are shown in
Table 3.

Table 3: Proportional distribution of births by age of mother, under the Booth standard 10–14 | 15–19 | 20–24 | 25–29 | 30–34 | 35–39 | 40–44 | 45–49
---|---|---|---|---|---|---|---
0.003 | 0.133 | 0.241 | 0.231 | 0.188 | 0.134 | 0.062 | 0.008

The total fertility rate (TFR) is the average number of births a woman could
expect to have over her life if she survived to the end of the reproductive
ages, under prevailing fertility rates. Values for the TFR are chosen by us,
and are 4, 5, and 6.

To connect deaths and births to population, classifying deaths and births by
age, sex, and period is not sufficient: we also need to classify them by
‘Lexis triangle’. Figure 1 depicts Lexis triangles for deaths for age groups
0–4 and 5–9. The triangles are defined by the solid lines marking out age
groups and periods, and the dashed lines on the diagonal. A death for age
group $a$ belongs to the upper Lexis triangle if the person dying was in age
group $a$ at the start of the period; otherwise the death belongs to the lower
Lexis triangle.

Figure 1: Disaggregating events by Lexis triangle. Event A belongs to the
lower Lexis triangle for age group 0–4, while event B belongs to the upper
Lexis triangle for age group 0–4. The grey area is in children aged less than
1 year.

With fertility, we assume that, for any given age, the rate in the lower Lexis
triangle is identical to the rate in the upper Lexis triangle. For mortality
we make the equivalent assumption for all age groups except the first. For
ages 0–4, assuming equality is too crude an approximation. The mortality rate
experienced by infants is much higher than the mortality rate experienced by
children aged 1–4, and values in the lower and upper Lexis triangles need to
be weighted accordingly. Appendix B describes our method for doing so.

###  2.2 Initial populations

Choosing an appropriate initial age-sex structure for the simulations is not
trivial. Unlike populations in which births and deaths are deterministic,
populations in which births and deaths are random do not have straightforward
stable population structures that might serve as a neutral starting point
(Caswell,, 2001, ch. 15). We instead use the stationary populations implied by
the associated life tables in a deterministic model. We start with the age
structure that would arise in a deterministic setting where births equalled
deaths, in which the share of the population in age group $a$ and sex $s$ is
proportional to $L_{as}$ (Keyfitz and Caswell,, 2005, ch. 5). We then scale
this population so that it equals $N$ and then round each cell to the nearest
integer.

###  2.3 Simulating death and survival

As a cohort ages, it moves diagonally across the Lexis diagram, alternating
between Lower and Upper triangles. We simulate death and survival one Lexis
triangle at a time. Within each Lexis triangle, we simulate individuals. To
simulate an individual, we go through the following steps:

  1. 1.

Draw a value for the amount of time left before the individual exits the
current Lexis triangle (which for a lower Lexis triangle means reaching time
$t+5$ and for an upper Lexis triangle means attaining age $a+5$.) We assume
that, at the start of each Lexis triangle, individuals are spaced evenly
across the time axis or the age axis, which allows us to derive the amount of
time left by drawing from a uniform distribution with minimum 0 and maximum 5.

  2. 2.

Draw a value for the amount of time before the individual dies, given the
mortality rates in the Lexis triangle. This is done by drawing from an
exponential distribution with the associated rate.

  3. 3.

If the second value is shorter than first, increment the number of deaths for
that Lexis triangle by one, and remove the individual from the simulated
population; otherwise include the individual in the population entering the
next Lexis triangle for that cohort.

###  2.4 Simulating births

Births are also simulated one triangle at a time. If the individual being
simulated is a female in the reproductive ages, then we multiply the amount of
time that the individual spent in the Lexis triangle by the fertility rate for
that triangle, to obtain an expected number of births, and then draw from a
Poisson distribution with that expected number. These births are added to the
population entering the lower Lexis triangle for age group 0–4.

###  2.5 Summary measures of population structure

We examine the effect of individual-level randomness on four summary measures
of population structure. We use these particular indicators on the grounds
that they are of substantive interest in a wide variety of demographic
studies, and can often be calculated from the type of data that is available
for historical populations, such as counts of children and adults. The
indicators are as follows:

  * •

Percent of population aged 0–14

  * •

Percent of population aged 60+

  * •

Child-woman ratio: The ratio of children aged 0–14 to women aged 15+

  * •

Dependency ratio: The ratio of people aged 0–14 and 60+ to people aged 15–59.

###  2.6 Implementation

The data preparation and the analysis of results is done in R, while the
simulation code is written in C++, for speed. We do 10,000 simulation runs for
each combination of inputs. The total time required for all simulation runs
combined is around 15 minutes. As noted above, all the code is available at
[addresses of github repositories].

###  2.7 Alternative approach to simulation

We also constructed simulations using the estimateAccount function in R
package demest111github.com/statisticsnz/demest. The results were very similar
to those from the microsimulation, but the calculations took 20 hours rather
than 15 minutes.

##  3 Results

###  3.1 Four simulated populations

(a) Birth and death counts for five-year periods.

(b) Population counts, by age and sex, at year 100

Figure 2: Results from four randomly-chosen simulations generated using a West
level 5 life table and a TFR of 5.

To give a sense of what the simulations contain, we first present output from
four randomly-chosen simulation runs, generated using a West level 5 life
table and a total fertility rate of 5. Figure 2(a) shows birth and death
rates. The populations in the top row of Figure 2(a) were generated using
initial population size $N=$ 100, and the populations in the bottom row were
generated using initial population size $N=$ 10,000. The two small populations
are much more variable over time than the two large populations, and look less
similar to each other. Figure 2(b) shows age-sex structure at year 100. The
small populations again differ much more than the large populations.

###  3.2 Distribution of summary indicators for total fertility rate 5 and
life table level 5

The results from the four simulation runs shown in Figure 2 are consistent
with the idea that smaller population size is associated with more
variability. With only four simulated populations, however, it is difficult to
assess the strength of the relationship. To do this, we need more simulation
runs.

Figure 3 shows distributions of summary measures over 10,000 runs. Values for
the summary measures evolve over time, so the figure shows values for each
year of the simulations. In all simulation runs in Figure 3, the total
fertility rate is 5 and the West life table level is 5. Each panel shows a
different combination of summary indicator and initial population size. The
top-left panel, for instance, shows the percent of the population aged 0–14 in
a population with starting size 100.

The white lines in Figure 3 show median values across the indicator across all
10,000 simulations. The lower limit of the dark grey bands represents the 25%
quantile: 25% of simulations, at each time point, have values for the
indicator that are lower than this limit. The upper limit of the dark grey
bands represent the 75% quantile: 75% of simulation runs have values that are
lower than this limit. The light grey bands represent the interval between the
2.5% and 97.5% quantiles. The wider the dark grey and light grey bands, the
more variability there is.

Figure 3: Distributions of summary indicators, over time, for populations with
TFR of 5 and West life table level 5. The light grey bands represent the
interval between the 2.5% and 97.5% quantiles; the dark grey bands represent
the interval between the 25% and 75% quantiles; and the white lines represent
medians.

For a given combination of settings, the only source of variability across
populations is random variation in demographic events. The width of the grey
bands therefore shows the impact on the population indicators of random
variation in demographic events. Comparing widths across columns shows how
sensitivity to random variation declines with population size.

The key finding from Figure 3 is that, with this particular set of parameters,
variability falls sharply between populations with initial size $100$ and
populations with initial size $1,000$, and falls less sharply between
populations with initial size $1,000$ and initial size $10,000$. Populations
of around $100$ can have vastly different values for the four summary
indicators, even when their underlying fertility and mortality rates are the
same. Populations of around $10,000$ show some differences in the summary
indicators, but these differences are relatively small. Populations of around
$1,000$ behave more like populations of $10,000$ than like populations of
$10,000$.

###  3.3 Distributions of summary indicators for all combinations of total
fertility rate and life table level

Figures 5–8 in the Appendix show results equivalent to those of Figure 3 for
other combinations of TFR and West life table levels. Comparing across the
five graphs is, however, difficult. To simplify comparisons, we restrict our
attention to a single post-transition year. Figure 4 focuses on year 100, and
gives results for five combinations of TFR and life table level. The tops and
bottoms of the boxes in Figure 4 correspond to the 25% and 75% quantiles in
Figure 3; the tops and bottoms of the vertical lines correspond to the tops
and bottoms of the 2.5% and 97.5% quantiles in Figure 3; and the horizontal
lines correspond to the medians.

Figure 4: Distributions of summary indicators at year 100 of the simulations,
for every combination of population size, fertility rates, and mortality
rates. The tips of the vertical lines represent 2.5% and 97.5% quantiles; the
tops and bottoms of the boxes represent 25% and 75% quantiles; and the
horizontal lines represent medians.

In Figure 4, like in Figure 3, there is a large decline in variability when
moving from initial population size $100$ to initial size population size
$1,000$. There is also a decline in variability when moving from initial
population size $1,000$ to initial population size $10,000$, but it is more
modest.

With initial population size $100$, random variability in the indicators is
sufficiently large, relative to differences across the five demographic
scenarios, that the distributions largely overlap. The distribution of the
percent aged 0–14 under a TFR of 5 and a life table level of 1, for instance,
overlaps with the distribution of the percent aged 0–14 under a TFR of 5 and a
life table level of 5.

With initial population size $10,000$, random variability is much less, and
the distributions of the indicators are accordingly much more distinct. There
is little overlap, for instance, between the distributions generated by a TFR
of 5 and a life table level of 1 and the distributions generated by a TFR of 5
and a life table level of 9.

The amount of overlap with initial population size $1,000$ is again
intermediate between the other two sizes. However, the results are again
closer to those for size $10,000$ than to size $100$.

##  4 Discussion

Our simulation results imply that, under the types of high-fertility, high-
mortality regimes typical of historical populations, demographic analyses of
very small populations, on the order of 100 people, need to be highly
attentive to random variation. It is risky to use differences in age structure
to make inferences about differences in underlying demographic rates, because
the differences in age structure could easily be due to random variation.

Demographic analyses of very small populations need to model this random
variation. One way of doing so when data from multiple small populations is
available is to use Bayesian hierarchical models that pool information across
the populations, and that allow information about plausible ranges to be
incorporated in a transparent and systematic way. Schmertmann and Hauer,
(2019) is an example of this approach. A second, complementary, approach is to
build population simulations that explicitly model the sources of variation.
The simulation methods presented in this paper are one way of doing so.

Our simulation results also imply that once population sizes reach 1,000 or
so, variation attributable to randomness in births and deaths has already
reached modest levels, not much greater than would be expected for much larger
populations. Demographers dealing with populations of this size do, of course,
need to be careful not to over-interpret small differences in age-sex
structure. But it appears that variability that is specifically due to
randomness in births and deaths is less of an obstacle to demographic analyses
of these populations than might be thought. Demographers working with broad
measures of age-sex structure in historical populations would probably be
better to devote scarce analytical resources to assessing the impact of issues
such as measurement error or variation in underlying demographic rates than to
individual-level random variation in demographic events.

Demographers working with broad measures of age-sex structure in populations
of 10,000 or high can, like demographers working with national-level
populations, safely ignore individual-level variation in demographic events.

## References

  * Booth, (1984) Booth, H. (1984).  Transforming Gompertz’s function for fertility analysis: The development of a standard for the relational Gompertz function.  Population Studies, 38(3):495–506.
  * Caldwell et al., (1987) Caldwell, J., Caldwell, P., Caldwell, B., Eversley, D., Fix, A. G., Howell, N., Nag, M., Taeuber, C., and van de Walle, F. (1987).  Anthropology and demography: The mutual reinforcement of speculation and research [and comments and reply].  Current anthropology, 28(1):25–43.
  * Caldwell, (1990) Caldwell, J. C. (1990).  Cultural and social factors influencing mortality levels in developing countries.  The Annals of the American Academy of Political and Social Science, 510(1):44–59.
  * Caswell, (2001) Caswell, H. (2001).  Matrix Population Models: Construction, Analysis, and Interpretation. Second Edition. Sinauer Associates, Sunderland.
  * Coale et al., (1983) Coale, A., Demeny, P., and Vaughan, B. (1983).  Regional model life tables and stable populations. 2nd ed.  Academic Press.
  * Coale and Watkins, (1986) Coale, A. J. and Watkins, S. C. (1986).  The Decline of Fertility in Europe.  Princeton University Press, Princeton.
  * Cook and Borah, (1971) Cook, S. F. and Borah, W. (1971).  Essays in Population History. Volume 1: Mexico and the Carribean.  University of California Press.
  * Covey et al., (2011) Covey, R. A., Childs, G., and Kippen, R. (2011).  Dynamics of indigenous demographic fluctuations: Lessons from sixteenth-century cusco, peru.  Current Anthropology, 52(3):335–360.
  * Dobyns, (1966) Dobyns, H. F. (1966).  Estimating aboriginal american population: An appraisal of techniques with a new hemispheric estimate.  Current Anthropology, 7(4):395–416.
  * Gelman and Vehtari, (2021) Gelman, A. and Vehtari, A. (2021).  What are the most important statistical ideas of the past 50 years?  Journal of the American Statistical Association, pages 1–29.
  * Henry, (1968) Henry, L. (1968).  Historical demography.  Daedalus, pages 385–396.
  * Jones, (2007) Jones, J. H. (2007).  demogR: A package for the construction and analysis of age-structured demographic models in R.  Journal of Statistical Software, 22(10):1–28.
  * Keyfitz and Caswell, (2005) Keyfitz, N. and Caswell, H. (2005).  Applied mathematical demography. Third Edition.  Springer.
  * Lee, (1998) Lee, R. D. (1998).  Probabilistic approaches to population forecasting.  Population and Development Review, 24:156–190.
  * Moultrie et al., (2013) Moultrie, T., Dorrington, R., Hill, A., Hill, K., Timæus, I., and Zaba, B. (2013).  Tools for Demographic Estimation.  International Union for the Scientific Study of Population, Paris.
  * Petersen et al., (1975) Petersen, W., Braidwood, R., Dobyns, H., Eberhard, W., Kennedy, R., Kurth, G., Meiklejohn, C., Onuki, Y., Weiss, K., and Wilkinson, P. (1975).  A demographer’s view of prehistoric demography.  Current Anthropology, 16(2):227–245.
  * Preston et al., (2001) Preston, S., Heuveline, P., and Guillot, M. (2001).  Demography: Modelling and Measuring Population Processes.  Blackwell, Oxford.
  * Roberts, (1989) Roberts, L. (1989).  Disease and death in the new world: Historians now agree that the european discovery of the americas touched off waves of epidemics, but a debate is raging over the size of pre-columbian populations.  Science, 246(4935):1245–1247.
  * Rosenthal, (1997) Rosenthal, P.-A. (1997).  Thirteen years of debate: From population history to French historical demography (1945-1958).  Population an English Selection, pages 215–241.
  * Schmertmann and Hauer, (2019) Schmertmann, C. P. and Hauer, M. E. (2019).  Bayesian estimation of total fertility from a population’s age–sex structure.  Statistical Modelling, 19(3):225–247.
  * Thornton, (1980) Thornton, R. (1980).  Recent estimates of the prehistoric california indian population.  Current Anthropology, 21(5):702–704.
  * Van Imhoff and Post, (1998) Van Imhoff, E. and Post, W. (1998).  Microsimulation methods for population projection.  Population: An English Selection, pages 97–138.
  * Wachter et al., (1978) Wachter, K. W., Hammel, E. A., and Laslett, P. (1978).  Statistical studies of historical social structure.  Elsevier.
  * Wrigley and Schofield, (1989) Wrigley, E. A. and Schofield, R. S. (1989).  The population history of England 1541-1871.  Cambridge University Press.
  * Zagheni, (2015) Zagheni, E. (2015).  Microsimulation in demographic research.  International encyclopedia of social and behavioral sciences, 15:343–346.

##  Appendix A Year-to-year changes in observed death rates in populations

Here we use some simple simulations to support the claim from the Introduction
that the number of deaths can easily half or double from one year in a small
population. We use a population size of 100, we assume that deaths follow a
Poisson distribution, and we ignore the effect of the annual number of deaths
on the annual population at risk, which is small in percentage terms. We
simulate deaths in pairs of years, and calculate the proportion of pairs in
which the ratio of deaths in one year to deaths in the other year is less than
0.5 or greater than 2.

    percent_halve_or_double <- function(popn_size, death_rate) {
        ## Calculate expected number of births
        lambda <- popn_size * death_rate
        ## Generate ’n’ births in first and second years
        deaths_first_year <- rpois(n = 1000000, lambda = lambda)
        deaths_second_year <- rpois(n = 1000000, lambda = lambda)
        ## Obtain the ratio of births in year 1 to births in year 2,
        ## excluding cases where births in both years are 0
        ratio <- deaths_second_year / deaths_first_year
        combined_deaths <- deaths_first_year + deaths_second_year
        ratio[combined_deaths == 0] <- 1
        ## Calculate the proportion of cases where births
        ## in year 2 where less than half, or more than twice,
        ## births in year 1
        ans <- mean((ratio < 0.5) | (ratio > 2))
        ## Convert to percent
        round(ans * 100)
    }

    set.seed(0)
    percent_halve_or_double(popn_size = 100, death_rate = 0.02)
    ## [1] 44
    percent_halve_or_double(popn_size = 100, death_rate = 0.03)
    ## [1] 37
    percent_halve_or_double(popn_size = 100, death_rate = 0.04)
    ## [1] 31

##  Appendix B Calculating mortality rates for Lexis triangles for age 0–4

We calculate life table mortality rates by Lexis triangle, for age group 0–4
and sex $s$, by first calculating person-years lived and deaths for each
triangle. The total number of person-years lived within age group 0 over a
5-year period is $5L_{0s}$. The total number of deaths is
$5L_{\text{0}s}m_{\text{0}s}$. Similarly, the total number of person-years
lived and deaths for age group 1–4 is $5L_{\text{1-4},s}$ and
$5L_{\text{1-4},s}m_{\text{1-4},s}$. We allocate person-years lived according
to the degree of overlap between each Lexis triangle and age group 0, depicted
as a grey rectangle in Figure 1. The lower Lexis triangle thus receives 4.5/5
of the person-years lived within age group 0, and the upper Lexis triangle
receives 0.5/5 person years. Deaths are allocated the same way. This yields
the formulas

$\displaystyle m_{\text{0--4},s,\text{Low}}$ $\displaystyle=\frac{\text{Deaths
in lower triangle}}{\text{Person-years lived in lower triangle}}$ (1)
$\displaystyle=\frac{(4.5/5)\times 5L_{\text{0},s}m_{\text{0},s}+(8/20)\times
5L_{\text{1--4},s}m_{\text{1--4},s}}{(4.5/5)\times
5L_{\text{0},s}+(8/20)\times 5L_{\text{1--4},s}}$ (2)
$\displaystyle=\frac{4.5L_{\text{0},s}m_{\text{0},s}+2L_{\text{1--
4},s}m_{\text{1--4},s}}{4.5L_{\text{0},s}+2L_{\text{1--4},s}}$ (3)
$\displaystyle m_{\text{0--4},s,\text{Up}}$ $\displaystyle=\frac{\text{Deaths
in upper triangle}}{\text{Person-years lived in upper triangle}}$ (4)
$\displaystyle=\frac{(0.5/5)\times 5L_{\text{0},s}m_{\text{0},s}+(12/20)\times
5L_{\text{1--4},s}m_{\text{1--4},s}}{(0.5/5)\times
5L_{\text{0},s}+(12/20)\times 5L_{\text{1--4},s}}$ (5)
$\displaystyle=\frac{0.5L_{\text{0},s}m_{\text{0},s}+3L_{\text{1--
4},s}m_{\text{1--4},s}}{0.5L_{\text{0},s}+3L_{\text{1--4},s}}.$ (6)

##  Appendix C Plots of simulations over time

Figure 5: Distributions of summary indicators, over time, for populations with
TFR 5 and West life table level 1. Figure 6: Distributions of summary
indicators, over time, for populations with TFR of 4 and West life table level
5. Figure 7: Distributions of summary indicators, over time, for populations
with TFR of 6 and West life table level 5. Figure 8: Distributions of summary
indicators, over time, for populations with TFR of 5 and West life table level
9.