11institutetext: Wellcome/EPSRC Centre for Interventional and Surgical
Sciences (WEISS), University College London, Charles Bell House, 43–45 Foley
Street, London W1W 7TY, UK
11email<EMAIL_ADDRESS>
# HUP-3D: A 3D multi-view synthetic dataset for assisted-egocentric hand-
ultrasound pose estimation

Manuel Birlo  Razvan Caramalau  Philip J. “Eddie” Edwards  Brian Dromey
Matthew J. Clarkson  Danail Stoyanov

###### Abstract

We present HUP-3D, a 3D multi-view multi-modal synthetic dataset for hand-
ultrasound (US) probe pose estimation in the context of obstetric ultrasound.
Egocentric markerless 3D joint pose estimation has potential applications in
mixed reality based medical education. The ability to understand hand and
probe movements programmatically opens the door to tailored guidance and
mentoring applications. Our dataset consists of over 31k sets of RGB, depth
and segmentation mask frames, including pose related ground truth data, with a
strong emphasis on image diversity and complexity. Adopting a camera
viewpoint-based sphere concept allows us to capture a variety of views and
generate multiple hand grasp poses using a pre-trained network. Additionally,
our approach includes a software-based image rendering concept, enhancing
diversity with various hand and arm textures, lighting conditions, and
background images. Furthermore, we validated our proposed dataset with state-
of-the-art learning models and we obtained the lowest hand-object keypoint
errors. The dataset and other details are provided with the supplementary
material. The source code of our grasp generation and rendering pipeline will
be made publicly available.

###### Keywords:

Egocentric 3D joint hand and tool pose estimation Synthetic datasets
Obstetrics ultrasound.

##  1 Introduction

The ability to infer hand and tool pose information from video data in
clinical setups opens the door to several potentially useful applications
aimed at assisting clinicians through context-specific evaluation of their
movement behavior, and respective novel methods focusing on marker-free video-
based clinical skill assessment were already proposed such as surgical hand
and tool pose estimation [24], surgical tool movement analysis [25] and skill
assessment in robotic surgery [26]. Furthermore, solutions for skill
assessment that rely on physical motion sensors have been developed in fields
such as hand motion analysis for endovascular procedures [24] and guidance for
using ultrasound (US) probes in obstetric US [25]. With the emerging trend of
mixed reality head mounted displays, dominated by the Microsoft HoloLens 2
111https://www.microsoft.com/en-us/hololens, egocentric pose estimation
methods arose, such as probe tracking for US-guided procedures [27].

In the context of obstetric US, we explore 3D joint hand-tool pose estimation
with an eye toward applications in mixed reality-based medical education,
where analysis of hand and probe movements could facilitate holographic
assistive probe guidance. Such innovations aim to help standardize clinical
training protocols, addressing the current lack of universally accepted
competence measures for effective US probe guidance, despite the existence of
standardized target scanning US planes [18]. By distinguishing between novice
and expert clinician movements, machine learning-based pose estimation emerges
as a powerful tool for developing standardized training approaches. This
technique not only aligns with but also supports established clinical
practices for fetal weight estimation through biometry, offering a pathway to
more uniform and effective clinical training methodologies [17].

Image datasets required for machine learning-based model training and
subsequent pose estimation can be categorized into real and synthetic images.
Real dataset generating methods often employ marker-free methods to capture
hand grasp information directly from image data [6, 7, 8, 9, 10]. While real
images offer the advantage of authentic contexts [10, 13], they pose
challenges for generating accurate ground truth due to the labor-intensive
nature of manual annotations and potential biases from sensor use [7]. Marker-
free approaches, by avoiding markers on the tool and/or hand, mitigate the
risk of pose prediction bias, making them a popular choice for reducing
inaccuracies associated with additional visible sensors or markers [8, 7].

Synthetic images, however, offer built-in ground truth from 3D models,
simulating realistic grasping scenarios with the benefits of easy scalability
and generalizability to real images [6]. Furthermore, synthetic ground truth
proves useful for addressing mutual occlusions resulting from hand-tool
interactions.

When creating training images for pose estimation, it’s crucial to account for
the dataset’s generalizability, particularly the expected camera location and
range of viewpoints. Applications vary, with some utilizing non-egocentric
views for capturing hand grasps [6, 10], while mixed and augmented reality
setups, especially those using head-mounted devices like the Microsoft
HoloLens 2222https://www.microsoft.com/en-us/hololens/, require considering
egocentric perspectives for pose estimation from device-recorded camera data
[7, 8]. When capturing clinical instruments like US probes in synthetic
images, the realism of hand grasps is constrained by specific hand-tool
contact areas and orientations. This requirement, previously noted in the
context of orthopedic surgical tools such as drills [7] and other instruments
like scalpels and diskplacers [8], poses a challenge. Traditional grasp
generation techniques, like those enabled by robotic grasping software [14]
and used in similar studies [6, 7], encounter difficulties due to the specific
dimensions and clinically relevant grasp positions of the US probe.
Consequently, we employed more flexible solutions such as a generative model
for machine learning-based grasp generation [9], which has been validated in
clinical environments [8].

To increase image diversity for egocentric applications in a scalable way, we
broadened the approach to multi-view by allowing camera movement around a
sphere’s surface, centered on the hand. This method supports both egocentric
distances and a mix of egocentric and non-egocentric viewpoints.

Our contributions can be summarized as follows:

  * •

A scalable synthetic multi-modal (RGB-D, segmentation maps) image generation
pipeline capable of producing a wide variety of realistic hand-ultrasound
probe grasp frames, independent of prior external data recording requirements

  * •

A novel sphere-based camera viewpoint concept that enhances frame
generalizability by combining egocentric head-hand distances with non-
egocentric camera viewpoints.

  * •

HUP-3D: A diverse multimodal synthetic dataset tailored for joint 3D hand and
tool pose estimation, featuring the Voluson™
C1-5-D333https://services.gehealthcare.com/gehcstorefront/p/5499513 ultrasound
probe commonly used in obstetrics, including a variety of hand poses,
textures, backgrounds, lighting, and camera angles

  * •

Lowest hand and object 3D pose estimation errors for a synthetic dataset with
a trained state-of-the-art model, HOPE-net [30]

##  2 Method

Figure 1: Grasp Generation (blue) and Rendering Pipeline (red): The process
begins with a MANO hand model initialization and a BPS-encoded Voluson model
point cloud. CoarseNet generates initial hand poses, further refined by
RefineNet for precise hand-probe alignment. In the rendering phase, the
optimized hand pose, model vertices, and a SMPL-H model are processed in
Blender. Using a multi-viewpoint camera via a spherical layout and centered on
the hand and arm, several textures and backgrounds are applied for diverse
RGB-D, segmentation maps, and annotations.

Focusing on potential medical education applications in the context of
obstetric US, but with a high degree of flexibility towards other use cases,
as described in section 1, our synthetic image generation pipeline is split
into two sections: grasp generation and grasp rendering. These are described
in detail in the following subsections. A graphical overview of our pipeline
is shown in Fig. 1.

###  2.1 Grasp generation

To streamline ground truth annotation, we adopted a strategy focused on
generating synthetic grasp images, avoiding the complexities associated with
annotating real images. This approach allowed us to maintain a clear and
manageable rendering workflow. The underlying idea in pursuing a purely
synthetic image generation approach is to explore the possibility of
generating a sufficiently large variety of training images that allow
acceptable generalizability to real images in joint 3D hand and tool pose
prediction. Our initial feasibility study incorporated the use of the robotic
grasping tool [14] already mentioned in section 1, which turned out to be
error-prone and did not produce a sufficiently large variety of plausible hand
grasps due to restrictions in terms of hand dexterity. Therefore, we adapted
the generative model proposed in [9] for joint 3D grasp generation to a more
clinical scenario. Our grasp generation process employs two sequential
networks based on the MANO hand model [28]: an encoder-decoder network that
generates initial coarse hand poses and a subsequent neural network dedicated
to fine-tuning these poses, specifically enhancing accuracy in hand-tool
interaction areas. The encoder which samples from a normal distributed 16
dimensional latent space, requires encoded point cloud representations [19] of
the probe model ($\Omega_{BPS}$), together with the MANO right hand model’s
initial translation $\gamma\in\mathbb{R}^{3}$ and hand wrist orientation
$\theta_{wrist}\in\mathbb{R}^{3}$. Defined Euler angles $\Theta_{Vol}$ for
probe meshes $\Omega_{BPS}$ were used for precise grasp pose control.
Originally, the model described in [9] was trained with ordinary objects (like
mugs, cameras etc.). However, we extended its capability to the Voluson US
probe The decoder outputs an initial hand pose
$[\gamma,\theta_{full\\_pose}]$, which is subsequently refined through a
neural network utilizing the vertices of the probe model $\Omega_{Vert}$ and
the vertex distances $\Delta^{hand}_{obj}$ between hand and probe. This
refined pose, expressed as $\Psi:=[\gamma^{**},\theta_{full\\_pose}]$, forms
the foundation for our grasp rendering approach detailed in Sec. 2.2.

