# Comparing Apples to Oranges: LLM-powered Multimodal Intention Prediction in
an Object Categorization Task

Hassan Ali, Philipp Allgeuer and Stefan Wermter1 1The research was supported
by Horizon Europe project TERAIS, under Grant agreement number 101079338 and
the DFG under the Crossmodal Learning (TRR-169) project. The authors are with
the Knowledge Technology Group, Department of Informatics, University of
Hamburg, Germany. Emails: {hassan.ali, philipp.allgeuer, stefan.wermter}@uni-
hamburg.de

###### Abstract

Intention-based Human-Robot Interaction (HRI) systems allow robots to perceive
and interpret user actions to proactively interact with humans and adapt to
their behavior. Therefore, intention prediction is pivotal in creating a
natural interactive collaboration between humans and robots. In this paper, we
examine the use of Large Language Models (LLMs) for inferring human intention
during a collaborative object categorization task with a physical robot. We
introduce a hierarchical approach for interpreting user non-verbal cues, like
hand gestures, body poses, and facial expressions and combining them with
environment states and user verbal cues captured using an existing Automatic
Speech Recognition (ASR) system. Our evaluation demonstrates the potential of
LLMs to interpret non-verbal cues and to combine them with their context-
understanding capabilities and real-world knowledge to support intention
prediction during human-robot interaction.

##  I Introduction

With the rapid advancements in AI, humanoid robots are set to play
increasingly prominent roles in human societies, spanning diverse sectors like
manufacturing, healthcare, and domestic environments. This technological shift
is expected to foster a natural and intuitive interaction between humans and
their robotic collaborators, hence resembling human-human interaction rather
than just relying on explicit commands [1]. Intention prediction is
fundamental in HRI and involves perceiving subtle cues from human actions,
gestures, and speech to infer underlying goals [2]; therefore, intention
prediction is key for enabling effective and seamless communication.
Nevertheless, interpreting and predicting human actions is challenging due to
the multimodal nature of human behavior, being complex and driven by
situational context. Also, human intentions are often conveyed with nuanced
non-verbal cues which are challenging to interpret.

Typically, intention-based systems require collecting a vast amount of data,
which is both impractical and restrictive, limiting the systems’ applicability
to specific use cases or a particular set of modalities and/or setups [3].
Recent research with LLMs in robotics have showcased their abstract reasoning
capabilities going beyond mere language applications. Besides their
proficiency in text comprehension and generation, LLMs have exhibited signs of
notable inferential skills, extensive real-world knowledge and robust
reasoning capabilities [4]. Additionally, they are performant at capturing
context, making them a candidate for advancing intention prediction in HRI.
However, LLMs heavily depend on textual input (user queries and prompts) to
generate contextually relevant responses. Thus, when used with robots, it
remains challenging for LLMs to replace human-aware models during the
interaction [5], even with the subsequent generation of multimodal LLMs, which
still face difficulties processing and understanding visual information and
basic visual patterns [6]. Thus, a reliable approach for equipping LLMs with
awareness of human non-verbal cues is still missing.

In this work, we provide the Neuro-Inspired COLlaborator [7] (NICOL) robot
(see Fig. 1) with the capacity to predict human intentions and anticipate
actions by integrating several verbal and non-verbal cues with an LLM-powered
intent reasoning framework. By following a two-layered processing approach,
referred to as perceptive reasoning and task reasoning, we process user non-
verbal cues, captured with our machine learning vision models. Then, we use a
text-based LLM’s inferential skills to identify current user actions relevant
for the collaborative environment and anticipate user goals behind the
inferred actions. We evaluate our system through a collaborative task for
object categorization. Our experiments confirm the potential of LLMs as viable
assets for building a competent intention prediction system. Also, we showcase
the importance of non-verbal communication for an intuitive and natural
collaboration with the robot. Overall, our work contributes to both research
in LLMs and robotics as well as human intention prediction in human-robot
collaboration, effectively bridging the gap between machine understanding and
the subtleties of human communication.

Figure 1: An overview of our intention prediction system. After inferring both
verbal and non-verbal cues using an LLM, the robot generates appropriate
actions, e.g., giving the bowl to the user.

##  II Related Work

Predicting user intention in HRI is a considerable challenge, given the
immense individual human differences and the complexity of dynamic
environments [3]. Also, efficient prediction of human intention requires
intensive comprehension of context and continuous awareness of the user [8].
Besides LLMs’ proficiency in language comprehension, they have shown promising
capabilities in contextual understanding [9]. They also exhibit adaptability
to user input and requirements, while possessing potential for processing
real-world knowledge [10]. Thus, LLMs are potential assets for tasks including
real-time coordination between a human and an AI agent [11], trust-based
human-aware zero-shot planning [5], belief-desire-intention driven planning
[12], and human-aware planning and action anticipation [13]. LLMs have also
been useful in anticipating human actions while tuned in a few-shot learning
fashion [14], and employed as tools for solving ambiguity in human
instructions through multiple dialogue rounds for clarification of intentions
[15].

While explaining intentions with verbal cues and conversation repair is
beneficial, incorporating non-verbal communication attributes to a seamless
and intuitive interaction, thereby increasing human awareness in AI systems.
Prior conducted work primarily utilized non-verbal cues for providing robots
with human-like capabilities for apparent enhancement of the HRI experience.
Some research included grounding LLMs in humanoid robots to generate movements
and poses [16] and improving the embodiment of a conversational agent with a
focus on engagement through features like facial expressions [17]. Other
approaches focused on improving the interaction liveliness by hand gestures
[18]. Work was also done on using the robot eye gaze for managing dialogue
while interacting with multiple users [19]. However, there remains a notable
gap in leveraging human non-verbal cues for the exact purpose of improving
intention prediction. In our work, we demonstrate that integrating non-verbal
cues with LLMs contributes to accurate predictions of human intentions,
ultimately enriching the interaction quality.

##  III Methodology

Prompt 1: Our system is based on a grounded LLM as the NICOL robot [20].
Input:
USER: Tell me about yourself.
Output:
SAY: Hello! I am NICOL, short for ”Neuro-Inspired Collaborator”. I am a
humanoid robot designed and built by the Knowledge Technology group of the
University of Hamburg, and I am here to collaborate with you!
0.4pt3pt 2pt Input:
USER: Give me the red object which is also a fruit.
Output:
SAY: I see that you’re referring to the apple. Sure I can give it to you.
ACTION: Give apple

We base our system on a robustly grounded ChatGPT LLM [20], mirroring the
behavior and interaction style of a physical robot, i.e., NICOL robot, rather
than interacting with users as an AI language model. This is achieved by
equipping the LLM with facts about the environment through the robot sensors,
and linking the LLM abstract reasoning with the robot physical actions in a
modular open architecture. An example response by the said LLM is in Prompt
III. The framework processes user queries and generates responses in first
person, encoding robot actions without explicit programming. These LLM-
generated responses are parsed into tokens recognized by NICOL’s ROS-based API
platform, then mapped into concrete actions (ex: SAY to generate speech or
ACTION to generate a robotic manipulation).

Our framework (see Fig. 2) utilizes the described concept [20] to implement
and operate through a dual process for intention prediction: perceptive
reasoning and task reasoning. Perceptive reasoning refers to the robot’s
ability to collect and analyze information about the surrounding environment
and interacting user, thus formulating a set of “beliefs” guiding its
behavior. In contrast, task reasoning involves incorporating user guidance,
explicit commands, and contextual factors to collaborate towards a mutual
goal. These inputs are integrated with the beliefs created with perceptive
reasoning to facilitate effective decision-making and goal attainment. By
seamlessly processing the information through the described layered hierarchy,
our framework enhances the robot’s adaptability and responsiveness towards
various non-verbal user cues, while effectively interacting with objects on
the table and generating useful robot actions.

Figure 2: Our proposed system hierarchy for intention prediction.

###  III-A Non-verbal Cues

Building a human-aware robotic system requires perceiving not only user speech
but also user non-verbal cues [21], conveying human intentions. Therefore, it
is essential to enable our system to detect and interpret human gestures and
body language. To meet the real-time demands of robotic systems, we leverage
MediaPipe [22] to integrate various perception components. MediaPipe offers a
machine learning pipeline that facilitates a GPU integration, thanks to the
recently released open-source Tasks API. Although it was originally developed
with a focus on mobile devices, our experiments showcase its versatility and
cross-platform suitability for a Python integration with a robot in a tabletop
setup. Essentially, MediaPipe constructs a computer vision pipeline consisting
of directed graphs of modular components. These graph nodes manage the data
flow through data streams, configured with increasing timestamps. Each
processing node performs a specific calculation like image transformation,
limb detection, limb tracking or a machine learning inferential task using
frameworks like Tensorflow and Caffe, transforming the sensory input into
meaningful outputs like body landmarks. We construct three different pipelines
to extract various non-verbal cues:

#### Hand Pipeline

this pipeline is based on a single-shot Convolutional Neural Network (CNN) for
identifying hand regions, and a regression-based hand tracker, which predicts
landmarks within a cropped area around the detected hand. The model predicts
21 keypoints, derived from the GHUM [23] model for human shape and poses,
including four for each finger and one for the wrist. We utilize these
landmarks to train a recognition model for hand gestures, including pointing,
holding an object, and resting. Although the detector is capable of locating
hands in various distances, we distinguish between non-proximity and proximity
hand gestures, happening within a defined physical area around the robot’s
workstation. We use the hand 3D coordinates to identify these hand gestures
before LLM prompting.

#### Human-Pose Pipeline

an estimation of user poses is done in a top-down way by first detecting the
human in the frame with an optimized MobileNetV2. Then, poses are estimated
with 33 landmarks representing the joints with a light-weight CNN, called the
BlazePose [24] model. The extracted skeleton keypoints are mapped into 3D
coordinates of the torso, arms, legs, and head. We utilize the generated 3D
keypoints to track the user’s location in real-time relative to the robot.
Thus, we enable our system to lively interact with the user based on their
location in the workspace.

#### Face Pipeline

a face detector outputs 478 keypoints representing 3D face landmarks. The
pipeline consists of a BlazeFace model based on a Single Shot Detector [25]
(SSD), which identifies the main face region representing the user’s eyes,
mouth, nose, converted into a face mesh model. Also, a Blendshape model
predicts 52 scores representing various facial features, used to recognize the
user’s facial expressions. Thus, we enable the system to identify basic facial
expressions, such as smiling, and consequently the robot to react based on the
given user emotional state.

Figure 3: Perceptive reasoning of user non-verbal cues with examples.

###  III-B Perceptive Reasoning

As previously mentioned, we base our intention understanding mechanism on a
grounded LLM [20]. When such a system is booted, a set of initial prompts are
executed, responsible for providing the LLM with information about the robot,
its capabilities, surrounding environment, and purpose. The LLM is also
provided with a real-time stream of information about the objects on the
robot’s table using the ViLD [26] open-vocabulary object detector as well as
user utterances with an ASR system based on Whisper [27]. A set of robot
actions are generated like object and head manipulation based on the prompts
output. We use this LLM backbone for predicting the user’s intention by
collecting user verbal commands, non-verbal cues, and the environment state
and converting them into texts used to prompt the LLM.

Essential to this reasoning step is the system’s ability to perceive human
actions and gestures, thus being constantly “human aware”. For that, we
develop a so-called perceptive reasoning method, in which the system
interprets non-verbal cues, namely: hand gestures, body poses, and facial
expressions, and combine them with the environment state to make an
appropriate decision leading to a collaborative action. We follow a
hierarchical approach for parsing and interpreting the cues before LLM
prompting. First, the camera input passes to the various detectors (ex: hand
detector), where joint keypoints are extracted and then passed into a
recognition module (ex: gesture recognition). The recognition output provides
basic textual tokens about the non-verbal cues, e.g., the detected limb and
recognized movement. Later, these tokens pass into the LLM as internal prompts
invisible to the user to be converted into meaningful textual representations
of the user’s state. Fig. 3 highlights the steps with concrete examples. For
this step, we consider proximity gestures only.

###  III-C Task Reasoning

During user-robot collaborations, the system considers the user states
extracted in the perceptive reasoning step while acting to achieve a mutual
goal. Since we evaluate our system using object categorization, this task will
be the ultimate goal steering the interaction. Tasks are described in the form
of task prompts, which are exchangeable and expandable in a modular way. In
task reasoning, the system combines the user states, user explicit commands,
i.e., user queries, as well as any given task prompts in order to prompt the
LLM. If an additional goal is defined by the user, then the system will use
the given information to achieve it. Otherwise, the system will respond to the
user based on the available information as the user intentions and goals are
often vague and require multiple iterations to clarify (thus the dotted line
in Fig. 2).

Prompt III-C displays an example of a defined task prompt for categorizing
objects. By following the proposed hierarchical architecture, the system can
infer user goals and intentions in a categorization task (see Prompt III-C)
and during object interactions (see Prompt III-C). Nevertheless, our system
exhibits various additional useful features, such as establishing a joint
object attention between the user and robot, utilizing the LLM’s language
comprehension and generation skills to perform conversation repairs, identify
object ambiguities (see Prompt III-C), all while interacting naturally and
intuitively with hand gestures. The robot can also identify user proximity,
i.e., their location relative to the robot, thus initiating the interaction
and playing both active and passive roles, while reacting socially to user
facial expressions (cf. Prompt III-C).

Prompt 2: An example of a defined task prompt for categorizing objects.
SYSTEM: In this demo, you and the user should collaboratively split the
objects on the table into two categories. The location of the object on the
table represents its category. There are three different possibilities (right,
centre, left).
Category 1 represents all the objects on the ride side of the table, and
Category 2 represents all the objects on the left side of the table. Objects
in the centre do not belong to neither categories yet and are still
uncategorized.
You will not be told the exact categories explicitly but you should try to
infer that from the actions of the user.
The categories will be based on known properties of the objects like shape,
purpose, material, and color. As you understand better the categories that the
user is trying to create, you should help the user decide in which category an
object belongs. Prompt 3: The robot infers the user’s goal when categorizing
objects. \- the user placed the can on the right side of the table and apple
on the opposite -
Input:
CUE: The user is pointing to the banana
USER: Considering the properties and usage of this object, which of category 1
or category 2 is more fitting?
Output:
Considering that the banana is a fruit, it would be more fitting in Category
2, which currently contains other fruits.
ACTION: Move yellow banana to the left side of the table
0.4pt3pt 2pt Input:
What is your inference about the user intention from the previous
interactions?
Output:
From our previous interactions, I infer that the user intended to categorize
the objects on the table into two groups: containers and fruits. Prompt 4: The
robot reads user intention and anticipates the next action. Input:
CUE: The user holds the tomato soup can
USER: I am hungry
Output
ACTION: Give bowl to user
SAY: Would you like me to pour the tomato soup into the bowl for you? Prompt
5: Object disambiguation and intuitive interaction using pointing gestures.
Input:
CUE: The user points to the banana
Output:
ACTION: Look at banana
SAY: The banana is a fruit that is usually yellow when ripe. It is a good
source of potassium, vitamin C, and dietary fiber.
0.4pt3pt 2pt Input:
USER: Give me the yellow object
Output:
SAY: Sure, I’d be happy to help! Could you please clarify which object you’re
referring to as the yellow one? Is it the lemon or the banana?
0.4pt3pt 2pt Input:
USER: The user points to the banana
Output:
SAY: Ah, I see! You’re referring to the banana. Sure, I can give you the
banana. Here you go!
ACTION: Give banana to user
0.4pt3pt 2pt Input:
USER: Give me another object
Output:
SAY: Certainly, which object would you like me to give you? Just let me know
and I’ll be happy to assist you.
0.4pt3pt 2pt Input:
CUE: The user points to the lemon
Output:
SAY: Certainly, I can give you the lemon.
ACTION: Give lemon to user  Prompt 6: The robot detects user proximity and
initiates an interaction accordingly. Input:
CUE: The user has both hands rested on the table
Output:
ACTION: Look at user
SAY: Hello there! How can I assist you today?
0.4pt3pt 2pt
Input:
CUE: The user is pointing to the table, The user is smiling
USER: Hi
Output
ACTION: Express happiness
SAY: Hello! How can I assist you today? Would you like me to look at or point
to a specific object on the table? Or perhaps you have a question for me?
Prompt 7: Prompts used for system evaluation in an object categorization task.
Phase 1 (Task Understanding Check):
USER: What are the objects in category 1?
USER: What are the objects in category 2?
USER: What do objects in category 1 have in common?
USER: What do objects in category 2 have in common?
USER: Which objects are currently uncategorized?
Phase 2 (Object Categorization Check):
the following prompt is repeated four times:
CUE: The user is pointing to the object_name
USER: Considering the properties of this object, which of category 1 or
category 2 is more fitting?
Phase 3 (Intention Understanding Check):
USER: What do you infer about the user intention from the previous
interactions?
USER: Can you explain in detail all the choices you made for categorization?

##  IV Experiments and Results

Figure 4: A concrete workflow example of the object categorization task
generated by our system.

###  IV-A NICOL Robot & Environment Overview

We use NICOL [7], which is a semi-humanoid robot developed at the Knowledge
Technology Group of the University of Hamburg. NICOL is an adult-sized robot
fixed on a tabletop with two arms and a humanoid robotic head. Each arm has a
five-fingered hand, is capable of bio-inspired motion planning [28], and can
interact with adult-sized objects. The design of the NICOL robot allows
collaboration with humans through precise physical interactions with the
environment. Nevertheless, it also features social interaction with robot cues
like facial expressions, gestures, and gaze.

###  IV-B Object Categorization Task

TABLE I: Results of the object categorization task.  |  Model |  Fruit container |  Fruit fruit |  Container fruit |  Container container
---|---|---|---|---|---
Phase 1 |  gpt-3.5 |  0.99 |  0.68 |  0.91 |  0.88
gpt-3.5-16k |  0.62 |  0.74 |  0.89 |  0.81
|  gpt-4 |  0.95 |  0.84 |  0.9 |  0.98
|  vicuna |  0.33 |  0.33 |  0.42 |  0.39
|  mistral |  0.32 |  0.18 |  0.21 |  0.33
Phase 2 |  gpt-3.5 |  0.82 |  0.52 |  0.8 |  0.61
gpt-3.5-16k |  0.7 |  0.73 |  0.98 |  0.78
|  gpt-4 |  0.96 |  0.93 |  0.97 |  0.99
|  vicuna |  0.54 |  0.48 |  0.73 |  0.53
|  mistral |  0.725 |  0.78 |  0.35 |  0.66
Phase 3 |  gpt-3.5 |  0.71 |  0.37 |  0.71 |  0.7
gpt-3.5-16k |  0.75 |  0.53 |  0.89 |  0.63
|  gpt-4 |  0.92 |  0.75 |  0.92 |  0.93
|  vicuna |  0.17 |  0.56 |  0.48 |  0.29
|  mistral |  0.75 |  0.53 |  0.44 |  0.4

In our task, the user and robot collaboratively categorize objects on the
table into two distinct groups. The robot assists the user in finding the fit
for each “uncategorized” object based on common properties and
characteristics, like shape, color, and purpose. Unlike typical autonomous
object categorization tasks [29][30][31], the robot is not given an explicit
definition of the categories but should infer the goal based on user actions.
The task starts as the user places one object on each side of the table, i.e.,
representing Category 1 and Category 2. Then, the robot responds to user
queries regarding where to place a specific uncategorized object. An
explanatory example of the categorization task is in Fig. 4.

We evaluate the task using six objects: red apple, yellow banana, red can,
yellow lemon, red bowl, and red cup. We select the objects due to their
various purposes (ex: fruit and containers), colors (ex: red and yellow), and
shapes (ex: round and cylindrical). For performance comparison, we start the
task with every possible pair of objects as category initiators. After
category initialization, the robot is prompted to categorize a randomly
selected object. We conduct 150 trials in the experiment with 10 for each
possible object pair, each with a reset chat history. Moreover, we evaluate
the system using OpenAI and open-source LLM backbones, namely: ChatGPT
gpt-3.5-turbo-0301, gpt-3.5-turbo-16k-0613, gpt-4-0613, vicuna-1.3-33b, and
mistral-7b-instruct-v0.2. We follow a three-phase evaluation scheme, with a
total of 11 prompts (see Prompt III-C): 1) Task Understanding Check: five
prompts testing the LLM’s understanding of the task and categories at hand, 2)
Object Categorization Check: for each “uncategorized” object, a prompt is
given in the form of a categorization request (four prompts in total), and 3)
Intention Understanding Check, in which the LLM’s understanding of user
intentions and goals is tested with two prompts. The system is also expected
to provide logical explanations about its decisions and actions. A summary of
the results is in Table I with values representing model accuracy, averaged by
the number of trials and object pairs.

TABLE II: Pearson correlation of categorization phases (with p-values).  Model | Phase 1–Phase 2 | Phase 2–Phase 3 | Phase 1–Phase 3
---|---|---|---
gpt-3.5 | 0.8 (0.0003) | 0.6 (0.018) | 0.53 (0.044)
gpt-3.5-16k | 0.53 (0.04) | 0.51 (0.05) | 0.4 (0.14)
gpt-4 | 0.44 (0.1) | 0.79 (0.0005) | 0.63 (0.012)
vicuna | 0.53 (0.043) | 0.75 (0.001) | 0.49 (0.062)
mistral | 0.56 (0.029) | 0.21 (0.44) | 0.26 (0.35)

Our system exhibited proficiency in predicting user goals and deriving
intended categories of the various fruits and containers, irrespective of the
comparison order, especially gpt-4, with highly consistent accuracy. The
performance of the gpt-3.5 models was comparable and dropped slightly when
comparing objects with similar function (fruit vs. fruit or container vs.
container), indicating these models gave more importance to “usage” over other
attributes. Nonetheless, gpt-4 demonstrated competence and flexibility when
the categorizing criteria depended on subtle features like color or shape. All
OpenAI models showed good task realisation in Phase 1 and the occasional
accuracy drop is explained by the accidental mixing of the category location
(swapping right and left), hinting that the LLM spatial reasoning and
orientation have not yet progressed to a performant level.

The open-source models came short in Phase 1 while having near to above-
average accuracy in Phases 2 and 3. Mistral showed higher competence in Phase
2 than Vicuna despite the smaller size (7B vs. 33B parameters), which we
hypothesize to be a result of its own training, i.e., Mistral is fine-tuned on
instructions data making it more flexible with user commands, while Vicuna is
fine-tuned on conversation data. However, our results showed a positive
correlation between the model’s task understanding in Phase 1 and its
performance in subsequent phases (cf. Table II). This becomes clearer by
observing the error rates (see Fig. 5), which are lower for “Category
Reasoning” than grounding factors (category definition and uncategorized
objects). Also, all LLMs showed high success rates at interpreting non-verbal
cues through our perceptive reasoning. From the OpenAI models, gpt-4 reached
the best performance at Phase 3, while the other OpenAI models, despite making
correct decisions, were less efficient at precise decision explaining.
Although ChatGPT models showed higher resilience in grounding to a new task,
the results showcased that the highest error rates were related to not fully
grasping some task-related concepts rather than the categorization decision
itself, suggesting that the model’s reasoning aligns with the objectives of
the task. As the model develops a better understanding of the task over time,
its performance is expected to improve accordingly.

##  V Discussion and Conclusion

We explored LLMs for predicting human intention, showing their ability to
infer user goals and generate discrete actions to fulfill them in a human-
robot collaborative task. We proposed a hierarchical process for intention
prediction with our perceptive and task reasoning modules. We illustrated the
high capacity of LLMs to interpret non-verbal cues through our perceptive
reasoning. Our system was also able to adapt to a given task, i.e., object
categorization, in a few-shot manner, highlighting its potential for rapid
task acquisition and deployment in diverse settings, thus laying the
groundwork for effective intention prediction in a robotic collaborative
scenario. Future research endeavors will assess the LLM scalability across a
broader range of interactive scenarios to ascertain the transferability of its
intention prediction skills across multiple tasks. Our future work will also
examine the applicability of multimodal LLMs for inferring user gestures and
poses in real-time interaction with the robot.

Figure 5: Sources of system errors in the object categorization task.