(a)

(b)

Figure 2: (a) Schematic grasp conversion from generative model to rendering
software, including probe offset ($\Delta z$) correction. (b) Grasp rendering
overview: (1) SMPL-H body model grasping the probe, showing egocentric and
non-egocentric views. (2) Right arm and sphere-based camera orientations with
remaining SMPL-H body parts hidden. (3) Camera angle sphere concept with views
at various latitudes, centered on hand mesh; defines sphere ($r_{sphr}$) and
circle ($r_{circ}$) radii. (4) Rendered hand-probe scene example from a sphere
camera position. 33footnotetext: https://www.meshlab.net/

###  2.2 Grasp rendering

Building on the concept of utilizing the Blender [15] open-source 3D graphics
software for grasp rendering, as demonstrated in [6, 7], we tailored our
rendering pipeline to accommodate the grasp poses $\Psi$ produced by the
generative model outlined in Sec. 2.1. Additionally, this rendering approach
incorporates a SMPL-H body model [16], a MANO right hand model $M_{Vert}$, and
the probe model’s vertex data $\Omega_{Vert}$. Details of the grasp rendering
part can be found in the lower part of Fig. 1. Our methodology positions the
probe at the origin of the rendering software’s internal world coordinate
system, yet the modified hand grasp pose $\Psi$ necessitates an offset for the
probe model’s translation $\gamma_{\omega}$, specifically a displacement of
$\Delta z$ along the positive z-axis, which was calculated programmatically
through polygon offset analysis. Subsequently, we adjusted for the probe’s
translation offset via $\gamma^{\Delta}_{\omega}=\gamma_{\omega}+(0,0,-\Delta
z)$. Fig. 2a schematically illustrates the probe’s z-offset correction. Fig.
2a schematically illustrates the probe’s z-offset correction.

To enhance the diversity of generated grasp images across different camera
perspectives, we transitioned from the purely egocentric viewpoint strategy to
a sphere-based methodology outlined in Sec. 2.2.1, illustrated in Figs. 1 and
2b. This helps a pose estimation model’s capability to accurately predict
grasp poses amid challenges like mutual occlusions between hands and tools.

For each grasp produced by our module, we generate a synthetic image for each
camera view angle $\Theta_{k}\in\\{\Theta_{1},\ldots,\Theta_{N}\\}$, covering
$N$ positions around the sphere. Our render scene setup utilizes two shades of
clinical gloves, varied scene lighting, and eight backgrounds (a lab with a
SPACE-FAN ultrasound fetus
model444https://www.kyotokagaku.com/en/products_data/us-7_en/, consultation
rooms, a white background, and real abdomens of pregnant women). The rendering
model outputs a comprehensive set of images for each grasp, including RGB-D
and segmentation maps, as well as ground truth annotations. Rendered sample
frames from the HUP-3D dataset are shown in Fig. 3.

Figure 3: Sample frames from the HUP-3D dataset, grouped columnwise, from left
to right: RGB, depth, segmentation map, and ground truth annotations.

####  2.2.1 Camera view angle sphere concept

Our methodology diverges from traditional egocentric viewpoints by
implementing a sphere-based camera view setup to capture both egocentric and
non-egocentric images, enhancing dataset diversity. This method, inspired by
[3], involves distributing camera positions around a sphere, creating a varied
perspective landscape around the right hand. To systematize camera placement,
the sphere is divided into horizontal segments, determined by latitude angles,
to evenly distribute viewpoints. Specifically, the number of latitude segments
$N_{\text{latitude\\_floors}}$ and circles per segment
$N_{\text{circles}}^{(i)}$ are calculated to ensure comprehensive coverage:

$N_{\text{latitude\\_floors}}=\left\lfloor\frac{\pi}{2\arcsin\left(\frac{r_{circ}}{r_{sph}}\right)}\right\rfloor,\quad
N_{\text{circles}}^{(i)}=\left\lfloor\frac{2\pi
r_{sph}\sin(\theta_{i})}{2r_{circ}}\right\rfloor$ (1)

For each segment $i$ and each circle $j$ within, camera locations are defined
by their spherical coordinates $(\theta_{i},\phi_{j}^{(i)})$, ensuring a
uniform spread of angles:

$(\theta_{i},\phi_{j}^{(i)})\quad\text{with}\quad\phi_{j}^{(i)}=j\cdot\frac{2\pi}{N_{\text{circles}}^{(i)}},\quad
j=0,1,\dots,N_{\text{circles}}^{(i)}-1$ (2)

This structured approach facilitates the generation of camera angles
$\Theta_{k}$, utilized in our subsequent rendering process. Figs. 1 and 2b
visually demonstrate this concept, showcasing the strategic camera placement
and the diverse grasp views it enables.

####  2.2.2 Dataset comparison

In Table 1, we enlist the top clinical and non-clinical datasets, together
with their properties. HUP-3D is the largest multi-view dataset for clinical
applications, presenting 3 possible modalities, RGB-DS (color, depth and
segmentation maps). Only POV-Surgery [8] contains a higher number, but with
less samples per tool (29k) and just firs-person view.

Dataset | # frames |  | Source
---
(Real/ Synth)
| Viewpoints
---
(Single/Multi/Ego)
Annotations | Modalities |  | Clinical
---
(no. of tools)
HO-3D [10] | 77.5k | Real | Single | automatic | RGB | -
ObMan [6] | 153k | Synth | Multi | automatic | RGB-DS | -
ContactPose [11] | 2.9M | Real | Multi | semi-automatic | RGB-D | -
Hein et al. [7] | 10.5k | Synth | Ego | automatic | RGB-DS | 1
POV-Surgery [8] | 88k | Synth | Ego | automatic | RGB-DS | 3
HUP-3D (ours) | 31680 | Synth | Multi | automatic | RGB-DS | 1

Table 1: Dataset comparison: HUP-3D outstands as the first multi-view 3D
hand-(clinical)object dataset.

##  3 Experiment

To support the utility of our proposed dataset HUP-3D, we deploy a deep
learning (DL) state-of-the-art model designed for other datasets like HO-3D
[10]. As mentioned before, our dataset consists of 31,680 image sets from 11
realistic hand-object grasps. In a supervised learning setting, we further
split the data as 7 grasps for training (20,160), 2 grasps for validation and
2 more for testing (5,760). This will ensure the generalisation capability of
the tested DL model.

###  3.1 3D hand-probe pose estimation

There has been extensive DL methods proposed for 3D hand-object pose
estimation in the computer vision community. One of these competitive
baselines is HOPE-net [30], originally tested on real data. HOPE-net extends
the capabilities of residual convolutional neural networks [29] with an
adaptive Graph U-Net module [31]. This module manages to reduce the highly
non-linear regression of the 3D hand and object coordinates.

Therefore, the task of the estimator is to map the RGB images to 3D world
coordinates of the hand skeleton and the object’s boundary corners. For this,
we minimise a mean square error loss during training on both 2D and 3D
coordinates. In terms of training, we follow the same settings as in the
original HOPE-net paper [30].