## References

  * [1] E. Hildt, “What sort of robots do we want to interact with? reflecting on the human side of human-artificial intelligence interaction,” _Frontiers in Computer Science_ , vol. 3, 2021.
  * [2] Y. Zhang and T. Doyle, “Integrating intention-based systems in human-robot interaction: a scoping review of sensors, algorithms, and trust,” _Frontiers in Robotics and AI_ , vol. 10, 2023.
  * [3] A. Lubitz, L. Gutzeit, and F. Kirchner, “CoBaIR: A python library for context-based intention recognition in human-robot-interaction,” in _2023 32nd IEEE RO-MAN_ , 2023, pp. 2003–2009.
  * [4] K. Valmeekam, S. Sreedharan, M. Marquez, A. O. Hernandez, and S. Kambhampati, “On the planning abilities of large language models,” _ArXiv_ , vol. abs/2302.06706, 2023.
  * [5] B. Zhang and H. Soh, “Large language models as zero-shot human models for human-robot interaction,” _IEEE/RSJ IROS_ , pp. 7961–7968, 2023.
  * [6] S. Tong, Z. Liu, Y. Zhai, Y. Ma, Y. LeCun, and S. Xie, “Eyes wide shut? exploring the visual shortcomings of multimodal LLMs,” _ArXiv_ , vol. abs/2401.06209, 2024.
  * [7] M. Kerzel, P. Allgeuer, E. Strahl, N. Frick, J.-G. Habekost, M. Eppe, and S. Wermter, “NICOL: A neuro-inspired collaborative semi-humanoid robot that bridges social interaction and reliable manipulation,” _IEEE Access_ , vol. 11, pp. 123 531–123 542, 2023.
  * [8] S. Veselic, C. Zito, and D. Farina, “Human-robot interaction with robust prediction of movement intention surpasses manual control,” _Frontiers in Neurorobotics_ , vol. 15, 2021.
  * [9] Y. Yao, J. Duan, K. Xu, Y. Cai, Z. Sun, and Y. Zhang, “A survey on large language model (LLM) security and privacy: The good, the bad, and the ugly,” _ArXiv_ , vol. abs/2312.02003, 2024.
  * [10] B. AlKhamissi, M. Li, A. Celikyilmaz, M. Diab, and M. Ghazvininejad, “A review on language models as knowledge bases,” _ArXiv_ , vol. abs/2204.06031, 2022.
  * [11] J. Liu, C. Yu, J. Gao, Y. Xie, Q. Liao, Y. Wu, and Y. Wang, “LLM-powered hierarchical language agent for real-time Human-AI coordination,” _ArXiv_ , vol. abs/2312.15224, 2024.
  * [12] M. Jang, Y. Yoon, J. Choi, H. Ong, and J. Kim, “A structured prompting based on belief-desire-intention model for proactive and explainable task planning,” in _Proceedings of the 11th International Conference on Human-Agent Interaction_ , 2023, p. 375–377.
  * [13] M. A. Graule and V. Isler, “GG-LLM: Geometrically grounding large language models for zero-shot human activity forecasting in human-aware task planning,” _ArXiv_ , vol. abs/2310.20034, 2023.
  * [14] Q. Zhao, S. Wang, C. Zhang, C. Fu, M. Q. Do, N. Agarwal, K. Lee, and C. Sun, “AntGPT: Can large language models help long-term action anticipation from videos?” _ArXiv_ , vol. abs/2307.16368, 2023.
  * [15] M. J. Q. Zhang and E. Choi, “Clarify when necessary: Resolving ambiguity through interaction with LMs,” _ArXiv_ , vol. abs/2311.09469, 2023.
  * [16] T. Yoshida, A. Masumori, and T. Ikegami, “From text to motion: Grounding GPT-4 in a humanoid robot ”Alter3”,” _ArXiv_ , vol. abs/2312.06571, 2023.
  * [17] N. Cherakara, F. Varghese, and S. S. et al., “FurChat: An embodied conversational agent using LLMs, combining open and closed-domain dialogue with facial expressions,” _ArXiv_ , vol. abs/2308.15214, 2023.
  * [18] Y. K. Lee, Y. Jung, G. Kang, and S. Hahn, “Developing social robots with empathetic non-verbal cues using large language models,” _ArXiv_ , vol. abs/2308.16529, 2023.
  * [19] G. J. Serfaty, V. O. Barnard, and J. P. Salisbury, “Generative facial expressions and eye gaze behavior from prompts for multi-human-robot interaction,” in _Adjunct Proceedings of the 36th Annual ACM Symposium on User Interface Software and Technology_ , NY, USA, 2023.
  * [20] P. Allgeuer, H. Ali, and S. Wermter, “When robots get chatty: Grounding multimodal human-robot conversation and collaboration,” 2024, submitted to ICANN 2024.
  * [21] H. Ali, D. Jirak, and S. Wermter, “Snapture—a novel neural architecture for combined static and dynamic hand gesture recognition,” _Cognitive Computation_ , vol. 15, no. 6, pp. 2014–2033, Nov 2023.
  * [22] C. Lugaresi, J. Tang, H. Nash, C. McClanahan, E. Uboweja, M. Hays, F. Zhang, C.-L. Chang, M. G. Yong, J. Lee, and W.-T. C. et al., “MediaPipe: A framework for building perception pipelines,” _ArXiv_ , vol. abs/1906.08172, 2019.
  * [23] H. Xu, E. G. Bazavan, A. Zanfir, B. Freeman, R. Sukthankar, and C. Sminchisescu, “GHUM & GHUML: Generative 3D human shape and articulated pose models,” in _IEEE/CVF Conference on Computer Vision and Pattern Recognition (Oral)_ , 2020, pp. 6184–6193.
  * [24] V. Bazarevsky, I. Grishchenko, K. Raveendran, T. Zhu, F. Zhang, and M. Grundmann, “BlazePose: On-device real-time body pose tracking,” _ArXiv_ , vol. abs/2006.10204, 2020.
  * [25] W. Liu, D. Anguelov, D. Erhan, C. Szegedy, S. Reed, C.-Y. Fu, and A. C. Berg, “SSD: Single shot multibox detector,” in _Computer Vision – ECCV 2016_. Springer International Publishing, 2016, pp. 21–37.
  * [26] X. Gu, T.-Y. Lin, W. Kuo, and Y. Cui, “Open-vocabulary object detection via vision and language knowledge distillation,” _ArXiv_ , vol. abs/2104.13921, 2022.
  * [27] A. Radford, J. W. Kim, T. Xu, and G. B. et al., “Robust speech recognition via large-scale weak supervision,” in _International Conference on Machine Learning_ , 2023, pp. 28 492–28 518.
  * [28] J.-G. Habekost, E. Strahl, P. Allgeuer, M. Kerzel, and S. Wermter, “Cycleik: Neuro-inspired inverse kinematics,” in _ICANN 2023_ , Sep 2023, p. 457–470.
  * [29] M. Koskinopoulou, F. Raptopoulos, and G. P. et al., “Robotic waste sorting technology: Toward a vision-based categorization system for the industrial robotic separation of recyclable waste,” _IEEE Robotics & Automation Magazine_, vol. 28, no. 2, pp. 50–60, 2021.
  * [30] Z. Pan and K. Hauser, “Decision making in joint push-grasp action space for large-scale object sorting,” in _IEEE ICRA_ , 2021, pp. 6199–6205.
  * [31] S. Ramadurai and H. Jeong, “Effect of human involvement on work performance and fluency in human-robot collaboration for recycling,” in _ACM/IEEE International Conference on HRI_ , 2022, pp. 1007–1011.


S. DILL, M. ROHR, mnmDTW FOR CAMERA-BASED MOVEMENT ERROR LOCALIZATION

# mnmDTW: An extension to Dynamic Time Warping for Camera-based Movement Error
Localization

Sebastian DILL1  Maurice ROHR1  1KIS*MED – AI Systems in Medicine, Technische
Universität Darmstadt, Merckstraße 25, Darmstadt, Germany<EMAIL_ADDRESS>darmstadt.de<EMAIL_ADDRESS>
###### Abstract

In this approach, we use Computer Vision (CV) methods to extract pose
information out of exercise videos. We then employ Dynamic Time Warping (DTW)
to calculate the deviation from a gold standard execution of the exercise.
Specifically, we calculate the distance between each body part individually to
get a more precise measure for exercise accuracy. We can show that exercise
mistakes are clearly visible, identifiable and localizable through this
metric.

###### keywords:

pose estimation, movement analysis, camera, DTW

##  1 Introduction

Physical therapy is a crucial step in the treatment of many injuries and
diseases. One example is hemophilia, where physiotherapy and rehabilitation
can help prevent disabilities and preserve a patient’s autonomy [1]. While
ideally, physical therapy is performed under supervision of a medical
professional who can offer individual and immediate feedback, most people do
not have the resources to visit a training session regularly. Furthermore,
home exercises have been shown to be beneficial to the healing process [2]
even without supervision through an expert. On the other hand, wrong
executions, misjudgement of one’s fitness level, and overexertion might lead
to an inefficient training or even worse, serious injuries [3]. To mitigate
these problems, an automated evaluation system can be applied to assess the
quality of exercise execution and lessen the need for human supervision.
Building on the advances in computer vision in recent years, significant
research has been conducted on video-based human pose estimation and motion
capture. One of the most commonly used tools to extract pose data from videos
is MediaPipe Pose based on the BlazePose model [4].
Dynamic Time Warping (DTW) was established as a fundamental method to estimate
the distance between two time-series. While originally conceived for the one-
dimensional case, it has also been extended to the multi-dimensional DTW
(mDTW) [5]. As such, it has been applied in the evaluation of human movements,
for example in the works of Sempena et al. [6] and Adistambha et al. [7].
However, their approaches only evaluate the movement as a whole without
incorporating the causes of the error. This information is critical when
giving the patient feedback on how they should improve their exercises. Liu
and Chu [8] address this shortcoming and propose a camera-based machine
learning exercise evaluation system that assesses how well an exercise is
performed based on posture data extracted from videos, where they not only
identified the overall correctness of the exercises but also which body part
was responsible for the wrong posture. However, they do not apply DTW to their
approach, but instead use domain knowledge to create metrics that they then
feed to a deep learning network.
In this paper, we present an approach that has both the advantages of the DTW-
based solution (low computational complexity, high speed) and the ability to
evaluate the error for each body part individually. We propose a new multi-
layer normalized multi-dimensional DTW (mnmDTW) approach that is capable of
assessing individual body parts and test it on an example exercise. Our
evaluation shows very promising results as we can not only correctly classify
predefined movement execution errors but also describe them qualitatively in
terms of localization and type.

(a) From left to right: Correct movement; Mistake 1: not going low enough and
Mistake 2: having feet too wide.

(b) visibility is color-coded from blue (1) to green (0).

Figure 1: Visualization of (a) the exercise, where the lowest point of
movement is shown for all performed variations and
(b) the MediaPipe Pose output, consisting of x-y-camera coordinates of 33
different landmarks.

##  2 Experiment

In this paper, we focused on the squat exercise, since it is a stationary
movement that can be done without special equipment and has some defined
mistakes. All exercises were recorded by a single camera with $1920\times
1080$ pixels and 30 frames per second. We recorded RGB-videos of a single
participant performing the exercise 18 times, distinguished into three
different categories, depending on the execution quality. Ten executions were
considered _correct_ , which we defined by the participant’s feet being about
as wide as their shoulders and the minimum knee angle being close to
$90^{\circ}$. Two common mistakes are also considered and recorded four times
each. Mistake 1 is defined by the participant not going low enough and their
knee angle staying well above $90^{\circ}$. Mistake 2 is defined by the
participant’s feet being further apart than shoulder width. Figure 1(a) shows
the lowest point for one execution each of all three variations. From the ten
videos showing the correct exercise, one was randomly selected as the “gold
standard”. The remaining 17 videos are considered test videos.

##  3 Methods

For each test video, joint positions are extracted, normalized and
synchronised to the gold standard with a first mDTW. After combining the
joints to groups representing a limb each, a second mDTW is calculated to
receive group-specific mnmDTW values. The overall process can be seen in Fig.
2 and is explained in the following chapter.

Figure 2: The general idea behind mnmDTW for movement error localization.
After extracting and normalizing pose data from videos, a first mDTW is done
to synchronize the test video data of length $N$ with the predefined gold
standard recording of length $M$ to receive a time signal of length
$P=\text{max}(M,N)$. Then, the synchronized data is ordered in $L$ limb groups
of dimension $t_{l}$ and a second mDTW is performed for each group to receive
group-specific distance metrics.

###  3.1 Extraction of Pose Information

All 18 videos are first processed with the MediaPipe Pose library. The
MediaPipe output consists of x-y-z-camera-coordinates of 33 different pose
landmarks as well as an estimate for visibility between 0 and 1 for every
landmark. For this work, we only used the x-y-camera-coordinates. Since these
are given in pixels, with the origin being in the upper left corner of the
image, we need to normalize the coordinates with z-normalization to remove the
bias. For an overview over the provided landmarks, see Fig. 1(b).

###  3.2 Dynamic Time Warping

The general idea of DTW is to measure similarity between two time-series
$\boldsymbol{x}=[x_{1},\dots,x_{i},\dots,x_{M}],\leavevmode\nobreak\
\boldsymbol{y}=[y_{1},\dots,y_{j},\dots,y_{N}]$, that can have different
speeds or lengths. DTW is a non-linear algorithm that disregards the exact
timestamps at which observations occur. Instead, it finds the optimal ordering
of timestamps by minimizing the Euclidean distance between the series under
all admissible temporal alignments. Each alignment is characterized by its
alignment path
$\pi=[(i_{1},j_{1}),\dots,(i_{p},j_{p}),\dots,(i_{P},j_{P})],\leavevmode\nobreak\
P=\text{max}(M,N)$ mapping indices from one series to the other. For an
alignment to be admissible, its path has to fulfill these constraints:

  1. 1.

Each point from one series must be matched to at least one point from the
other, in a monotonically increasing way: $i_{p-1}\leq i_{p}\leq i_{p-1}+1$
and $j_{p-1}\leq j_{p}\leq j_{p-1}+1$.

  2. 2.

The first index from one series must be matched to the first index of the
other series: $\pi_{1}=(1,1)$. The same applies to the last indices:
$\pi_{P}=(M,N)$.

All mapping paths $\pi$ that satisfy these constraints span a set of possible
paths $\mathcal{A}$. The overall DTW error is then calculated as

$d_{\text{DTW}}(\boldsymbol{x},\boldsymbol{y})=\text{min}_{\pi\in\mathcal{A}(\boldsymbol{x},\boldsymbol{y})}\sqrt{\sum_{(i,j)\in\pi}d_{ij}^{2}},$
(1)

with the Euclidean distance $d_{ij}=\sqrt{(x_{i}-y_{j})^{2}}$.
In the case of K-dimensional time series, $\boldsymbol{x}$ and
$\boldsymbol{y}$ become matrices
$\boldsymbol{X}=[\boldsymbol{x}_{1},\dots,\boldsymbol{x}_{i},\dots,\boldsymbol{x}_{M}]$,
with $\boldsymbol{x}_{i}=[x_{i,1},\dots,x_{i,k},\dots,x_{i,K}]$ and
$\boldsymbol{Y}=[\boldsymbol{y}_{1},\dots,\boldsymbol{y}_{j},\dots,\boldsymbol{y}_{N}]$,
with $\boldsymbol{y}_{j}=[y_{j,1},\dots,y_{j,k},\dots,y_{j,K}]$. The distance
$d_{ij}$ then can be calculated as

$d_{ij}=\sqrt{\sum_{k}(x_{i,k}-y_{j,k})^{2}}.$ (2)

###  3.3 Multi-layer Normalized Multi-Dimensional Dynamic Time Warping

We calculate the multi-dimensional DTW (mDTW) as defined in equations 1 and 2
over all 66 dimensions (x-y-coordinates for 33 landmarks) to align all test
recordings to the gold standard. To evaluate which body part contributes the
most to the error, we then combine the landmarks into limb groups. The mapping
of landmarks to limb group can be seen in Table 1. Then, the mDTW distance
between the gold standard and the aligned recordings are calculated separately
both for each limb group and for x- and y-coordinates. The resulting mnmDTW
values are a metric for how similar specific body parts move in comparison to
the gold standard. Furthermore, the separation into x- and y-coordinates gives
us more information on the type of error.

landmark | limb group | dimension $t_{l}$
---|---|---
0-10 | head | 11
11, 12, 23, 24 | torso | 4
13, 15, 17, 19, 21 | left_arm | 5
14, 16, 18, 20, 22 | right_arm | 5
25, 27, 29, 31 | left_leg | 4
26, 28, 30, 32 | right_leg | 4
Table 1: Mapping of landmarks (see fig. 1(b)) to limb group.

For better comparability and interpretability, another five correct exercises
were selected as a control group. The corresponding videos were used to
calculate baseline mnmDTW values for each limb that were then averaged. All
other mnmDTW distances were normalized by dividing by the baseline average
value of the corresponding limb. This makes the mnmDTW values more intuitive.
Generally speaking, all values of 1 and below are regarded as _good_. On the
other hand, the less the movement matches the gold standard the greater the
mnmDTW value.

##  4 Results

Figures 3(a), 3(b) and 3(c) show mnmDTW values for three representative
example exercises, one from each class. As expected, for the _correct_
exercise, the values are approximately one for all limb groups. For mistake 1,
“having your feet too wide”, the values for the limb groups “left_leg” and
“right_leg” are far higher, especially in the x-axis, while the values
corresponding to the upper body parts remain close to one. Since this mistake
is defined solely by a horizontal offset of the feet, this observation matches
our expectations. It can also be seen that the mnmDTW value is even higher for
the left leg. This can be explained through the the projection of the movement
into camera-coordinates. As can be seen from Fig. 1(a), the right leg is
further away from the camera and therefore, the displacement is smaller in
camera-coordinates. The projection also explains why the left leg has a high
metric for the y-coordinate. For mistake 2, “not going deep enough”, all limbs
have an increased mnmDTW value, mainly in the y-axis. Furthermore, it can be
observed that the values are increasing with the limb’s height and the maximum
error is achieved for the head. Again, this matches the human observations
when looking at Fig. 1(a): Since the error is caused by not bending the knees
far enough, there is little difference in the positions of the feet and lower
legs. The further up the limbs are, the higher the offset gets, with
especially high errors for arms and head.

(a) correct

(b) Mistake 1: having your feet too wide.

(c) Mistake 2: not going deep enough.

Figure 3: mnmDTW values for three example exercises.

##  5 Discussion

For all of our recordings, the mnmDTW values accurately describe the visible
deviation from the gold standard movement. This not only enables a
classification based on predefined mistakes, but also, our results indicate
the method’s potential to qualitatively and quantitatively evaluate and
describe movements. It is possible not only to tell which limbs are
responsible for the error, but also, by looking at the coordinates separately,
to estimate where the error might be located and what it might look like. This
way, even unknown movement errors could be identified and corrected.

##  6 Conclusion and Outlook

This small-scale test shows very promising results for the mnmDTW approach. As
discussed, the metric not only allows for simple classifications into
predefined classes of movement quality, but also enables us to make
qualitative and quantitative statements about what the cause for low exercise
quality might be. However, our experiment was very limited in size and
variety. Overall, only 18 recordings were made of a single person from only
one camera angle. A larger experiment would be the logical next step. Here,
the influence of different looks, camera angles, fitness levels and exercises
can be researched. There are also several aspects how the mnmDTW metric could
be improved in the future. First, the metric could be extended to
3D-coordinates, which would reduce the expected dependency on the camera
properties and recording situation. Also, the degree to which the coordinates
are grouped into limbs, is a parameter of interest. Making the groups smaller
enables even more precise localization of the error, but is also expected to
introduce more noise to the metric.

## Acknowledgements

The research described in the paper was supervised by Prof. C. Hoog Antink, TU
Darmstadt.

## References

  * [1] HEIJNEN, L; BUZZARD, B. B. The role of physical therapy and rehabilitation in the management of hemophilia in developing countries. _Seminars in thrombosis and hemostasis_ , 2005, pp. 513-517.
  * [2] PROFFITT, R. Home exercise programs for adults with neurological injuries: A survey. _The American Journal of Occupational Therapy_ , 2016, pp. 1-8.
  * [3] JONES, B. H. et al. Exercise, training and injuries. _Sports Medicine_ , 1994, pp. 202-214.
  * [4] BAZAREVSKY, V. et al. Blazepose: On-device real-time body pose tracking. _arXiv preprint arXiv:2006.10204_ , 2020.
  * [5] SHOKOOHI-YEKTA, M. et al. Generalizing DTW to the multi-dimensional case requires an adaptive approach. _Data mining and knowledge discovery_ , 2017, pp. 1-31.
  * [6] SEMPENA, S. et al. Human action recognition using Dynamic Time Warping, _Proceedings of the 2011 International Conference on Electrical Engineering and Informatics_ , Bandung, Indonesia, 2011, pp. 1-5.
  * [7] ADISTAMBHA, K. et al. Motion classification using Dynamic Time Warping, _2008 IEEE 10th Workshop on Multimedia Signal Processing_ , Cairns, Australia, 2008, pp. 622-627.
  * [8] LIU, A; CHU, W. A posture evaluation system for fitness videos based on recurrent neural network. _2020 International Symposium on Computer, Consumer and Control (IS3C)_ , IEEE, 2020, pp. 185-188.

Sebastian DILL was born in München, Germany, in 1996, finished his Abitur in
2014 and then began his studies in Darmstadt. In 2021 he finished his master
degree in electrical engineering at the Technical University Darmstadt where
he is currently working towards a Ph.D. degree in electrical engineering at
the AI Systems in Medicine group. His research interests include non-obtrusive
pose estimation, movement analysis and physical therapy.

Maurice ROHR was born in Mannheim, Germany, in 1995, finished his Abitur in
2014 and then began his studies in Darmstadt. In 2020 he finished his master
degree in electrical engineering at the Technical University Darmstadt where
he is currently working towards a Ph.D. degree in electrical engineering at
the AI Systems in Medicine group. His research interests include medical
sensor fusion, imaging technologies and simulation.


# On ordered beta distribution and the generalized incomplete beta function

Mayad Al-Saidi111 Dept. of Mathematics and Statistics, York University, 4700
Keele Street, Toronto, ON, M3J 1P3, Canada. Email<EMAIL_ADDRESS>Alexey Kuznetsov222Dept. of Mathematics and Statistics, York University, 4700
Keele Street, Toronto, ON, M3J 1P3, Canada. Email<EMAIL_ADDRESS>and
Mikhail Nediak333School of Business, Queen’s University, 143 Union St.,
Kingston, ON, K7L 3N6, Canada. Email:
<EMAIL_ADDRESS>
###### Abstract

Motivated by applications in Bayesian analysis we introduce a multidimensional
beta distribution in an ordered simplex. We study properties of this
distribution and connect them with the generalized incomplete beta function.
This function is crucial in applications of multidimensional beta
distribution, thus we present two efficient numerical algorithms for computing
the generalized incomplete beta function, one based on Taylor series expansion
and another based on Chebyshev polynomials.

Keywords: Beta distribution, beta function, incomplete beta function, ordered
simplex, Chebyshev polynomials
2020 Mathematics Subject Classification : Primary 60E05, Secondary 65D15

##  1 Introduction

Before we introduce the main obect of our study – the ordered beta
distribution – let us consider the following motivating example. Let $S$ be a
random variable with binomial distribution with parameters $N$ and $X$, where
$N$ is a positive integer and $X$ is a random variable having
${\textnormal{Beta}}(\alpha,\beta)$ distribution. Then it is well-known that
the posterior distribution of $X$ given $S=m$ (where $0\leq m\leq N$) is
${\textnormal{Beta}}(\alpha+m,\beta+N-m)$. The ordered beta distribution
arises naturally in the following multidimensional generalization of the above
example. Let $\\{S_{i}\\}_{1\leq i\leq n}$ be random variables such that
$S_{i}$ has binomial $(N_{i},X_{i})$ distribution, and $\\{X_{i}\\}_{1\leq
i\leq n}$ are random variables such that

$0\leq X_{1}\leq X_{2}\leq\dots\leq X_{n}\leq 1$

almost surely. What prior distribution can we assign to the vector

${\mathbf{X}}=(X_{1},X_{2},\dots,X_{n})$

such that the posterior distribution of ${\mathbf{X}}$ given $S_{i}=m_{i}$ for
$i=1,2,\dots,n$ has the same form as the prior?

To answer this question, we first define the main object of our investigation.

######  Definition 1.

Let $n\in{\mathbb{N}}$ and $a_{i}>0$ and $b_{i}>0$ for $1\leq i\leq n$. We say
that a random vector ${\mathbf{X}}=(X_{1},X_{2},\dots,X_{n})$ has ordered beta
distribution with parameters $\\{a_{i}\\}_{1\leq i\leq n}$ and
$\\{b_{i}\\}_{1\leq i\leq n}$ if

${\mathbb{P}}(X_{1}\in{\textnormal{d}}x_{1},X_{2}\in{\textnormal{d}}x_{2},\dots,X_{n}\in{\textnormal{d}}x_{n})=C^{-1}{\mathbf{1}}_{\\{0\leq
x_{1}\leq x_{2}\leq\cdots\leq x_{n}\leq
1\\}}\prod\limits_{i=1}^{n}x_{i}^{a_{i}-1}(1-x_{i})^{b_{i}-1}{\textnormal{d}}x_{i}$
(1)

where $C$ is the normalization constant.

Before we proceed, we would like to explain the notation that will be used
everywhere in this paper. We denote random variables by capital letters (such
as $X,Y,Z$), real numbers by lower case letters (such as $x,y,z$) and vectors
by bold font (such as ${\mathbf{X}}$ or ${\mathbf{x}}$).