Quantitatively, once trained, we measure the error in millimeters between the
predicted joints and the ground-truth. On the testing set, we obtain a total
error of 8.65 mm as a mean per joint position error (MPJPE) with 5.33 mm from
the hand and 17.05 mm from the object. The testing error is the lowest
compared to the other clinical datasets like POV-Surgery [8] (14.35 mm) and
Hein et al. [7] (17.02 mm) where even more advanced DL models were used. This
proves that a multi-viewpoint dataset helps to estimate the object and the
hand location with higher precision. In our experiments, we have also tested a
more simple baseline composed just of ResNet-50 [29], where the error was
higher at 9.69 mm. In Fig. 4, we visually confirm the accurate predicted
keypoints from our HUP-3D test set.

Figure 4: Qualitative results, shown with 4 test images from HUP-3D: image
columns from left to right: RGB, predicted hand joints, predicted probe
corners, predicted joints and corners, ground truth of joints and corners

##  4 Conclusion and future work

We introduced HUP-3D, a pioneering multi-view 3D hand-object dataset tailored
for obstetrics ultrasound hand-probe grasps. HUP-3D aims to enhance research
in clinical movement analysis via egocentric camera and mixed reality
applications. Our data generation process leverages a versatile model for
grasp generation and an efficient automated rendering pipeline, illustrating
the benefits of our multi-view camera sphere approach. A baseline model
evaluation confirmed our method’s effectiveness, even with significant hand-
probe occlusions. Future efforts will focus on improving real-world
applicability by incorporating automatically annotated real images and
developing more sophisticated grasp generation techniques that incorporate
temporal sequences for better manual interaction and predictions.

## References

  * [1] Yana Hasson, Bugra Tekin, Federica Bogo, Ivan Laptev, Marc Pollefeys, and Cordelia Schmid.  Leveraging Photometric Consistency Over Time for Sparsely Supervised Hand-Object Reconstruction.  In Proceedings of the IEEE/CVF CVPR, June 2020.
  * [2] H. Jiang, S. Liu, J. Wang, and X. Wang, “Hand-Object Contact Consistency Reasoning for Human Grasps Generation,” in Proceedings of the ICCV, 2021.
  * [3] Abdulkadir Akin, E. Erdede, Hossein Afshari, Alexandre Schmid, and Yusuf Leblebici.  Enhanced Omnidirectional Image Reconstruction Algorithm and Its Real-Time Hardware.  In Proceedings - 15th Euromicro Conference on Digital System Design, DSD 2012, Sep 2012.  10.1109/DSD.2012.52.  ISBN 978-1-4673-2498-4.
  * [4] Shaowei Liu, Hanwen Jiang, Jiarui Xu, Sifei Liu, and Xiaolong Wang.  Semi-Supervised 3D Hand-Object Poses Estimation With Interactions in Time.  In Proceedings of the IEEE/CVF CVPR, pages 14687–14697, June 2021.
  * [5] Markus Oberweger, Paul Wohlhart, and Vincent Lepetit.  Generalized Feedback Loop for Joint Hand-Object Pose Estimation.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 42(8):1898–1912, 2020.  10.1109/TPAMI.2019.2907951.
  * [6] Y. Hasson, G. Varol, D. Tzionas, I. Kalevatykh, M. J. Black, I. Laptev, and C. Schmid, "Learning joint reconstruction of hands and manipulated objects," in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2019.
  * [7] Hein, J., Seibold, M., Bogo, F., Farshad, M., Pollefeys, M., Fürnstahl, P. and Navab, N. (2021). Towards markerless surgical tool and hand pose estimation. International Journal of Computer Assisted Radiology and Surgery, 16(5), 799–808. https://doi.org/10.1007/s11548-021-02369-2
  * [8] R. Wang, S. Ktistakis, S. Zhang, M. Meboldt, and Q. Lohmeyer, "POV-Surgery: A Dataset for Egocentric Hand and Tool Pose Estimation During Surgical Activities," in Proceedings of the International Conference on Medical Image Computing and Computer-Assisted Intervention, 2023, pp. 440–450.
  * [9] O. Taheri, N. Ghorbani, M. J. Black, and D. Tzionas, “GRAB: A Dataset of Whole-Body Human Grasping of Objects,” in Proceedings of the European Conference on Computer Vision (ECCV), 2020. [Online]. Available: https://grab.is.tue.mpg.de
  * [10] Shreyas Hampali, Mahdi Rad, Markus Oberweger, and Vincent Lepetit.  HOnnotate: A Method for 3D Annotation of Hand and Object Poses.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2020.
  * [11] Samarth Brahmbhatt, Chengcheng Tang, Christopher D. Twigg, Charles C. Kemp, and James Hays.  ContactPose: A Dataset of Grasps with Object Contact and Hand Pose.  In Andrea Vedaldi, Horst Bischof, Thomas Brox, and Jan-Michael Frahm, editors, Computer Vision – ECCV 2020, pages 361–378. Springer International Publishing, Cham, 2020.  ISBN 978-3-030-58601-0.
  * [12] Bardia Doosti, Shujon Naha, Majid Mirbagheri, and David J. Crandall.  HOPE-Net: A Graph-Based Model for Hand-Object Pose Estimation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), June 2020.
  * [13] T. Kwon, B. Tekin, J. Stühmer, F. Bogo, and M. Pollefeys, “H2O: Two hands manipulating objects for first person interaction recognition,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pp. 10138–10148, 2021.
  * [14] Andrew T. Miller and Peter K. Allen.  Graspit! A Versatile Simulator for Robotic Grasping.  IEEE Robotics & Automation Magazine, 11(4):110–122, 2004.  IEEE.
  * [15] Blender Online Community.  Blender - a 3D modelling and rendering package.  Stichting Blender Foundation, Amsterdam, 2018.
  * [16] Javier Romero, Dimitrios Tzionas, and Michael J. Black.  Embodied hands: modeling and capturing hands and bodies together.  ACM Transactions on Graphics, 36(6):245, November 2017.  Association for Computing Machinery, New York, NY, USA.  ISSN 0730-0301.
  * [17] B. P. Dromey, S. Ahmed, F. Vasconcelos, E. Mazomenos, Y. Kunpalin, S. Ourselin, J. Deprest, A. L. David, D. Stoyanov, and D. M. Peebles, “Dimensionless squared jerk: An objective differential to assess experienced and novice probe movement in obstetric ultrasound,” _Prenatal Diagnosis_ , vol. 11, 2020.
  * [18] Y. Cai, R. Droste, H. Sharma, P. Chatelain, L. Drukker, A. T. Papageorghiou, and J. A. Noble, “Spatio-temporal visual attention modelling of standard biometry plane-finding navigation,” _Medical Image Analysis_ , vol. 65, 2020.
  * [19] Sergey Prokudin, Christoph Lassner, and Javier Romero.  Efficient learning on point clouds with basis point sets.  In _Proceedings of the IEEE/CVF ICCV_ , 2019.
  * [20] G. Varol, J. Romero, X. Martin, N. Mahmood, M.J. Black, I. Laptev, and C. Schmid.  Learning from synthetic humans.  In: _Proceedings of the IEEE CVPR_ , 2017.
  * [21] D. P. Azari, Y. H. Hu, B. L. Miller, B. V. Le, and R. G. Radwin, “Using surgeon hand motions to predict surgical maneuvers,” _Human Factors_ , vol. 61, 2019, SAGE Publications Sage CA: Los Angeles, CA.
  * [22] X.-H. Zhou, G.-B. Bian, X.-L. Xie, Z.-G. Hou, X. Qu, and S. Guan, “Analysis of Interventionalists’ Natural Behaviors for Recognizing Motion Patterns of Endovascular Tools During Percutaneous Coronary Interventions,” _IEEE Transactions on Biomedical Circuits and Systems_ , vol. 13, 2019.
  * [23] R. Droste, L. Drukker, A. T. Papageorghiou, and J. A. Noble, “Automatic probe movement guidance for freehand obstetric ultrasound,” in _MICCAI 2020: 23rd International Conference_ , Springer, 2020.
  * [24] E. D. Goodman, K. K. Patel, Y. Zhang, W. Locke, C. J. Kennedy, R. Mehrotra, S. Ren, M. Y. Guan, M. Downing, H. W. Chen, et al., “A real-time spatiotemporal AI model analyzes skill in open surgical videos,” _arXiv preprint arXiv:2112.07219_ , 2021.
  * [25] A. Jin et al., “Tool Detection and Operative Skill Assessment in Surgical Videos Using Region-Based Convolutional Neural Networks,” in _2018 IEEE WACV_ , 2018.
  * [26] G. Lajkó, R. Nagyné Elek, and T. Haidegger, “Endoscopic image-based skill assessment in robot-assisted minimally invasive surgery,” _Sensors_ , vol. 21, MDPI.
  * [27] T. Nguyen, W. Plishker, A. Matisoff, K. Sharma, and R. Shekhar, “HoloUS: Augmented reality visualization of live ultrasound images using HoloLens for ultrasound-guided procedures,” _International Journal of Computer Assisted Radiology and Surgery_ , Springer vol. 17, 2022.
  * [28] J. Romero, D. Tzionas, and M. J. Black, “Embodied Hands: Modeling and Capturing Hands and Bodies Together,” _ACM Transactions on Graphics, (Proc. SIGGRAPH Asia)_ , Nov. 2017.
  * [29] Kaiming He and Xiangyu Zhang and Shaoqing Ren and Jian Sun, “Deep Residual Learning for Image Recognition”, 2015, Tech Report, eprint=1512.03385.
  * [30] B. Doosti, S. Naha, M. Mirbagheri and D. Crandall, “HOPE-Net: A Graph-based Model for Hand-Object Pose Estimation”, (CVPR), June, 2020.
  * [31] H. Gao and S. Ji, “Graph U-Nets", ICML, 2019