Definition 2 could be stated in the following equivalent way. Take $n$
independent random variables $Y_{i}\sim{\textrm{Beta}}(a_{i},b_{i})$ and
denote by $\Delta_{*}^{n}$ the ordered simplex

$\Delta_{*}^{n}:=\\{{\mathbf{x}}=(x_{1},\dots,x_{n})\in{\mathbb{R}}^{n}\;:\;0\leq
x_{1}\leq x_{2}\leq\dots\leq x_{n}\leq 1\\}.$ (2)

Then the ordered beta distribution can be defined as the distribution of the
vector ${\mathbf{Y}}=(Y_{1},\dots,Y_{n})$ conditioned on the event that it
lies in the ordered simplex $\Delta_{*}^{n}$:

${\mathbb{P}}({\mathbf{X}}\in{\textnormal{d}}{\mathbf{x}})={\mathbb{P}}({\mathbf{Y}}\in{\textnormal{d}}{\mathbf{x}}|{\mathbf{Y}}\in\Delta_{*}^{n}).$
(3)

The ordered beta distribution was introduced in [3] (also subsequently used in
[1]) in the study of dynamic pricing and demand learning problems. Let us
summarize the setup of these problems. Consider a retailer who wants to find a
price level at which selling a particular product would bring the highest
revenue. To solve this problem the retailer needs to know the demand for this
product as a function of price. One strategy to learn the demand is to fix a
set of $n$ ordered prices $p_{n}<p_{n-1}<\dots<p_{2}<p_{1}$ and then to
estimate the proportion of population that would buy the product at price
$p_{j}$ (for each $j$). Let $X_{j}\in(0,1)$ be the probability that a randomly
sampled customer would buy the product at price $p_{j}$. The demand is a
decreasing function of price, which implies that the probabilities $X_{j}$
must satisfy $0\leq X_{1}\leq X_{2}\leq\dots\leq X_{n}\leq 1$. We observe how
many customers made the decision to purchase or not to purchase the product
priced at $p_{j}$, and we denote by $m_{j}$ the number of successful sales at
price $p_{j}$ while $k_{j}$ is the number of times the product was offered at
price $p_{j}$ but the customer decided not to purchase it. We assume that the
customers make decisions independently of each other, thus the number of
customers who decided to purchase the product at price level $p_{j}$ had
binomial distribution with parameters $m_{j}+k_{j}$ and $X_{j}$. The question
is how the retailer can use the information contained in vectors
${\mathbf{m}}=(m_{1},\dots,m_{n})$ and ${\mathbf{k}}=(k_{1},\dots,k_{n})$ to
learn more about the probabilities $X_{j}$ and to increase the revenue? The
answer to this question can be found in [3], and it is based on the following
easily-verified property: if the prior distribution of purchase probabilities
${\mathbf{X}}=(X_{1},\dots,X_{n})$ is ordered beta with parameters
${\mathbf{a}}=(a_{1},\dots,a_{n})$ and ${\mathbf{b}}=(b_{1},\dots,b_{n})$,
then the posterior distribution of ${\mathbf{X}}$ given ${\mathbf{m}}$ and
${\mathbf{k}}$ is also ordered beta with parameters
${\mathbf{a}}+{\mathbf{m}}$ and ${\mathbf{b}}+{\mathbf{k}}$. This result can
also be found in [1] [Proposition 1, page 13].

The aim of this paper is to study the ordered beta distribution and to state
its properties and also to present several numerical algorithms which will
make it easier to use this distribution in applications. The paper is
organized as follows. In Section 2 we present various probabilistic properties
of ordered beta distribution (such as marginal distributions). To state these
properties, we need to introduce a new special function, which we call the
generalized incomplete beta function. The generalized incomplete beta function
is indispensable for applications of the ordered beta distribution, in
particular it is important to be able to compute this function numerically. In
Section 3 we study analytic properties of the generalized incomplete beta
function and in Section 4 we provide two efficient numerical algorithms for
its computation. We demonstrate the efficiency of these algorithms by
presenting the results of several numerical experiments in Section 5.

##  2 Properties of the ordered beta distribution

The following function will be needed for describing the properties of the
ordered beta distribution:

######  Definition 2.

Let $n\in{\mathbb{N}}$ and $a_{i}>0$, $b_{i}>0$ for $1\leq i\leq n$. The
generalized incomplete beta function (with parameters $\\{a_{i}\\}_{1\leq
i\leq n}$ and $\\{b_{i}\\}_{1\leq i\leq n}$) is defined for $z\in[0,1]$ as
follows:

$\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}=\int\limits{\mathbf{1}}_{\\{0\leq
x_{1}\leq x_{2}\leq\dots\leq x_{n}\leq
z\\}}\Big{[}\prod\limits_{i=1}^{n}x_{i}^{a_{i}-1}(1-x_{i})^{b_{i}-1}\Big{]}{\textnormal{d}}x_{1}\dots\textnormal{d}x_{n}.$
(4)

We will call $\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{)}:=\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}1\Big{)}$ the generalized beta function.

It is clear that the normalization constant in (1) must be given by
$C=\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{)}.$ When $n=1$ we recover the classical
incomplete beta function

$\mathrm{B}\Big{(}\begin{matrix}a_{1}\\\
b_{1}\end{matrix}\Big{|}z\Big{)}=\int_{0}^{z}x^{a_{1}-1}(1-x)^{b_{1}-1}{\textnormal{d}}x,$

and the classical beta function

$\mathrm{B}\Big{(}\begin{matrix}a_{1}\\\
b_{1}\end{matrix}\Big{)}=\int_{0}^{1}x^{a_{1}-1}(1-x)^{b_{1}-1}{\textnormal{d}}x=\frac{\Gamma(a_{1})\Gamma(b_{1})}{\Gamma(a_{1}+b_{1})}.$

The above two functions are usually denoted as $\mathrm{B}_{z}(a,b)$ and
$\mathrm{B}(a,b)$, see [5][Section 8.38].

Everywhere in this section we assume that $n\geq 2$ and that ${\mathbf{X}}$
has ordered beta distribution with parameters $\\{a_{i}\\}_{1\leq i\leq n}$
and $\\{b_{i}\\}_{1\leq i\leq n}$. To simplify the presentation of results, in
the rest of the paper we follow the convention that
$\mathrm{B}\Big{(}\begin{matrix}-\\\ -\end{matrix}\Big{|}z\Big{)}=1$ for all
$z$.

In the next result we collect several properties of the ordered beta
distriubution.

######  Theorem 1.

Let ${\mathbf{X}}=(X_{1},X_{2},\dots,X_{n})$ have ordered beta distribution
with parameters $\\{a_{i}\\}_{1\leq i\leq n}$ and $\\{b_{i}\\}_{1\leq i\leq
n}$ and let $C=\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{)}$. The following statements are true:

  * (i)

The random vector $\hat{\mathbf{X}}:=(1-X_{n},1-X_{n-1},\dots,1-X_{1})$ has
ordered beta distribution with parameters $(b_{n},b_{n-1},\dots,b_{1})$ and
$(a_{n},a_{n-1},\dots,a_{1})$.

  * (ii)

For $1\leq k\leq n-1$ and $z\in(0,1)$

${\mathbb{P}}(X_{k}\leq
z<X_{k+1})=C^{-1}\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{k}\\\
b_{1},\dots,b_{k}\end{matrix}\Big{|}z\Big{)}\mathrm{B}\Big{(}\begin{matrix}b_{n},b_{n-1},\dots,b_{k+1}\\\
a_{n},a_{n-1},\dots,a_{k+1}\end{matrix}\Big{|}1-z\Big{)}.$ (5)

  * (iii)

For $1\leq k\leq n$ and $x\in(0,1)$

${\mathbb{P}}(X_{k}\in{\textnormal{d}}x)=C^{-1}x^{a_{k}-1}(1-x)^{b_{k}-1}\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{k-1}\\\
b_{1},\dots,b_{k-1}\end{matrix}\Big{|}x\Big{)}\mathrm{B}\Big{(}\begin{matrix}b_{n},b_{n-1},\dots,b_{k+1}\\\
a_{n},a_{n-1},\dots,a_{k+1}\end{matrix}\Big{|}1-x\Big{)}{\textnormal{d}}x.$
(6)

  * (iv)

For $1\leq k\leq n$ and $z\in[0,1]$

$\displaystyle{\mathbb{P}}(X_{k}\leq
z)=\sum\limits_{j=k}^{n}C^{-1}\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{j}\\\
b_{1},\dots,b_{j}\end{matrix}\Big{|}z\Big{)}\mathrm{B}\Big{(}\begin{matrix}b_{n},b_{n-1},\dots,b_{j+1}\\\
a_{n},a_{n-1},\dots,a_{j+1}\end{matrix}\Big{|}1-z\Big{)},$ (7)

and

$\displaystyle{\mathbb{P}}(z<X_{k})$
$\displaystyle=\sum\limits_{j=0}^{k-1}C^{-1}\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{j}\\\
b_{1},\dots,b_{j}\end{matrix}\Big{|}z\Big{)}\mathrm{B}\Big{(}\begin{matrix}b_{n},b_{n-1},\dots,b_{j+1}\\\
a_{n},a_{n-1},\dots,a_{j+1}\end{matrix}\Big{|}1-z\Big{)}.$ (8)

  * (v)

If $\alpha>-a_{k}$ and $\beta>-b_{k}$ then

${\mathbb{E}}[X_{k}^{\alpha}(1-X_{k})^{\beta}]=C^{-1}\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{k}+\alpha,\dots
a_{n}\\\ b_{1},\dots,b_{k}+\beta,\dots,b_{n}\end{matrix}\Big{)},$ (9)

and, more generally, if $\alpha_{i}>-a_{i}$ and $\beta_{i}>-b_{i}$ for
$i=1,2,\dots,n$ then

${\mathbb{E}}\Big{[}\prod\limits_{i=1}^{n}X_{i}^{\alpha_{i}}(1-X_{k})^{\beta_{i}}]=C^{-1}\mathrm{B}\Big{(}\begin{matrix}a_{1}+\alpha_{1},\dots,a_{n}+\alpha_{n}\\\
b_{1}+\beta_{1},\dots,b_{n}+\beta_{n}\end{matrix}\Big{)}.$ (10)

###### Proof.

The result in item (i) follows directly from Definition 2. To prove the result
in item (ii) we first note the following identity: for $z\in(0,1)$

$\int_{[0,1]^{n}}{\mathbf{1}}_{\\{z\leq x_{1}\leq x_{2}\leq\cdots\leq
x_{n}\leq
1\\}}\prod\limits_{i=1}^{n}x_{i}^{a_{i}-1}(1-x_{i})^{b_{i}-1}{\textnormal{d}}x_{i}=\mathrm{B}\Big{(}\begin{matrix}b_{n},b_{n-1}\dots,b_{1}\\\
a_{n},a_{n-1},\dots,a_{1}\end{matrix}\Big{|}1-z\Big{)}.$ (11)

Formula (11) can be derived by applying a change of variables $x_{i}\mapsto
1-y_{i}$ and using (4). Now we compute with the help of (4) and (11)

$\displaystyle{\mathbb{P}}(X_{k}\leq z<X_{k+1})$
$\displaystyle=\int_{[0,1]^{n}}C^{-1}{\mathbf{1}}_{\\{0\leq x_{1}\leq
x_{2}\leq\cdots\leq x_{n}\leq
1\\}}{\mathbf{1}}_{\\{x_{k}<z<x_{k+1}\\}}\prod\limits_{i=1}^{n}x_{i}^{a_{i}-1}(1-x_{i})^{b_{i}-1}{\textnormal{d}}x_{i}$
$\displaystyle=C^{-1}\int_{[0,1]^{k}}{\mathbf{1}}_{\\{0\leq x_{1}\leq
x_{2}\leq\cdots\leq
x_{k}<z\\}}\prod\limits_{i=1}^{k}x_{i}^{a_{i}-1}(1-x_{i})^{b_{i}-1}{\textnormal{d}}x_{i}$
$\displaystyle\;\;\;\;\;\;\times\int_{[0,1]^{n-k}}{\mathbf{1}}_{\\{z<x_{k+1}\leq\cdots\leq
x_{n}\leq
1\\}}\prod\limits_{i=k+1}^{n}x_{i}^{a_{i}-1}(1-x_{i})^{b_{i}-1}{\textnormal{d}}x_{i}$
$\displaystyle=C^{-1}\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{k}\\\
b_{1},\dots,b_{k}\end{matrix}\Big{|}z\Big{)}\mathrm{B}\Big{(}\begin{matrix}b_{n},b_{n-1},\dots,b_{k+1}\\\
a_{n},a_{n-1},\dots,a_{k+1}\end{matrix}\Big{|}1-z\Big{)}.$

The derivation of the marginal distribution of $X_{k}$ given in (6) follows in
the same way from (1) by integrating over variables all $x_{i}$ with $i\neq k$
and using (4) and (11) .

Formula (7) follows from (5) by noting that for an ordered beta distributed
vector ${\mathbf{X}}$ the event $\\{X_{k}\leq z\\}$ is the disjoint union of
events $\\{X_{j}\leq z<X_{j+1}\\}$, $k\leq j\leq n$. An identical argument is
used to prove (8).

Formulas (9) and (10) follow directly from (1) and (4).
$\sqcap\kern-8.0pt\hbox{$\sqcup$}$

##  3 Properties of the generalized incomplete Beta function

In the next proposition we collect several properties of the generalized
incomplete beta function, which will be useful later.

######  Proposition 1.

Assume that $z\in(0,1)$ and $a_{i}>0$, $b_{i}>0$ for $1\leq i\leq n$. Then the
following statements are true:

  * (i)

$\mathrm{B}\Big{(}\begin{matrix}a_{1},a_{2},\dots,a_{n}\\\
b_{1},b_{2},\dots,b_{n}\end{matrix}\Big{)}=\mathrm{B}\Big{(}\begin{matrix}b_{n},b_{n-1},\dots,b_{1}\\\
a_{n},a_{n-1},\dots,a_{1}\end{matrix}\Big{)},$ (12)

  * (ii)

$\displaystyle\sum\limits_{k=0}^{n}\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{k}\\\
b_{1},\dots,b_{k}\end{matrix}\Big{|}z\Big{)}\mathrm{B}\Big{(}\begin{matrix}b_{n},b_{n-1},\dots,b_{k+1}\\\
a_{n},a_{n-1},\dots,a_{k+1}\end{matrix}\Big{|}1-z\Big{)}=\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{)},$ (13)

  * (iii)

$\displaystyle\sum\limits_{k=0}^{n}(-1)^{k}\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{k}\\\
b_{1},\dots,b_{k}\end{matrix}\Big{|}z\Big{)}\mathrm{B}\Big{(}\begin{matrix}a_{n},a_{n-1},\dots,a_{k+1}\\\
b_{n},b_{n-1},\dots,b_{k+1}\end{matrix}\Big{|}z\Big{)}=0,$ (14)

  * (iv)

$\mathrm{B}\Big{(}\begin{matrix}a_{1},a_{2},\dots,a_{n}\\\
b_{1},b_{2},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}=\int_{0}^{z}x^{a_{n}-1}(1-x)^{b_{n}-1}\mathrm{B}\Big{(}\begin{matrix}a_{1},a_{2},\dots,a_{n-1}\\\
b_{1},b_{2},\dots,b_{n-1}\end{matrix}\Big{|}x\Big{)}{\textnormal{d}}x,$ (15)

  * (v)

for $1\leq k\leq n$

$\mathrm{B}\Big{(}\begin{matrix}a_{1},a_{2},\dots,a_{n}\\\
b_{1},b_{2},\dots,b_{n}\end{matrix}\Big{)}=\int_{0}^{1}x^{a_{k}-1}(1-x)^{b_{k}-1}\mathrm{B}\Big{(}\begin{matrix}a_{1},a_{2},\dots,a_{k-1}\\\
b_{1},b_{2},\dots,b_{k-1}\end{matrix}\Big{|}x\Big{)}\mathrm{B}\Big{(}\begin{matrix}b_{n},b_{n-1},\dots,b_{k+1}\\\
a_{n},a_{n-1},\dots,a_{k+1}\end{matrix}\Big{|}1-x\Big{)}{\textnormal{d}}x.$
(16)

###### Proof.

Formula (12) follows from (4) by setting $z=1$ and changing the variables of
integration $x_{i}\mapsto 1-y_{i}$. The identity (13) is obtained from (7) and
(8) by noting that ${\mathbb{P}}(X_{k}\leq z)+{\mathbb{P}}(z<X_{k})=1$.

To prove formula (14) one would need to multiply both sides of the following
identity

$\displaystyle{\mathbf{1}}_{\\{x_{1}\leq x_{2}\leq\cdots\leq
x_{n}\\}}-{\mathbf{1}}_{\\{x_{1}\leq x_{2}\leq\cdots\leq
x_{n-1}\\}}+{\mathbf{1}}_{\\{x_{1}\leq x_{2}\leq\cdots\leq
x_{n-2}\\}}{\mathbf{1}}_{\\{x_{n}<x_{n-1}\\}}$ (17)
$\displaystyle-{\mathbf{1}}_{\\{x_{1}\leq x_{2}\leq\cdots\leq
x_{n-3}\\}}{\mathbf{1}}_{\\{x_{n}<x_{n-1}<x_{n-2}\\}}+\dots+(-1)^{n}{\mathbf{1}}_{\\{x_{n}<x_{n-1}<\cdots<x_{1}\\}}=0,$

by

$\prod\limits_{i=1}^{n}x_{i}^{a_{i}-1}(1-x_{i})^{b_{i}-1}{\textnormal{d}}x_{i}$

and integrate over $[0,z]^{n}$. The identity (17) can be established easily by
induction; the case $n=2$ is equivalent to an obvious statement

${\mathbf{1}}_{\\{x_{1}\leq x_{2}\\}}+{\mathbf{1}}_{\\{x_{2}<x_{1}\\}}=1.$

Formula (15) follows from the definition (4) and formula (16) is a corollary
of (6). $\sqcap\kern-8.0pt\hbox{$\sqcup$}$

For $a_{i}>0$, $b_{i}>0$ and $z\in(0,1]$ we define

$\beta\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}:=z^{-a_{1}-a_{2}-\dots-
a_{n}}\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)},$ (18)

This function has better analytical properties, compared with the generalized
incomplete Beta function.

######  Proposition 2.

The function $z\mapsto\beta\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}$ is analytic in
${\mathbb{C}}\setminus[1,\infty)$ and satisfies

$\beta\Big{(}\begin{matrix}a_{1},a_{2},\dots,a_{n}\\\
b_{1},b_{2},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}=\int_{0}^{1}y^{a_{1}+a_{2}+\dots+a_{n}-1}(1-yz)^{b_{n}-1}\beta\Big{(}\begin{matrix}a_{1},a_{2},\dots,a_{n-1}\\\
b_{1},b_{2},\dots,b_{n-1}\end{matrix}\Big{|}yz\Big{)}{\textnormal{d}}y.$ (19)

If $b_{i}\in{\mathbb{N}}$ for all $1\leq i\leq n$, then
$\beta\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}$ is a polynomial in $z$.

###### Proof.

For $z\in(0,1)$ formula (19) follows from (15) by a change of variables
$x=zy$. The fact that $\beta\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}$ is analytic in
${\mathbb{C}}\setminus[1,\infty)$ follows from (19) by induction on $n$ (note
that for every $y\in(0,1)$ the function $z\mapsto(1-yz)^{b_{n}-1}$ is analytic
in ${\mathbb{C}}\setminus[1,\infty)$). The fact that $\beta$-function is a
polynomial in $z$ when the coefficients $b_{i}$ are integers follows from (19)
by induction in $n$. $\sqcap\kern-8.0pt\hbox{$\sqcup$}$

It is clear from Proposition 2 that the generalized incomplete beta function
$\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}$ is a polynomial in $z$ if
$a_{1}+a_{2}+\dots+a_{n}$ and all coefficients $b_{i}$ are positive integers
and that in the general case it is an analytic function in the complex plane
with two cuts along $(-\infty,0]$ and $[1,\infty)$.

##  4 Computing the generalized incomplete Beta function

Now we turn our attention to the question of computing the generalized
incomplete beta function for arbitrary $z\in(0,1]$. Our first observation is
that it is enough to be able to compute this function for $z\in(0,1/2]$. To
show this, we rewrite the identity (13) in the form

$\displaystyle\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{)}$
$\displaystyle=\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}+\mathrm{B}\Big{(}\begin{matrix}b_{n},b_{n-1},\dots,b_{1}\\\
a_{n},a_{n-1},\dots,a_{1}\end{matrix}\Big{|}1-z\Big{)}$ (20)
$\displaystyle+\sum\limits_{k=1}^{n-1}\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{k}\\\
b_{1},\dots,b_{k}\end{matrix}\Big{|}z\Big{)}\mathrm{B}\Big{(}\begin{matrix}b_{n},b_{n-1},\dots,b_{k+1}\\\
a_{n},a_{n-1},\dots,a_{k+1}\end{matrix}\Big{|}1-z\Big{)}.$

We see that we can compute the value of
$\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{)}$ by setting $z=1/2$ in (20) and when
$z\in(1/2,1)$ we can compute
$\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}$ by induction, since all other
terms in (20) either depend on $1-z\in(0,1/2)$ or have fewer parameters. Thus,
one can use identity (20) coupled with induction on $n$ to compute the value
of any generalized incomplete beta function (with arbitrary number of
parameters) for any value of $z\in[0,1]$.

Instead of computing the generalized incomplete beta function
$B\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}$, it is more convenient to solve
an equivalent problem of computing
$\beta\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}$, as this function is analytic in
a wider domain. Next we present two algorithms for computing this function:
the first algorithm is based on Taylor series expansion and the second
algorithm is based on Chebyshev series.

###  4.1 The Taylor expansion method

The computation proceeds by iteration of (15). First of all, we note that when
$n=1$ we have a series representation

$\beta\Big{(}\begin{matrix}a_{1}\\\
b_{1}\end{matrix}\Big{|}z\Big{)}=\int_{0}^{1}y^{a_{1}-1}(1-yz)^{b_{1}-1}{\textnormal{d}}y=\sum\limits_{k\geq
0}\frac{(1-b_{1})_{k}}{(a_{1}+k)k!}z^{k},\;\;\;|z|<1,$ (21)

which is easily obtained by integrating the binomial series expansion of the
term $(1-yz)^{b_{1}}$. Here $(a)_{k}=a(a+1)\dots(a+k-1)$ denotes the
Pochhammer symbol. Computing the coefficients of Taylor series of
$\beta\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}$ when $n\geq 2$ can be done
iteratively with the help of the following result.

######  Proposition 3.

For $n\geq 1$ and $k\geq 0$ denote by $c_{k}^{(n)}$ the coefficients in the
Taylor series expansions

$\beta\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}=\sum\limits_{k\geq
0}c^{(n)}_{k}z^{k},\;\;\;|z|<1.$ (22)

Then for $n\geq 2$ and $k\geq 0$

$c^{(n)}_{k}=\frac{1}{k+a_{1}+\dots+a_{n}}\sum\limits_{l=0}^{k}c^{(n-1)}_{k-l}\times\frac{(1-b_{n})_{l}}{l!}.$
(23)

###### Proof.

Using the Binomial series and (15) we obtain for $|z|<1$

$\displaystyle\beta\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}$
$\displaystyle=\int_{0}^{1}y^{a_{1}+\dots+a_{n}-1}(1-yz)^{b_{n}-1}\beta\Big{(}\begin{matrix}a_{1},\dots,a_{n-1}\\\
b_{1},\dots,b_{n-1}\end{matrix}\Big{|}yz\Big{)}{\textnormal{d}}y$
$\displaystyle=\int_{0}^{1}y^{a_{1}+\dots+a_{n}-1}\Big{[}\sum\limits_{l\geq
0}\frac{(1-b_{n})_{l}}{l!}(yz)^{l}\Big{]}\times\Big{[}\sum\limits_{k\geq
0}c_{k}^{(n-1)}(yz)^{k}\Big{]}{\textnormal{d}}y$
$\displaystyle=\int_{0}^{1}y^{a_{1}+\dots+a_{n}-1}\sum\limits_{k\geq
0}\Big{[}\sum\limits_{l=0}^{k}c^{(n-1)}_{k-l}\times\frac{(1-b_{n})_{l}}{l!}\Big{]}(yz)^{k}{\textnormal{d}}y$
$\displaystyle=\sum\limits_{k\geq
0}\frac{z^{k}}{k+a_{1}+\dots+a_{n}}\Big{[}\sum\limits_{l=0}^{k}c^{(n-1)}_{k-l}\times\frac{(1-b_{n})_{l}}{l!}\Big{]}.$

$\sqcap\kern-8.0pt\hbox{$\sqcup$}$

Now we present the algorithm for computing
$\beta\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}$ for $z\in(0,1/2]$.

Step 1: Fix a large integer $N$ and compute $\\{c_{k}^{(1)}\\}_{0\leq k\leq
N}$ via (21).