##  1 Supplementary material

####  1.0.1 HUP-3D Dataset

Our HUP-3D dataset can be downloaded via the following google drive link:
HUP-3D

####  1.0.2 Grasp generation and rendering pipeline configuration

Table 1 outlines the different parameters our frame generation pipeline
utilizes, which influence the diversity and quantity of frames produced:

Table 1: Frame generation configuration Parameter |  Description |  Values
---|---|---
$(\Theta_{Vol})_{j}$ |  k-th ($k\in[1,..,N]$) probe rotation Euler angles, used within GrabNet |  Manually selected set of three Euler angles $[\alpha,\beta,\gamma|\alpha,\beta,\gamma\in[0,...,360]]$
$(G_{[\Theta_{Vol}]})_{k}$ |  k-th manually selected GrabNet-based grasp out of all grasps generated by $\Theta_{Vol}$ |  11 manually selected plausible grasps via visual inspection in MeshLab
$r_{sph}$ |  Sphere radius |  0.8m
$r_{circ}$ |  Sphere surface circle radius |  0.15m
$\Theta_{k}$ |  k-th ($k\in[1,..,92]$) camera view point Euler angles that depend on $r_{sph}$ and $r_{circ}$ |  92 - 2 excluded = 90 values. Concrete Euler angles can be derived from applying the sphere concept described in section 2.2.1
$z_{ego}$ |  Camera distance (egocentric
head to hand distance) |  $[0.5,0.8]$
$\mathcal{G}_{bkgr}$ |  Background image for rendered
grasp frames |  1 x plain white, 1 x consultation room, 3 x SPACE-FAN phantom, 3 x real pregnant mother belly (white / brown / black)
$H_{rgb}$ |  RGB values of glove texture |  $\left[\begin{array}[]{c}0.5647058824\\\ 0.5921568627\\\ 0.768627451\end{array}\right]$, $\left[\begin{array}[]{c}0.38039215686\\\ 0.61960784314\\\ 0.8666666667\end{array}\right]$
$A_{rgb}$ |  RGB values of the arm |  $[1.0,0.6784313725,0.3764705882]$

Using the variable notations of Table 1, the frame generation in our dataset
can be expressed as in Algorithm 1 below:

Algorithm 1 Combinatorics for frame generation based on a selected grasp

1:Input: Selected grasp $(G_{[\Theta_{Vol}]})_{k}$, egocentric viewpoint
$z_{ego}$, background image $\mathcal{G}_{bkgr}$, hand texture $H_{rgb}$,
camera view Euler angles $\Theta_{k}$, sphere radius $r_{sph}$, circle radius
$r_{circ}$

2:Output: Set of frames $F$

3:Initialize $F\leftarrow\emptyset$

4:for all
$(G_{[\Theta_{Vol}]})_{k},z_{ego},\mathcal{G}_{bkgr},H_{rgb},\Theta_{k}$ do

5: $frame\leftarrow
GenerateFrame((G_{[\Theta_{Vol}]})_{k},z_{ego},\mathcal{G}_{bkgr},H_{rgb},\Theta_{k},r_{sph},r_{circ})$

6: $F\leftarrow F\cup\\{frame\\}$

7:end for


# Analysis and application of multiplicative stochastic process with a sample-
dependent lower bound

Ken Yamamoto  Department of Physics and Earth Sciences, Faculty of Science,
University of the Ryukyus, Senbaru, Nishihara, Okinawa 903-0213, Japan
Yoshihiro Yamazaki  Department of Physics, Faculty of Advanced Science and
Engineering, Waseda University, Okubo, Shinjuku, Tokyo 169-8555, Japan

###### Abstract

A multiplicative stochastic process with the lower bound lognormally
distributed is investigated. For the process, the model is constructed, and
its distribution function (involving four parameters) and the related
statistical properties are derived. By adjusting the parameters, it is
confirmed that the theoretical distribution is consistent with empirical
distributions of some real data.

##  I Introduction

Probability distributions provide useful information on several complex
phenomena and help achieve effective modeling. Even if the precise dynamics or
underlying mechanisms are unknown, a probability distribution can be simply
constructed by observed data. As such, the analysis of the empirical
probability distributions is regarded as a possible first step toward
understanding complex systems.

The heavy-tailed probability distribution, which decays slower than
exponential distributions, is an important class of probability distributions
[1]. Power-law distribution [2] is representative of the heavy-tailed
distributions. The appearance of a heavy-tailed distribution implies that
there exist elements having a far larger size than the characteristic size and
that a mechanism for generating large-sized elements is embedded in the
system.

Meanwhile, lognormal distribution is another representative heavy-tailed
distribution that appears in various complex systems [3, 4]. Here, the random
variable $X$ is said to follow a lognormal distribution if its logarithm $\ln
X$ is normally distributed. The probability density $f(x)$ and cumulative
distribution $F(x)$ of the lognormal distribution can be described as follows:

$\displaystyle f(x)$ $\displaystyle=\frac{1}{\sqrt{2\pi}\sigma
x}\exp\left(-\frac{(\ln x-\mu)^{2}}{2\sigma^{2}}\right),$ (1a) $\displaystyle
F(x)$
$\displaystyle=\int_{x}^{\infty}f(x)dx=\frac{1}{2}\left[1-\operatorname{erf}\left(\frac{\ln
x-\mu}{\sqrt{2}\sigma}\right)\right],$ (1b)

where $\mu$ and $\sigma$ correspond to the mean and standard deviation of $\ln
X$, respectively, and $\operatorname{erf}$ is the Gauss error function defined
by the following:

$\operatorname{erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-y^{2}}dy.$

The lognormal distribution is generated by the following multiplicative
stochastic process

$X_{n+1}=B_{n}X_{n},$ (2)

where $B_{0},B_{1},\ldots$ represent positive random variables that are
independently and identically distributed. In fact, according to the central
limit theorem, the distribution of $\ln X_{n}$ converges to a normal
distribution. Qualitatively, because the lognormal distribution is attained
via numerous factors that are accumulated multiplicatively, the lognormal
distribution is observed in a variety of situations. This model was first
introduced by Kolmogorov [5] to explain a fragment-size distribution in a
cascade fracture.

The multiplicative stochastic process (2) can produce a power-law distribution
when an additional effect is introduced, one example being the introduction of
a lower bound $a$ for $X_{n}$ [6] such that $X_{n}$ is restricted to not less
than $a$:

$X_{n+1}=\begin{cases}B_{n}X_{n}&(B_{n}X_{n}\geq a)\\\
X_{n}&(B_{n}X_{n}<a).\end{cases}$ (3)

In the $n\to\infty$ limit, $X_{n}$ follows a stationary power-law
distribution:

$\Phi_{a}(x)=P(X_{n}>x)=\left(\frac{x}{a}\right)^{-\beta}\quad(x\geq a),$ (4)

where the exponent $\beta>0$ is given by the positive solution of

$E[B_{t}^{\beta}]=1.$

The corresponding probability density function is expressed as

$\phi_{a}(x)=-\frac{d\Phi_{a}(x)}{dx}=\frac{\beta}{a}\left(\frac{x}{a}\right)^{-\beta-1}\quad(x\geq
a).$ (5)

Historically, this model dates back to the work of Champernowne [7], who
applied it to the power-law income distribution. A continuous version that
involves geometric Brownian motion with a reflection boundary has been
investigated [8]. Other mechanisms for the power-law distributions from the
multiplicative process (2) have been also studied in terms of the introduction
of additive noise [9], a reset event [10], random stopping [11], and the sum
of variables [12].

In this paper, we investigate a stochastic process based on business-related
empirical data. Here, the empirical distribution of the data follows a
lognormal distribution up to the mid-sized region, with a tail fatter than the
usual lognormal distribution. We propose that this tail behavior can be
modeled via a multiplicative stochastic process with a specific lower bound,
in which the initial value is drawn from a lognormal distribution and the
lower bound is equal to this initial value. After formulating the model, we
derive the distribution function and its properties.

##  II Modeling

Figure 1 shows the cumulative number of troubleshooting costs in a specific
Japanese company. The cumulative number is defined as the number of cases with
costs larger than $x$. The absolute values are substituted with the costs
divided by the average as shown in Fig. 1 [13]. The total number of cases is
$N=2571$. The solid curve corresponds to the lognormal distribution with
$\mu=-0.84$ and $\sigma=1.24$. (Note that the cumulative number is given by
$NF(x)$.) The empirical distribution deviates from the lognormal distribution
in the tail ($x\gtrsim 10$).

Figure 1:  The cumulative number of troubleshooting costs (normalized by the
average cost) in the company (the number of cases is $N=2571$). The solid
curve is the lognormal distribution with $\mu=-0.84$ and $\sigma=1.24$, and
the dashed line indicates the power-law slope $x^{-1.53}$ in the tail.

While only a few (the order of $10^{1}$) points deviate from the lognormal
distribution, they are highly expensive cases, ranging from approximately $10$
to $100$ times greater than the average value. Hence, the precise description
of the tail behavior based on a mathematical mechanism will be financially
important in practical situations. Moreover, it has been widely observed that
the main part of the data follows a lognormal distribution, but the tail
deviates from the lognormal [14, 15, 16]. Such tail behavior is generally
difficult to determine and can often be controversial. For example, while the
population distribution of American cities was reported to be lognormally
distributed [17], its tail was later noted to exhibit a power-law decay [18].
This issue has not been resolved yet because the measurement of the tail is
problematic [19, 20]. Furthermore, the population distribution varies
according to city, town, and village [21], and a numerical model incorporating
the natural and social effects on population increase or decrease in a
simplified form has been proposed [21, 22]. Therefore, the tail of an
empirical distribution is worth modeling and analyzing in theoretical terms.

Assuming the tail of the empirical distribution shown in Fig. 1 follows power-
law decay, we obtain the power-law exponent $-1.53\pm 0.08$ using the method
by Clauset, Shalizi, and Newman [23]. However, this power-law decay (shown by
the dashed line in Fig. 1) does not seem to express the tail of the empirical
distribution well. The power-law decay cannot always be justified
theoretically, and it should be regarded as a rough assumption or a working
hypothesis. A more appropriate function may be required to understand the
system accurately.

We thus propose the following phenomenological model. Based on the lognormal
distribution in the low- and mid-cost regions, it is assumed that the cost is
determined by several multiplicatively affecting factors. Thus, this can be
naively modeled using the stochastic process described in Eq. (2), which
produces the lognormal distribution. Additionally, for expensive cases,
additional tasks, such as individual inspection and customization, tend to be
involved. These additional tasks may increase the cost but never reduce it to
less than the initial cost. This operation can be expressed by a
multiplicative process wherein the initial cost acts as the lower bound. Due
to this additional process, the price can further increase, resulting in a
deviation from the lognormal distribution.

The above model can be simulated as follows. First, let us define $x_{*}$ as
the boundary cost where the empirical data begins to deviate from the
lognormal distribution. A random number $x$ is then drawn from the lognormal
distribution. If $x\geq x_{*}$, the following multiplicative process, wherein
the lower bound is $x$ itself, is performed, where the initial value is
$X_{0}=x$ and the evolution is given by

$X_{n+1}=\begin{cases}B_{n}X_{n}&(B_{n}X_{n}\geq x)\\\
X_{n}&(B_{n}X_{n}<x)\end{cases}$

until $X_{n}$ becomes stationary. Here, $B_{0},B_{1},\ldots$ are independently
and identically distributed positive random numbers. Unlike in Eq. (3), where
the lower bound $a$ is constant, the lower bound $x$ is distributed
lognormally and varies for each sample in this model.

Therefore, $x$ is lognormally distributed for $x<x_{*}$, whereas the
distribution changes for $x\geq x_{*}$ due to the additional process. In the
next section, we analyze the proposed model and derive the distribution
function for $x\geq x_{*}$.

##  III Model analysis

As mentioned earlier, the proposed model involves a multiplicative process
with a varying lower bound. First, we derive the corresponding probability
density $g(x)$ and cumulative distribution $G(x)$.

For $x<x_{*}$, the effect of the lower bound does not come into play. Hence,
$g(x)$ and $G(x)$ are identical to the lognormal distribution $f(x)$ and
$F(x)$ described in Eq. (1), respectively.

Figure 2:  Illustration of the interchange of the integration in Eq. (7). The
shaded area composed of semi-infinite triangular (corresponding to $G_{1}(x)$)
and rectangular ($G_{2}(x)$) areas represent the domain of the iterated
integral in the $y$-$a$ plane.

For $x\geq x_{*}$, $g(x)$ is obtained by superposing the stationary
distributions corresponding to the different values of the lower bound. When
the lower bound is $a$, the stationary distribution follows the power law
$\phi_{a}(x)$ described in Eq. (5). In the proposed model, the lower bound $a$
has a lognormal weight $f(a)$ [see Eq. (1a)], such that

$g(x)=\int_{x_{*}}^{x}\phi_{a}(x)f(a)da.$ (6)

The cumulative distribution can be expressed as follows:

$G(x)=\int_{x}^{\infty}g(y)dy=\int_{x}^{\infty}dy\int_{x_{*}}^{y}\phi_{a}(y)f(a)da.$

The order of integration can be interchanged as follows (see Fig. 2 for
reference):

$\int_{x}^{\infty}dy\int_{x_{*}}^{y}da\phi_{a}(y)f(a)=\left(\int_{x}^{\infty}da\int_{a}^{\infty}dy+\int_{x_{*}}^{x}da\int_{x}^{\infty}dy\right)\phi_{a}(y)f(a)=G_{1}(x)+G_{2}(x).$
(7)

By the normalization $\int_{a}^{\infty}\phi_{a}(y)dy=1$, we can easily obtain

$G_{1}(x)=\int_{x}^{\infty}f(a)da\int_{a}^{\infty}\phi_{a}(y)dy=\int_{x}^{\infty}f(a)da=F(x).$
(8)

The function $G_{2}(x)$ can be written as

$\displaystyle G_{2}(x)$
$\displaystyle=\int_{x_{*}}^{x}daf(a)\int_{x}^{\infty}\phi_{a}(y)dy=\int_{x_{*}}^{x}f(a)\Phi_{a}(x)da$
$\displaystyle=\int_{x_{*}}^{x}\frac{1}{\sqrt{2\pi}\sigma
a}\exp\left(-\frac{(\ln
a-\mu)^{2}}{2\sigma^{2}}\right)\left(\frac{x}{a}\right)^{-\beta}da$
$\displaystyle=x^{-\beta}\exp\left(\mu\beta+\frac{\sigma^{2}\beta^{2}}{2}\right)\int_{x_{*}}^{x}\frac{1}{\sqrt{2\pi}\sigma
a}\exp\left(-\frac{(\ln a-\mu-\beta\sigma^{2})^{2}}{2\sigma^{2}}\right)da.$
(9)