Step 2: For $m=2,3,\dots,n$, once we have the values of
$\\{c_{k}^{(m-1)}\\}_{0\leq k\leq N}$ we compute the values of
$\\{c_{k}^{(m)}\\}_{0\leq k\leq N}$ via (23). This computation involves a
convolution of two sequences, thus each step can be made more efficient with
the help of the Fast Fourier Transform.

Step 3: Compute the value of $\beta\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}$ by truncating the series in (21)
after $N$ terms.

Several remarks are in order. First of all, we need to discuss the error due
to truncation of series in (21). Since the function
$\beta\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}$ is analytic in
${\mathbb{C}}\setminus[1,\infty)$ and is continuous as $z\to 1$, for every $n$
we have $c_{k}^{(n)}\to 0$ as $k\to+\infty$. This fact can be seen by
indentifying these coefficients with the coefficients of Fourier series of a
continuous function

$t\mapsto\beta\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}e^{2\pi{\textnormal{i}}t}\Big{)}.$

Thus after truncating the series in (21) after $N$ terms we will have an error
of $O(|z|^{N})$, and since $z\in(0,1/2]$ this error is less than $C2^{-N}$,
for some constant $C$, which may depend on $a_{i}$ and $b_{i}$. Second, let us
discuss the computational complexity of this algorithm. Using the Fast Fourier
Transform, each Step 2 can be computed in $O(N\ln(N))$ arithmetic operations.
Thus, the computational complexity of evaluating
$\beta\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}$ is $O(nN\ln(N))$ arithmetic
operations.

The previous two comments show that the algorithm based on Taylor series is
fast-convergent and very efficient. However, there is one potential problem
with this algorithm. When one or more of parameters $a_{i}$ (or $b_{i}$) is
large, the coefficients $(1-a_{i})_{l}$ (or $(1-b_{i})_{l}$) that will appear
in (23) will also become very large and will have alternating sign, causing
lot of cancellation in the sum in (23) and resulting in loss of precision.
Therefore, when some of the parameters $a_{i}$, $b_{i}$ are large one has to
be mindful of this potential loss of precision. We resolved this problem by
using a multi-precision arithmetic when computing the coefficients
$c_{k}^{(n)}$ via (23).

###  4.2 The Chebyshev expansion method

Let $\\{\xi^{(n)}_{k}\\}_{0\leq k\leq N}$ be the coefficients in the Chebyshev
expansion

$\beta\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{|}z\Big{)}=\frac{1}{2}\xi^{(n)}_{0}+\sum\limits_{k\geq
1}\xi^{(n)}_{k}T_{k}(4z-1).$ (24)

Note that we use scaled Chebyshev polynomials $T_{k}(4z-1)$, since we are
interested in computing the $\beta$-function only for $z\in(0,1/2]$ and the
function $4z-1$ maps the interval $(0,1/2]$ onto $(-1,1)$. Below we present an
algorithm that computes the coefficients $\\{\xi^{(m)}_{k}\\}_{0\leq k\leq N}$
using the previously computed values of $\\{\xi^{(m-1)}_{k}\\}_{0\leq k\leq
N}$.

Initial step: Choose a large integer $N$ and set $m=1$. Let
$\hat{\xi}^{(0)}_{0}=2$ and $\hat{\xi}^{(0)}_{k}=0$ for $1\leq k\leq N$.

Step 1: For $j=0,1,\dots,N$ compute (using the Fast Fourier Transform)

$v_{j}=\frac{1}{2}\hat{\xi}^{(m-1)}_{0}+\sum\limits_{k=1}^{N}\hat{\xi}^{(m-1)}_{k}\cos\Big{(}\frac{\pi
k(j+1/2)}{N+1}\Big{)}.$ (25)

Step 2: For $k=0,1,\dots,N$ compute (using the Fast Fourier Transform)

$\eta_{k}=\frac{2}{N+1}\sum\limits_{j=0}^{N}(1-(1+x_{j})/4)^{b_{m}-1}v_{j}\cos\Big{(}\frac{\pi
k(j+1/2)}{N+1}\Big{)},$ (26)

where

$x_{j}:=\cos\Big{(}\frac{\pi(j+1/2)}{N+1}\Big{)}.$

Step 3: We set $\mu_{N}=\mu_{N+1}=0$ and $A_{m}:=a_{1}+a_{2}+\dots+a_{m}$ and
compute recursively

$\mu_{k-1}=\frac{1}{1+A_{m}/k}\Big{(}8\eta_{k}-2\mu_{k}-\mu_{k+1}(1-A_{m}/k)\Big{)},\;\;\;k=N,N-1,\dots,1,$
(27)

Step 4: We compute

$\hat{\xi}^{(m)}_{0}=\frac{1}{A_{m}}\Big{(}\eta_{0}-\frac{1}{4}(\mu_{0}+\mu_{1})\Big{)}$
(28)

and

$\hat{\xi}^{(m)}_{k}=\frac{\mu_{k-1}-\mu_{k+1}}{8k},\;\;\;k=1,2,\dots,N.$ (29)

After running the above algorithm for $m=1,2,\dots,n$, we obtain the
coefficients $\\{\hat{\xi}^{(n)}_{k}\\}_{0\leq k\leq N}$ (which also depend on
$N$) that serve as approximations to the desired coefficients $\xi^{(n)}_{k}$.
It is clear that the computational complexity of this algorithm is also
$O(nN\ln(N))$ – the same as for Taylor expansion method. The next theorem
gives us a bound for the error of the approximation.

######  Theorem 2.

There exists a constant $C$ (depending on $a_{i}$ and $b_{i}$) such that for
all $N\geq 1$

$\max\limits_{0\leq k\leq N}|\hat{\xi}^{(n)}_{k}-\xi^{(n)}_{k}|<C3^{-8N/5}.$
(30)

######  Remark 1.

The constant in (30) can be slightly improved: as will be clear from the proof
of Theorem 2, it is true that

$\max\limits_{0\leq k\leq N}|\hat{\xi}^{(n)}_{k}-\xi^{(n)}_{k}|=O(\rho^{-N})$

for every $\rho\in(0,3+8\sqrt{8})$, where the constant in the big-O notation
depends on $a_{i}$, $b_{i}$ and $\rho$. Note that $3^{8/5}=5.7995...$ and
$3+\sqrt{8}=5.8284...$, thus the bound given in (30), while suboptimal, should
be good enough for all practical purposes.

The proof of Theorem 2 relies on the following

Fact: Assume the function $g(z)$ is analytic in
${\mathbb{C}}\setminus[1,\infty)$. Then the coefficients in the Chebyshev
expansion

$g(z)=w_{0}/2+\sum\limits_{k\geq 1}w_{k}T_{k}(4z-1),\;\;\;0\leq z\leq 1/2,$

satisfy $w_{k}=O(\rho^{-k})$ for any $\rho\in(0,3+\sqrt{8})$.

The above result follows easily from the classical fact that the Chebyshev
series expansion of a function $f$ converges in the largest Bernstein ellipse
in which the function $f$ is analytic (see [4][Theorem 5.16]). More precisely,
the function $f((1+z)/4)$ is analytic in ${\mathbb{C}}\setminus[3,\infty)$,
thus it is analytic inside the Bernstein ellipse

$E_{\rho}:=\\{z\in{\mathbb{C}}:\frac{\textnormal{Re}(z)^{2}}{a_{\rho}^{2}}+\frac{\textnormal{Im}(z)^{2}}{b_{\rho}^{2}}=1\\}$

with $a_{\rho}=(\rho+\rho^{-1})/2$ and $b_{\rho}=(\rho-\rho^{-1})/2$ with
$\rho=3+\sqrt{8}$.

Proof of Theorem 2: The choice of $\hat{\xi}^{(0)}_{k}$ in the initial step
corresponds to our convention that

$\beta\Big{(}\begin{matrix}-\\\ -\end{matrix}\Big{|}z\Big{)}\equiv
1=\frac{1}{2}\hat{\xi}^{(0)}_{0}+\sum\limits_{k\geq
1}\hat{\xi}^{(0)}_{k}T_{k}(4z-1),$

thus at this level there is no approximation error:
$\xi^{(0)}_{k}=\hat{\xi}^{(0)}_{k}$ for all $k$.

Let us now justify Steps 1-4. We denote

$f(z):=\beta\Big{(}\begin{matrix}a_{1},\dots,a_{m-1}\\\
b_{1},\dots,b_{m-1}\end{matrix}\Big{|}z\Big{)},\;\;\;F(z):=\beta\Big{(}\begin{matrix}a_{1},\dots,a_{m}\\\
b_{1},\dots,b_{m}\end{matrix}\Big{|}z\Big{)}.$

According to (15) and (18), these two functions are related by the equation

$z^{A_{m}}F(z)=\int_{0}^{z}x^{A_{m}-1}(1-x)^{b_{m}-1}f(x)dx,\;\;\;z\in(0,1),$

where $A_{m}:=a_{1}+a_{2}+\dots+a_{m}$. We take derivative of both sides of
the above equation and obtain a differential equation

$A_{m}F(z)+zF^{\prime}(z)=(1-z)^{b_{m}-1}f(z),\;\;\;z\in(0,1).$ (31)

Assume we know the coefficients $\theta_{k}$ in the Chebyshev expansion

$(1-z)^{b_{m}-1}f(z)=\theta_{0}/2+\sum\limits_{k\geq
1}\theta_{k}T_{k}(4z-1),\;\;\;0\leq z\leq 1/2.$ (32)

Let us denote by $\nu_{k}$ the coefficients in the Chebyshev expansion

$F^{\prime}(z)=\nu_{0}/2+\sum\limits_{k\geq 1}\nu_{k}T_{k}(4z-1),\;\;\;0\leq
z\leq 1/2.$ (33)

Chebyshev polynomials satisfy

$\int T_{n}(x)dx=\frac{T_{n+1}(x)}{2(n+1)}-\frac{T_{n-1}(x)}{2(n-1)},$

which can be established from formulas 8.941 and 8.949.1 in [5]. From here we
obtain

$\int
T_{k}(4z-1)dz=\frac{1}{8}\Big{(}\frac{T_{k+1}(4z-1)}{k+1}-\frac{T_{k-1}(4z-1)}{k-1}\Big{)},\;\;\;k\geq
2.$

We also check directly (using the facts $T_{1}(x)=x$ and $T_{2}(x)=2x^{2}-1$)
that

$\int 1dz=T_{1}(4z-1)/4,\;\;\;\int T_{1}(4z-1)dz=T_{2}(4z-1)/16.$

Combining the above formulas, we obtain, after integrating both sides of (33)
and rearranging the terms of the series, the following expansion

$F(z)=Q+\sum\limits_{k\geq
1}\frac{\nu_{k-1}-\nu_{k+1}}{8k}T_{k}(4z-1),\;\;\;0\leq z\leq 1/2,$ (34)

for some constant $Q$. Next we use formulas

$\displaystyle z=\frac{1}{4}(T_{1}(4z-1)+1),$ $\displaystyle
zT_{k}(4z-1)=\frac{1}{8}(T_{k+1}(4z-1)+2T_{k}(4z-1)+T_{k-1}(4z-1)),\;\;\;k\geq
1,$

and equations (31), (32), (33) and (34) and conclude that

$\displaystyle A_{m}\Big{[}Q+\sum\limits_{k\geq
1}\frac{\nu_{k-1}-\nu_{k+1}}{8k}T_{k}(4z-1)\Big{]}+\frac{\nu_{0}}{8}\Big{(}T_{1}(4z-1)+1\Big{)}$
$\displaystyle+\sum\limits_{k\geq
1}\frac{\nu_{k}}{8}\Big{(}T_{k+1}(4z-1)+2T_{k}(4z-1)+T_{k-1}(4z-1)\Big{)}=\theta_{0}/2+\sum\limits_{k\geq
1}\theta_{k}T_{k}(4z-1).$

Comparing the coefficients in front of $T_{k}$ in the above equation gives us
the following identities:

$\displaystyle 2A_{m}Q+\frac{1}{4}(\nu_{0}+\nu_{1})=\theta_{0},$ (35)
$\displaystyle
A_{m}\frac{\nu_{k-1}-\nu_{k+1}}{8k}+\frac{1}{8}(\nu_{k-1}+2\nu_{k}+\nu_{k+1})=\theta_{k},\;\;\;k\geq
1.$ (36)

The equation (36) can be written in an equivalent form

$\nu_{k-1}=\frac{1}{1+\frac{A_{m}}{k}}\Big{(}8\theta_{k}-2\nu_{k}-\nu_{k+1}\Big{(}1-\frac{A_{m}}{k}\Big{)}\Big{)}$
(37)

Now we have a preliminary form of the algorithm for computing $\xi_{k}^{(m)}$:
First we compute recursively $\nu_{k}$ via (37) and then evaluate

$\displaystyle\xi_{0}^{(m)}$
$\displaystyle=2Q=\frac{1}{A_{m}}\Big{(}\theta_{0}-\frac{1}{4}(\nu_{0}+\nu_{1})\Big{)},$
(38) $\displaystyle\xi_{k}^{(m)}$
$\displaystyle=\frac{\nu_{k-1}-\nu_{k+1}}{8k},\;\;\;1\leq k\leq N.$ (39)

There are two problems that we need to overcome to make this work. First of
all, when performing backward iteration in (37) we need to start with some
values $\nu_{N}$ and $\nu_{N+1}$ in order to compute $\nu_{k}$ for
$k=N-1,N-2,\dots,1,0$, however we do not know the values of $\nu_{N}$ and
$\nu_{N+1}$. Second, we do not know the values of coefficients $\theta_{k}$ in
the Chebyshev expansion (32), thus we will need to approximate them and
control the resulting error. Let us first consider the former problem.

We define $\\{\hat{\nu}_{k}\\}_{0\leq k<N}$ to be the values obtained by
recursion (37) starting from $\hat{\nu}_{N}=\hat{\nu}_{N+1}=0$. Let us
estimate the difference $y_{k}=\nu_{k}-\hat{\nu}_{k}$ for $0\leq k<N$. Note
that $y_{k}$ satisfy the homogeneous recurrence equation

$y_{k-1}=\frac{1}{1+\frac{A_{m}}{k}}\Big{(}-2y_{k}-y_{k+1}\Big{(}1-\frac{A_{m}}{k}\Big{)}\Big{)},\;\;\;0\leq
k<N$ (40)

started with values $y_{N}=\nu_{N}$ and $y_{N+1}=\nu_{N+1}$. We denote ${\bf
u}_{k}=[y_{k-1},y_{k}]^{T}\in{{\mathbb{R}}^{2}}$ and

$Q_{k}:=\begin{bmatrix}-\frac{2}{1+A_{m}/k}&-\frac{1-A_{m}/k}{1+A_{m}/k}\\\
1&0\end{bmatrix}.$

Then recursion (40) can be rewritten in vector-matrix form as follows

${\bf u}_{k}=Q_{k}{\bf u}_{k+1}$

and we find

${\bf u}_{k}=Q_{k}Q_{k+1}\dots Q_{N}{\bf u}_{N+1}.$ (41)

Let us also define $R_{k}=V^{-1}Q_{k}V$, where

$V:=\begin{bmatrix}-2&1\\\ 2&0\end{bmatrix}.$

One can check that

$R_{k}=\begin{bmatrix}-1&\frac{1}{2}\\\
0&-\frac{1-A_{m}/k}{1+A_{m}/k}\end{bmatrix},$

so that $R_{k}$ belongs to class ${\mathcal{D}}$ of matrices of the form

$\begin{bmatrix}-1&\frac{1}{2}\\\ 0&d_{2,2}\end{bmatrix}$

with $|d_{2,2}|<1$. It is easy to show by induction that a product of $l$
matrices from class ${\mathcal{D}}$ has the form

$\begin{bmatrix}(-1)^{l}&d_{1,2}\\\ 0&d_{2,2}\end{bmatrix}$

where $|d_{1,2}|<l/2$ and $|d_{2,2}|<1$, thus the $\lVert\cdot\rVert_{\infty}$
norm of this matrix is less than $1+l$. From (41) we conclude

${\bf u}_{k}=VR_{k}R_{k+1}\dots R_{N}V^{-1}{\bf u}_{N+1}$

so that $\lVert{\bf u}_{k}\rVert_{\infty}<6(2+N-k)\lVert{\bf
u}_{N+1}\rVert_{\infty}$ for $0\leq k\leq N$ (since $\lVert
V\rVert_{\infty}=3$ and $\lVert V^{-1}\rVert_{\infty}=2$). Since $\nu_{k}$ are
coefficients in the Chebyshev expansion (33) of the function $F^{\prime}(z)$
that is analytic in ${\mathbb{C}}\setminus[-1,\infty)$, we know that
$\nu_{k}=O(\rho^{-k})$ for any $\rho\in(0,3+\sqrt{8})$. Thus $\lVert{\bf
u}_{N+1}\rVert_{\infty}=O(\rho^{-N})$ and for every $0\leq k\leq N$ we have
$\lVert{\bf u}_{k}\rVert_{\infty}=O(N\rho^{-N})$ for any
$\rho\in(0,3+\sqrt{8})$. Since $3^{8/5}<(3+\sqrt{8})$, for
$\rho\in(3^{8/5},3+\sqrt{8})$ we have $N\rho^{-N}=O(3^{-8N/5})$ and we
conclude that $y_{k}=\nu_{k}-\hat{\nu}_{k}=O(3^{-8N/5})$.

Thus we have shown that the backward recursion (37) is stable and if we start
it from $\nu_{N}=\nu_{N+1}=0$ we will contribute an error of size
$O(3^{-8N/5})$.

Now we need to address the second issue, namely, that we do not know the
coefficients $\theta_{k}$ in the Chebyshev expansion (32). However, we can
approximate $f$ by its truncated Chebyshev series, which will give us an
approximation

$(1-z)^{b_{m}-1}f(z)\approx
g(z):=(1-z)^{b_{m}-1}\Big{[}\frac{1}{2}\xi^{(m-1)}_{0}+\sum\limits_{k=1}^{N}\xi^{(m-1)}_{k}T_{k}(4z-1)\Big{]},$

and the difference $(1-z)^{b_{m}-1}f(z)-g(z)$ will be bounded by
$O(\rho^{-N})$ for any $\rho\in(0,3+\sqrt{8})$, uniformly on $0\leq z\leq 1/2$
(this follows from the fact that the coefficients $\xi^{(m-1)}_{k}$ decay at
rate $O(\rho^{-k})$ as $k\to+\infty$). The numbers $\eta_{k}$ computed in (26)
are precisely the Chebyshev coefficients of $g(z)$ on $0\leq z\leq 1/2$. Thus
$|\eta_{k}-\theta_{k}|=O(\rho^{-N})$ for $0\leq k\leq N$ and any
$\rho\in(0,3+\sqrt{8})$. Since we have already demonstrated that the backward
iteration (37) is stable, we see that when we replace $\theta_{k}$ by
$\eta_{k}$ we contribute an overall error of the size $O(3^{-8N/5})$.
$\sqcap\kern-8.0pt\hbox{$\sqcup$}$

##  5 Numerical examples

In this section we present several examples that will illustrate the
performance of the two algorithms for computing the generalized beta function

$\mathrm{B}\Big{(}\begin{matrix}a_{1},\dots,a_{n}\\\
b_{1},\dots,b_{n}\end{matrix}\Big{)}.$

We chose to compute this particular function because it is the normalizing
constant for the ordered beta distribution in (1) and it is also used to
compute various moments of the generalized beta distribution as in (10). These
two computations were instrumental in [3] and it was the question of
efficiently computing the values of the generalized beta function that
motivated this research. However, as we explained in section 4, the
computational complexity of computing the generalized incomplete beta function
is the same and the algorithms are very similar.

For our first example we set $n=3$ and

$[a_{1},a_{2},a_{3}]=[0.8,\;0.3,\;1.5],\;\;\;[b_{1},b_{2},b_{3}]=[0.4,\;1.7,\;0.8].$

We compute the value of

$\mathrm{B}\Big{(}\begin{matrix}a_{1},a_{2},a_{3}\\\
b_{1},b_{2},b_{3}\end{matrix}\Big{)}=0.4868940470437834231542713481277\dots$
(42)

using a Taylor series expansion approach described in section 4.1. The code
was written in Fortran and we used a multiprecision module MPFUN (written by
David Bailey [2]). The Taylor series were truncated at $N=500$ and we did all
computations with precision of 500 digits and then repeated it with the
precision of 1000 digits – this helped us to ensure that all digits displayed
in (42) are correct.

Next we take the value in (42) as the “exact” value and compute approximations
using Taylor’s method and Chebyshev method with with different values of $N$.
Here we have implemented the code in Matlab in standard double precision. We
define ${\mathcal{E}}^{T}(N)$ to be (the absolute value of) the difference
between the exact value in (42) and the value computed using the algorithm
presented in section 4.1 with $N$ terms in the Taylor series expansion.
Similarly we define ${\mathcal{E}}^{C}(N)$ to be (the absolute value of) the
difference between the exact value in (42) and the value computed using the
algorithm presented in section 4.2 with $N$ terms in the Chebyshev series
expansion. On Figure 1a we present the plots of the two errors versus $N$. We
see that both errors converge to zero exponentially fast and that the
Chebyshev approximation converges faster. This confirms our theoretical
results in section 4.

(a)

(b)

Figure 1: The errors of the approximation
${\mathcal{E}}^{{\textnormal{T}}}(N)$ (solid line) and
${\mathcal{E}}^{{\textnormal{C}}}(N)$ (line with circles) versus $N$ on the
$x$-axis.

For our second example we consider parameters

$[a_{1},a_{2},a_{3}]=[50.8,\;0.3,\;1.5],\;\;\;[b_{1},b_{2},b_{3}]=[0.4,\;1.7,\;0.8].$

Note that we have changed the value of $a_{1}$ to a relatively large value of
$50.8$, while keeping all other parameters the same as in the first example.
Again we compute the “exact” value

$\mathrm{B}\Big{(}\begin{matrix}a_{1},a_{2},a_{3}\\\
b_{1},b_{2},b_{3}\end{matrix}\Big{)}=10^{-6}\times
9.9752436394601281551585749018468\dots$

using Fortran code with a multiprecision module. On figure 1b we show the
convergence rate of Taylor series and Chebyshev approximations. We see that
the Chebyshev method still converges exponentially fast and is stable, whereas
the method based on Taylor expansion struggles in this example: the error
first increases to a large value of around $10^{6}$ and only then decreases.
Even with large value of $N=200$ Taylor series method gives us a value of
$10^{-6}\times 9.975628\dots$, thus capturing only four correct digits. This
confirms our observation at the end of section 4.1, that the Taylor series
algorithm is not appropriate in the situation when some of the parameters
$a_{i}$, $b_{i}$ are large, as there is loss of precision arising due to
subtraction of large numbers in (23). In this case one should use either the
Chebyshev method of Taylor series method with high precision.

Finally we consider an example with large number of parameters. Here we
consider an example with 100 parameters and we set

$a_{i}=\frac{2i-1}{200},\;\;\;b_{i}=1-a(i),$

for $1\leq i\leq 100$. Again we compute the “exact” value in Fortran using
precision of 500 digits and 500 terms of Taylor series:

$\mathrm{B}\Big{(}\begin{matrix}a_{1},a_{2},\dots,a_{100}\\\
b_{1},b_{2},\dots,b_{100}\end{matrix}\Big{)}=10^{-33}\times
4.2217553528914884124401921234246\dots$

Then we try to compute the same value using our Matlab code. The method based
on Taylor expansion with $N=50$ terms gives us

$10^{-33}\times 4.221755352891131\dots$

and the method based on Chebyshev expansion with $N=20$ terms gives us a
comparably accurate result

$10^{-33}\times 4.221755352891663\dots.$

We see that both methods work well in this case.

We provide both Matlab programs (Taylor and Chebyshev methods) for computing
the generalized beta function
$\mathrm{B}\Big{(}\begin{matrix}a_{1},a_{2},\dots,a_{n}\\\
b_{1},b_{2},\dots,b_{n}\end{matrix}\Big{)}$ on
https://kuznetsov.mathstats.yorku.ca/code/.

## Acknowledgements

Research of A.K. was supported by the Natural Sciences and Engineering
Research Council of Canada.

## References

  * [1] G. Achtari.  Applications of Machine Learning in Revenue Management and Routing.  PhD thesis, Queen’s University, Canada, ISBN 979-8-209-57148-3, 2018
  * [2] D. Bailey.  High-Precision Software Directory.  https://www.davidhbailey.com/dhbsoftware/.
  * [3] T. Levina, Y. Levin, J. McGill, and M. Nediak.  Strategic consumer cooperation in a name-your-own-price channel.  Production and Operations Management, 24(12):1883–1900, 2015.  https://doi.org/10.1111/poms.12412.
  * [4] J. C. Mason and D. C. Handscomb.  Chebyshev Polynomials.  Chapman and Hall/CRC, 2002.  https://doi.org/10.1201/9781420036114.
  * [5] D. Zwillinger and A. Jeffrey.  Table of integrals, series and products.  Academic Press, 7th edition, 2007.