Since the integrand represents the lognormal probability density (as a
function of $a$) where the parameter $\mu$ is replaced with
$\mu+\beta\sigma^{2}$, we have

$G_{2}(x)=\frac{x^{-\beta}}{2}\exp\left(\mu\beta+\frac{\sigma^{2}\beta^{2}}{2}\right)\left[\operatorname{erf}\left(\frac{\ln
x-\mu-\beta\sigma^{2}}{\sqrt{2}\sigma}\right)-\operatorname{erf}\left(\frac{\ln
x_{*}-\mu-\beta\sigma^{2}}{\sqrt{2}\sigma}\right)\right].$ (10)

Thus, for $x\geq x_{*}$, we can obtain

$G(x)=\frac{1}{2}\left[1-\operatorname{erf}\left(\frac{\ln
x-\mu}{\sqrt{2}\sigma}\right)\right]+\frac{1}{2}\exp\left(\mu\beta+\frac{\sigma^{2}\beta^{2}}{2}\right)x^{-\beta}\left[\operatorname{erf}\left(\frac{\ln
x-\mu-\beta\sigma^{2}}{\sqrt{2}\sigma}\right)-\operatorname{erf}\left(\frac{\ln
x_{*}-\mu-\beta\sigma^{2}}{\sqrt{2}\sigma}\right)\right].$ (11)

Here, the first term on the right-hand side represents the lognormal
distribution, while the second term is regarded as the attendant correction
term.

By differentiating Eq. (11), we obtain

$\displaystyle g(x)$ $\displaystyle=-\frac{dG(x)}{dx}$
$\displaystyle=\frac{1}{2}\beta
x^{-\beta-1}\exp\left(\mu\beta+\frac{\sigma^{2}\beta^{2}}{2}\right)\left[\operatorname{erf}\left(\frac{\ln
x-\mu-\beta\sigma^{2}}{\sqrt{2}\sigma}\right)-\operatorname{erf}\left(\frac{\ln
x_{*}-\mu-\beta\sigma^{2}}{\sqrt{2}\sigma}\right)\right]$

for $x\geq x_{*}$.

Figure 3 shows graphs of $g(x)$ and $G(x)$ with fixed $(\mu,\sigma)=(0,1.5)$
and $x_{*}=10$, and varying $\beta=0.5,1,2,4$, and $8$. The tail decays
rapidly for large $\beta$, and the probability density $g(x)$ is not
continuous at $x=x_{*}$: $g(x)\to 0$ as $x\searrow x_{*}$. In accordance with
the discontinuity of $g(x)$, the $G(x)$ graph is not smooth at $x=x_{*}$.

(a) (b)

Figure 3:  Graphs of (a) the probability density $g(x)$ and (b) the cumulative
distribution $G(x)$ for $(\mu,\sigma)=(0,1.5)$, $x_{*}=10$, and
$\beta=0.5,1,2,4$, and $8$. The decay is fast for large $\beta$ (indicated by
the arrow in each graph).

If a random variable $X$ has the probability density $g(x)$, the $k$th moment
of $X$ can be calculated as

$\displaystyle E[X^{k}]=$
$\displaystyle\int_{0}^{x_{*}}x^{k}f(x)dx+\int_{x_{*}}^{\infty}x^{k}g(x)dx$
$\displaystyle=\frac{1}{2}\exp\left(k\mu+\frac{k^{2}\sigma^{2}}{2}\right)\left[\frac{2\beta-k}{\beta-k}-\frac{k}{\beta-k}\operatorname{erf}\left(\frac{\ln
x_{*}-\mu-k\sigma^{2}}{\sqrt{2}\sigma}\right)\right]$

for $k<\beta$, and the moment diverges for $k\geq\beta$. The calculation of
the integral is straightforward but cumbersome. This criterion for the
divergence of moments can be derived simply using the asymptotic power law
$g(x)\sim x^{-\beta-1}$ or $G(x)\sim x^{-\beta}$.

##  IV Data analysis

The function $G(x)$ involves four adjustable parameters: $\mu,\sigma,x_{*}$,
and $\beta$. The optimal value of these parameters for fitting the given data
can be determined as follows. For fixed $x_{*}$, we can obtain $\mu$ and
$\sigma$ according to the lognormal fitting in the region $x<x_{*}$. Using
these values $\mu$ and $\sigma$, we can obtain $\beta$ by fitting Eq. (11) in
the region $x>x_{*}$. By calculating $\mu,\sigma,\beta$, and the fitting error
for different $x_{*}$, we can determine the optimal set of the parameters as
one minimizing the fitting error.

Figure 4 shows the result of the curve fitting of $G(x)$ to the empirical
distribution shown in Fig. 1. Here, the parameter values are
$(\mu,\sigma)=(-0.84,1.24)$ (same as those in Fig. 1), $x_{*}=4.21$, and
$\beta=5.28$. $G(x)$ is consistent with the empirical distribution almost over
the entire range, which indicates the validity of the proposed model i.e., the
multiplicative stochastic process with the sample-dependent lower bound.

The naive power-law fitting shown in Fig. 1 gives the power-law exponent
$-1.53(>-2)$ for a cumulative distribution, implying that the variance
diverges if the same power law extends to infinity. In contrast, we obtained
$\beta=5.28$ when using $G(x)$, which results in finite variance. Here, it can
be expected to obtain the correct statistical properties of the data when
using $G(x)$, ensuring appropriate and effective cost management in practical
situations.

Figure 4:  The cumulative number of normalized trouble handling costs (circle
points) and $G(x)$ (solid curve). The fitting parameters are
$(\mu,\sigma)=(-0.84,1.24)$, $x_{*}=4.21$, and $\beta=5.28$. The vertical
dashed line indicates the position of $x_{*}$.

As outlined in Sec. III [see Eqs. (8) and (10)], $G(x)$ for $x\geq x_{*}$ can
be decomposed into two parts $G(x)=G_{1}(x)+G_{2}(x)$, where $G_{1}(x)$ is a
lognormal function, and $G_{2}(x)$ is a power-law-like function (not an exact
power law due to the $\operatorname{erf}$ factor). Figure 5(a) shows
$G_{1}(x)$ and $G_{2}(x)$ for $(\mu,\sigma)=(-0.84,1.24)$, $x_{*}=4.21$, and
$\beta=5.28$ (same values as in Fig. 4). In this case, the magnitude of
$G_{2}(x)$ is greater than that of $G_{1}(x)$ for $x>25.8$, and the power-law-
like part $G_{2}$ is dominant in the tail. Specifically, as Fig. 5(b) shows,
the ratio $G_{2}(x)/G(x)$ represents an increasing function. Furthermore,
using the asymptotic expansion
$1-\operatorname{erf}(z)\simeq\exp(-z^{2})/(z\sqrt{\pi})$, $G_{2}(x)/G(x)\to
1$ can be presented as $x\to\infty$ for any parameter value.

(a) (b)

Figure 5:  (a) Decomposition of $G(x)$ into the lognormal part $G_{1}(x)$ and
power-law-like part $G_{2}(x)$. The parameter values are the same as those in
Fig. 4: $(\mu,\sigma)=(-0.84,1.24)$, $x_{*}=4.21$, and $\beta=5.28$. (b) The
ratio $G_{2}(x)/G(x)$ increases with increasing $x$, which suggests that the
power-law-like part $G_{2}(x)$ is dominant in the tail.

In addition to the troubleshooting costs, we present the empirical
distributions of other social data. Figure 6(a) shows the distribution of the
number of employees aggregated according to a municipality ($N=1737$) in Japan
[24] in the form of the cumulative number. Here, the solid curve represents
$G(x)$ for $(\mu,\sigma)=(9.17,1.59)$, $x_{*}=3.09\times 10^{5}$, and
$\beta=4.20$. We believe that the lognormal behavior in the $x<x_{*}$ region
is related to the lognormal population distribution [21], and that the heavy
tail likely reflects the extreme concentration of companies in metropolises,
such as Tokyo and Osaka. Meanwhile, Fig. 6(b) shows the cumulative number
(points) of the management expenses grants for 90 national universities in
Japan in 2019 [25] and $G(x)$ for $(\mu,\sigma)=(2.11,0.82)$,
$x_{*}=$21\text{\times}{10}^{9}\text{\,}\mathrm{y}\mathrm{e}\mathrm{n}$$, and
$\beta=3.80$ (solid curve). Here, $12$ universities have $x$ greater than
$x_{*}$. While it is not clear whether these data are generated by the same
statistical mechanism as that outlined in this study in terms of the
troubleshooting costs, we believe that $G(x)$ effectively expresses the
empirical distribution; thus, providing some insight into these phenomena.

(a) (b)

Figure 6:  (a) Cumulative number of employees in municipalities ($N=1737$) in
Japan in 2014. The solid curve represents $G(x)$ for
$(\mu,\sigma)=(9.17,1.59)$, $x_{*}=3.09\times 10^{5}$, and $\beta=4.20$. (b)
Cumulative number of management expenses grants for the national universities
($N=90$) in Japan in 2019. The solid curve corresponds to $G(x)$ for
$(\mu,\sigma)=(2.11,0.82)$,
$x_{*}=$21\text{\times}{10}^{9}\text{\,}\mathrm{y}\mathrm{e}\mathrm{n}$$, and
$\beta=3.80$.

Next, we compare the proposed distribution $G(x)$ with other distribution
functions. We selected the following three distributions that have been
frequently used in the analysis of empirical data: the lognormal distribution
$F(x)$ [shown in Eq. (1b)], the lognormal distribution with power-law tail,
and the double lognormal distribution (a mixture of two lognormal
distributions) given by

$\frac{\alpha}{2}\left[1-\operatorname{erf}\left(\frac{\ln
x-\mu_{1}}{\sqrt{2}\sigma_{1}}\right)\right]+\frac{1-\alpha}{2}\left[1-\operatorname{erf}\left(\frac{\ln
x-\mu_{2}}{\sqrt{2}\sigma_{2}}\right)\right],$

where $\mu_{1},\sigma_{1},\mu_{2},\sigma_{2}$, and $\alpha$ are the fitting
parameters. Figure 7 shows the cumulative number of employees in
municipalities (a) and management expenses grants for the Japanese national
universities (b) (the datasets are the same as with Fig. 6), and these three
distributions.

(a)
(b)

Figure 7:  Results of the curve fitting for (a) the number of employees and
(b) the management expenses grants for the Japanese national universities by
the lognormal distribution (left), the lognormal distribution and power-law
tail (center), and the double lognormal distribution (right). The vertical
dashed line in the center graphs indicates the border between the lognormal
and power-law distributions.

To obtain the center graphs of Fig. 7, we performed the following procedure.
First, the power-law exponent and the lower bound of the power-law
distribution were determined using the method outlined in Ref. [23], and we
regarded this lower bound as the crossover point of the power-law and
lognormal distributions. Next, the lognormal fitting was performed in the
range up to the crossover point. Finally, the coefficient of the power-law
distribution was computed so that the power-law and lognormal functions were
continuously connected.

To evaluate the quality of curve fitting using the four distributions (the
lognormal distribution, the lognormal distribution with power-law tail, the
double lognormal distribution, and the proposed distribution), we measured the
distance between the empirical distribution and each candidate distribution.
Among a wide variety of distances or metrics between probability
distributions, we selected two simple ones: the Kolmogorov-Smirnov statistic
$D$ and the first Wasserstein distance $W_{1}$. If the cumulative
distributions of the empirical data and the fitting result are respectively
written as $F_{\mathrm{emp}}$ and $F_{\mathrm{fit}}$, these two distances are
defined as

$\displaystyle D(F_{\mathrm{fit}},F_{\mathrm{emp}})$
$\displaystyle=\sup_{0<x<\infty}|F_{\mathrm{fit}}(x)-F_{\mathrm{emp}}(x)|,$
$\displaystyle W_{1}(F_{\mathrm{fit}},F_{\mathrm{emp}})$
$\displaystyle=\int_{0}^{\infty}|F_{\mathrm{fit}}(x)-F_{\mathrm{emp}}(x)|dx.$

These distances take non-negative values, and the small distances indicate
that the estimated distribution is close to the empirical one.

Table 1:  Distances between the empirical data and the fitting distributions.
The datasets are (a) number of employees, (b) management expenses grants, and
(c) troubleshooting costs. The relative values to the case of the proposed
distribution $G$ are shown.

(a) Number of employees
Lognormal Lognormal Double +power lognormal Relative $D$ 1.49 1.00 1.00
Relative $W_{1}$ 1.78 1.78 0.95
(b) Management expenses grants
Lognormal Lognormal Double +power lognormal Relative $D$ 1.21 1.65 1.12
Relative $W_{1}$ 1.87 2.62 1.28
(c) Troubleshooting costs
Lognormal Lognormal Double +power lognormal Relative $D$ 1.49 1.00 1.00
Relative $W_{1}$ 1.78 1.78 0.95

The calculation result of $W_{1}$ and $D$ for the lognormal distribution,
lognormal distribution with power-law tail, and double lognormal distribution
in Fig. 7, in comparison with the proposed distribution $G$ in Figs. 6 and 4,
is presented in Table 1. The relative values (the distance between each
distribution and the empirical distribution divided by the case of the
proposed distribution $G$) are shown in the table; the value greater than
unity represents that the distribution $G$ is closer to the empirical
distribution than the corresponding distribution with respect to the
corresponding distance.

In the data analysis of complex systems including social ones, the tail of
empirical distributions has often been focused on. However, the above
statistical distances are not suitable for the evaluation of the tail
behavior. In fact, Fig. 8 shows the absolute difference
$|F_{\mathrm{fit}}(x)-F_{\mathrm{emp}}(x)|$ of the empirical distribution for
the number of employees and its double-lognormal fitting result [Fig. 7(a)
right] as a function of $x$, indicating that large differences occur in the
region $x\lesssim 10^{5}$. (A similar tendency is found in the other two
datasets and other distributions.) Hence, $W_{1}$ and $D$ are dominantly
affected by the body part, and the contribution of differences in the tail is
relatively small. The four distributions (lognormal, lognormal with power-law
tail, double lognormal, and $G$) have lognormal forms for small $x$, and
similar $W_{1}$ and $D$ values in Table 1 are considered to be attributed to
this common lognormality.

Figure 8:  The absolute difference $|F_{\mathrm{fit}}(x)-F_{\mathrm{emp}}(x)|$
for the case of the distribution of the number of employees and the double
lognormal distribution [corresponding to the right panel of Fig. 7(a)]. The
differences in the tail ($x\gtrsim 10^{5}$) are small compared with those in
the body ($x\lesssim 10^{5}$).

To magnify the difference in the tail, we propose two distances:

$\Delta_{1}(F_{\mathrm{fit}},F_{\mathrm{emp}})=\int_{x_{\mathrm{min}}}^{x_{\mathrm{max}}}|\ln
F_{\mathrm{fit}}(x)-\ln F_{\mathrm{emp}}(x)|dx$

and

$\Delta_{2}(F_{\mathrm{fit}},F_{\mathrm{emp}})=\int_{x_{\mathrm{min}}}^{x_{\mathrm{max}}}\frac{|F_{\mathrm{fit}}(x)-F_{\mathrm{emp}}(x)|}{F_{\mathrm{emp}}(x)}dx,$

where $x_{\mathrm{min}}$ and $x_{\mathrm{max}}$ are the minimum and maximum
values of the empirical data, respectively. (The upper limit of the integral
is set to $x_{\mathrm{max}}$ in order to avoid the singularity
$F_{\mathrm{emp}}(x)=0$ for $x>x_{\mathrm{max}}$.) Like other statistical
distances, $\Delta_{1}$ and $\Delta_{2}$ are always nonnegative and they
become smaller as the estimated distribution $F_{\mathrm{fit}}$ is closer to
the empirical one $F_{\mathrm{emp}}$.

Table 2:  Distances $\Delta_{1}$ and $\Delta_{2}$ between the empirical data
and the fitting distributions. The datasets are (a) number of employees, (b)
management expenses grants, and (c) troubleshooting costs. The relative values
to the case of the proposed distribution $G$ are shown.

(a) Number of employees
Lognormal Lognormal Double +power lognormal Relative $\Delta_{1}$ 4.58 2.96
1.13 Relative $\Delta_{2}$ 3.88 2.56 1.09
(b) Management expenses grants
Lognormal Lognormal Double +power lognormal Relative $\Delta_{1}$ 1.42 2.38
1.26 Relative $\Delta_{2}$ 1.35 1.95 1.27
(c) Troubleshooting costs
Lognormal Lognormal Double +power lognormal Relative $\Delta_{1}$ 1.87 3.06
1.08 Relative $\Delta_{2}$ 1.71 4.32 1.12

The quality of curve fitting based on $\Delta_{1}$ and $\Delta_{2}$ is
presented in Table 2. As in Table 1, the relative values of the three
distributions in Fig. 7 against the distribution $G$ in Figs. 6 and 4 are
shown in this table. The proposed distribution $G$ describes the tail of
empirical distributions better than (or at least comparable to) the other
three distributions (lognormal, lognormal with power-law tail, and double
lognormal).

In addition to the quantitative aspect, another advantage of $G$ is that it
possesses a mechanism (stochastic process) different from the other three
distributions. That is, the agreement of the empirical dataset with the
proposed distribution implies that the data was generated under the effect of
the proposed mechanism, namely a multiplicative stochastic process with a
sample-dependent lower bound. This conclusion is based on only the three
specific datasets shown in Figs. 4 and 6, and general discussion on the
quality of each candidate distribution is the subject for a future study.

##  V Discussion

The empirical data of troubleshooting costs outlined in this paper is
normalized in terms of the average. Here, we analyze the effect of this
normalization. If we set $m$ as the average cost, the raw value $\tilde{x}$
and the normalized value $x$ satisfy $x=\tilde{x}/m$. The threshold $x_{*}$ is
transformed into $\tilde{x_{*}}=mx_{*}$, since $x\lessgtr x_{*}$ changes to
$\tilde{x}\lessgtr mx_{*}=\tilde{x_{*}}$, while the parameter $\mu$ in $G(x)$
appears in the form of $\ln x-\mu$ or $\ln x_{*}-\mu$ at all times [the factor
$\exp(\mu\beta)x^{-\beta}$ is reduced to $\exp(-\beta(\ln x-\mu))$].
Therefore, the parameter $\mu$ becomes $\tilde{\mu}=\mu+\ln m$. It should be
noted that the parameters $\sigma$ and $\beta$ in $G$ do not change due to the
normalization. In conclusion, writing the parameters explicitly as
$G(x;\mu,\sigma,x_{*},\beta)$, we have
$G(x;\mu,\sigma,x_{*},\beta)=G(\tilde{x};\mu+\ln m,\sigma,mx_{*},\beta)$. The
same dependence of the parameters holds even if the constant $m$ is not equal
to the average value.

The distributions $g(x)$ and $G(x)$ are derived based on the assumption that
the lower bound of the multiplicative process is lognormally distributed. We
can generalize these results by replacing the lognormal weight $f(a)$
described in Eq. (6) with alternative probability distributions. Here, when a
probability density $h$ and its cumulative distribution $H$ are used instead
of the lognormal distribution $f$ and $F$, we obtain $G_{1}(x)=H(x)$ and

$G_{2}(x)=x^{-\beta}\int_{x_{*}}^{x}a^{\beta}h(a)da,$

which are similar to Eqs. (8) and (9). The asymptotic power law $G(x)\simeq
G_{2}(x)\sim x^{-\beta}$ generally holds if $h$ decays faster than
$x^{-\beta-1}$.

The tail of an empirical probability distribution usually contains only a
small amount of data. For example, while the total number of cases is $N=2571$
in Fig. 4, only $85$ cases ($3.3\text{\,}\mathrm{\char 37\relax}$ of the
total) lie in $x\geq x_{*}$. Due to statistical fluctuations, the distribution
tail is difficult to specify. From Figs. 4 and 6, the proposed distribution
$G$ fits well even up to the tail, which appears to increase reliability for
the estimation of the asymptotic power-law exponent $\beta$. By its
construction, the proposed model, a multiplicative process with a sample-
dependent lower bound, will be useful for a lognormal distribution with a fat
tail, especially for a system whose dynamics or evolution rules change at a
certain threshold value. In fact, the troubleshooting cost, the number of
employees, and the university management expenses grants analyzed in this
study are real examples which follow the proposed distribution $G$, and this
result suggests that the proposed mechanism may be embedded in the generation
process of them. Moreover, we expect that there exist other systems that can
be well described by the proposed model.

We think that the proposed distribution $g$ or $G$ is an effective
distribution but not the definitive one. In general, the distribution function
of a given dataset should be determined by examining several candidate
distributions, and we hope that the distribution $G$ outlined in this paper
becomes one of the possible choices.

## Acknowledgments

The authors are very grateful to Mr. Hirokazu Yanagawa and Mr. Takashi Bando
for constructive discussion and permission to use data. KY was supported by a
Grant-in-Aid for Scientific Research (C) 19K03656 from Japan Society for the
Promotion of Science.

## References

  * [1] S. Foss, D. Korshunov, and S. Zachary, An Introduction to Heavy-Tailed and Subexponential Distributions, Springer Science & Business Media, New York, 2013.
  * [2] M. E. J. Newman, Contemporary Phys. 46, 323 (2005).
  * [3] N. Kobayashi, H. Kuninaka, J. Wakita, and M. Matsushita, J. Phys. Soc. Jpn. 80, 072001 (2011).
  * [4] E. Limpert, W. A. Stahel, and M. Abbt, BioScience 51, 341 (2001).
  * [5] A. N. Kolmogorov, Dokl. Akad. Nauk SSSR 31, 99 (1941).
  * [6] M. Levy and S. Solomon, Int. J. Mod. Phys. C 7, 595 (1996).
  * [7] D. G. Champernowne, Econ. J. 63, 318 (1953).
  * [8] X. Gabaix, Quarterly J. Econ. 114, 739 (1999).
  * [9] H. Takayasu, A. Sato, and M. Takayasu, Phys. Rev. Lett. 79, 966 (1997).
  * [10] S. C. Manrubia and D. H. Zanette, Phys. Rev. E 59, 4945 (1999).
  * [11] K. Yamamoto and Y. Yamazaki, Phys. Rev. E 85, 011145 (2012).
  * [12] K. Yamamoto, Phys. Rev. E 89, 042115 (2014).
  * [13] At the company’s request, its name cannot be revealed, and the use of the absolute cost values was prohibited.
  * [14] F. Clementi and M. Gallegati, Physica A 350, 427 (2005).
  * [15] M. Golosovsky and S. Solomon, EPJ Special Topics 205, 303 (2012).
  * [16] Á. Corral and Á. González, Earth and Space Sci. 6, 673 (2019).
  * [17] J. Eeckhout, Amer. Econ. Rev. 94, 1429 (2004).
  * [18] M. Levy, Amer. Econ. Rev. 99, 1672 (2009).
  * [19] J. Eeckhout, Amer. Econ. Rev. 99, 1676 (2009).
  * [20] Y. Ioannides and S. Skouras, SSRN 1481254 (2009). http://dx.doi.org/10.2139/ssrn.1481254
  * [21] Y. Sasaki, H. Kuninaka, and M. Matsushita, J. Phys. Soc. Jpn. 76, 074801 (2007).
  * [22] Y. Yamazaki and K. Takamura, J. Phys. Soc. Jpn. 82, 065003 (2013).
  * [23] A. Clauset, C. R. Shalizi, and M. E. J. Newman, SIAM Rev. 51, 661 (2009).
  * [24] Japan National Statistics Center, Standardized Statistical Data Set for Education https://www.nstac.go.jp/SSDSE/pastdata.html [in Japanese]
  * [25] Japan Ministry of Education, Culture Sports, Science and Technology. (2020, November 17). Management Expenses Grants for National Universities in Japan in 2019. https://www.mext.go.jp/kaigisiryo/content/20201117-mxt_hojinka-000011109_3.pdf [in Japanese]