# Unknown Script: Impact of Script on Cross-Lingual Transfer

Wondimagegnhue Tsegaye Tufa, Ilia Markov, Piek Vossen
Vrije Universiteit Amsterdam
De Boelelaan 1105, 1081 HV Amsterdam, The Netherlands
{w.t.tufa, i.markov<EMAIL_ADDRESS>

###### Abstract

Cross-lingual transfer has become an effective way of transferring knowledge
between languages. In this paper, we explore an often-overlooked aspect in
this domain: the influence of the source language of the base language model
on transfer performance. We conduct a series of experiments to determine the
effect of the script and tokenizer used in the pre-trained model on the
performance of the downstream task. Our findings reveal the importance of the
tokenizer as a stronger factor than the sharing of the script, the language
typology match and the model size.

##  1 Introduction

The dominant Natural Language Processing (NLP) approach now involves cross-
lingual transfer using monolingual and multilingual pre-trained language
models, especially for lesser-resourced languages. In line with this trend,
numerous monolingual models have been released for various languages,
including but not limited to English, Spanish, and Arabic Devlin et al.
(2019); Cañete et al. (2023); Antoun et al. (2020). Multilingual versions of
these models are trained for over 100 languages, such as xlm-R Conneau et al.
(2019) and mBert Devlin et al. (2019).

Despite the advancements in language models for high-resourced languages, many
world languages remain excluded from these models. Out of over 6,500 spoken
languages globally, less than 2% are represented in the current models,
effectively leaving the rest unseen and unaccounted for in the current
language processing technology Muller et al. (2021); Hammarström (2016); Joshi
et al. (2020).

Training a monolingual model for each of the world’s 6,500+ spoken languages
is impractical due to substantial data and compute resource requirements.
Several alternative approaches have been proposed, such as language adapters
Pfeiffer et al. (2020) and multilingual models Conneau et al. (2019); Devlin
et al. (2019). However, these methods are either limited to the script that is
part of the pre-trained model in the case of language adapters or extending
existing multilingual models to accommodate new languages with these languages
creates interference between languages Lauscher et al. (2020); Conneau et al.
(2019); Wang et al. (2020).

Zero-shot transfer, where a pre-trained model is fine-tuned on labeled data in
a resource-rich source language and then directly applied to a resource-poor
target language, has also been considered Pires et al. (2019); Conneau et al.
(2019). While this method is appealing, it often yields sub-optimal results
Lauscher et al. (2020); Pfeiffer et al. (2021), although using even a small
amount of labeled data in the target language has been shown to be effective
Lauscher et al. (2020).

The effectiveness of cross-lingual transfer is influenced by several factors,
such as language similarity between the source and the target language, the
target language’s script, the pre-trained model’s size, and the amount of the
pre-training data Muller et al. (2021); Pires et al. (2019); Cao et al.
(2020); Wu and Dredze (2019).

In this paper, we explore which factors determine cross-lingual transfer
performance for a new language with a unique script. We consider a challenging
case where the target language is not part of the pre-trained model, and the
script is also not seen in the pre-training. Our analysis targets three main
factors: language similarity, tokenization methods, and script.

We design an experiment by varying these factors. We evaluate a range of
existing monolingual and multilingual models, specifically choosing those
trained in typologically related or typologically distant from our target
language. Furthermore, we select models trained on various tokenizers,
allowing us to assess how these choices impact model performance for a
language with unique characteristics. In this way, we want to measure the
impact of different factors on transfer to a new language with a unique
script.

Our analysis shows that the strongest factor for transfer is the tokenizer
used, more than sharing of the script, the language typology match and the
model size. Using a Byte Piece Encoder was most effective to leverage the
pretrained model knowledge during fine-tuning on the target language.

##  2 Related Work

Languages not typically represented in language models display considerable
variation in several aspects Joshi et al. (2020). Firstly, the amount of data
for these languages varies greatly. Secondly, many of these languages use
scripts different from the Latin script Muller et al. (2021); Joshi et al.
(2020). Finally, regarding linguistic characteristics, these languages often
have distinct morphological and syntactical properties, especially when
compared to high-resourced Indo-European languages. Recently, cross-lingual
transfer has become an effective method for extending the capabilities of pre-
trained monolingual and multilingual models for various languages. In this
section, we present an overview of works exploring cross-lingual transfer.

#### Multilingual Models and Adapters

Multilingual models enable transfer between high-resource and low-resource
languages Conneau et al. (2019); Devlin et al. (2019). However, they suffer
from the ‘curse of multilingualism’ and interference between languages Conneau
et al. (2019); Wang et al. (2020), where the model’s effectiveness decreases
as the number of languages increases due to the parameter limit of the model.
Language adapters address these challenges by storing language-specific
knowledge of each language in dedicated parameters Pfeiffer et al. (2020).
This increases the capacity of the multilingual model without introducing
interference between languages. These methods are however not directly
applicable to languages that use scripts not covered in the training data of
these models Pfeiffer et al. (2021).

#### Zero-shot and Few-shot Transfer

In zero-shot transfer, a fine-tuned model on resource-rich source language is
directly applied to a resource-poor target language Pires et al. (2019);
Conneau et al. (2019). While this method is appealing, it often yields sub-
optimal results Lauscher et al. (2020); Pfeiffer et al. (2021). Alternatively,
using even a small amount of labeled data in the target language (few-shot
transfer) has shown to be more effective Lauscher et al. (2020). It remains
unclear what factors determine this effect and to what extent.

#### Factor analysis

The success of cross-lingual transfer is impacted by various factors Muller et
al. (2021); Pires et al. (2019); Cao et al. (2020); Wu and Dredze (2019).
Muller et al. (2021) demonstrated that the performance of transfer can
significantly differ based on factors such as the script of the language, the
amount of available data, and the relationship between the language used for
pre-training and the target language. However, the literature concerning the
effect of script and tokenizer is mixed. For example, Muller et al. (2021)
identifies script as the most crucial factor affecting transfer performance
and shows that transliterating the script to Latin enhances the performance of
cross-lingual transfer. Contrary to this Artetxe et al. (2019) and Karthikeyan
et al. (2020) show that script and lexical overlap are relevant and that large
monolinqual models learn semantic abstractions that are generalizable to other
languages. A similar mixed result has been reported on the effect of
tokenizers in cross-lingual transfer. Rust et al. (2021) show that tokenizers
are a crucial factor in the success of cross-lingual transfer for multilingual
models, while Wu et al. (2022) report it as less important. The analysis of
Muller et al. (2021) covers multilingual models while Wu et al. (2022) focuses
only on an English model.

While similar to the approach of Muller et al. (2021), our study takes a
different direction. Instead of analyzing the performance of multiple
languages with a single multilingual model, we focus on one language that is
unique in its script and not covered by existing models. We select Amharic as
our target language. Amharic is a Semitic language and is classified as
morphologically complex. It has a unique script, distinct from the Latin and
Arabic alphabet, with no shared characters. Additionally, it is categorized as
a Class 2 under-resourced language, indicating a significant lack of data and
tools for language processing Joshi et al. (2020); Adelani et al. (2021). We
evaluate a range of existing monolingual and multilingual models, specifically
choosing those that are trained in languages either closely related or quite
distinct from Amharic. Furthermore, we select models trained on various
tokenizers, allowing us to assess how these choices impact model performance
for a language with a unique script. In this way, we want to measure the
impact of different factors on transfer performance to this language.

We focus on two main research questions:

  1. 1.

To what extent does the script of a source language influence cross-lingual
transfer to a new language in monolingual and multilingual models?

  2. 2.

To what extent does the tokenizer influence cross-lingual transfer to a new
language in monolingual and multilingual models?

##  3 Methodology

We experimented with a few shot settings where we fine-tuned pre-trained
models on target language tasks where the target language and script were not
part of the pre-training. Our few-shot settings follow a standard setting
where we take an existing base model, fine-tune it, and test it on the
original target language.

#### Language

For our test language, we have selected Amharic. As the source languages for
our study, we have chosen English, which is typologically distant from
Amharic, and Arabic, which is typologically related. But both of these
languages do not share their script with Amharic. We use the original Fidel
and Romanized scripts for our fine-tuning and evaluation.

###  3.1 Task and model

Table 1 shows a summary of pre-trained models, the tokenizer, and the language
used in the pre-trained models. We select English and Arabic monolingual and
various multilingual models.

Model | Tokenizer | Model-Type | Language | Model Size
---|---|---|---|---
BERT Devlin et al. (2019) | WordPiece | Monolingual | English | 110M
RoBERTa Liu et al. (2021) | BPE | Monolingual | English | 125M
ALBERT Devlin et al. (2019) | SentencePiece | Monolingual | English | 12M
Bert-base-arabic Antoun et al. (2020) | WordPiece | Monolingual | Arabic | 110M
CANINE-c Clark et al. (2022) | Character | Multilingual | Multiple | 121M
m-BERT Devlin et al. (2019) | WordPiece | Multilingual | Multiple | 110M
Table 1: Overview of models we use, their tokenizers, model types, languages
and model size

#### Task and Datasets

We conducted our experiment on two tasks: Named Entity Recognition (NER) and
part-of-speech (POS) tagging. For the NER task, we use MasakhaNER Adelani et
al. (2021), and for POS, we use the Amharic-POS dataset from Gezmu et al.
(2021). The MasakhaNER dataset is a publicly accessible resource for NER tasks
in ten African languages. This dataset has four types of entity labels: PER
(Person), ORG (Organization), LOC (Location), and DATE (Date). It has 1750 for
training, 250 for validation, and 500 for the test set. The POS dataset
contains 218K sentences with 18 POS tagsets. We sample 2.5K examples, with
1750 sentences for training, 250 for validation, and 500 for testing.

#### Tokenizers

Our model selection encompasses the most widely used tokenizers: SentencePiece
Kudo and Richardson (2018), BPE Sennrich et al. (2016), character-based
tokenizer Clark et al. (2022) and WordPiece Schuster and Nakajima (2012).

#### Fine-tunning

We use a standard fine-tuning setup for all models on two tasks. Our few-shot
setting differs from a standard setting in that all base language models do
not see our target language or the script during pre-training. We think this
test evaluates two capabilities. First, we evaluate whether the models we are
testing enable cross-lingual transfer to a new language and script without
pre-training on these specific linguistic elements. This will help us answer
whether a base language model matters when we target a new language and unseen
script. Second, our evaluation tests how well these models adapt to a new
language and script with limited fine-tuning data.

Model | Fidel-NER | Latin-NER | Fidel-POS | Latin-POS
---|---|---|---|---
| F1 | P | R | F1 | P | R | F1 | P | R | F1 | P | R
Roberta-base | 0.57 | 0.59 | 0.55 | 0.55 | 0.53 | 0.56 | 0.82 | 0.79 | 0.85 | 0.86 | 0.85 | 0.86
CANINE-c | 0.12 | 0.13 | 0.1 | 0.09 | 0.12 | 0.07 | 0.21 | 0.33 | 0.15 | 0.19 | 0.31 | 0.13
Bert-base | - | - | - | 0.57 | 0.6 | 0.55 | 0.24 | 0.49 | 0.16 | 0.82 | 0.81 | 0.84
M-Bert | - | - | - | 0.59 | 0.62 | 0.57 | 0.23 | 0.48 | 0.15 | 0.84 | 0.83 | 0.85
Bert-arabic | - | - | - | 0.59 | 0.6 | 0.58 | 0.24 | 0.45 | 0.16 | 0.81 | 0.8 | 0.82
Albert-base | - | - | - | 0.5 | 0.5 | 0.49 | 0.24 | 0.47 | 0.16 | 0.76 | 0.76 | 0.79
Table 2: Results of few-shot experiments where we fine-tune different models
in NER and POS tasks with Fidel and romanized Latin script representations.
The empty cell shows that we don’t observe a decrease in the loss. We fine-
tune all our models for 25 epochs. The f1-score is averaged over five runs and
has a standard deviation between 0.003 and 0.009. The highest f1-score for
each script is highlighted in bold

##  4 Results and Analysis

Table 2 shows the result for the few-shot setting on NER and POS fine-tuned in
the original Amharic script (Fidel) and its Romanized version. Roberta-base
model stands out, showing robust performance across both scripts, with a
marginal preference for the Fidel script for the NER task and the Latin script
for the POS task. Models such as Bert-base, m-Bert, bert-Arabic, and Albert-
base fail to recognize entities in the Fidel script entirely, reflected in
zero F1 scores, but show some capability with the Latin script. This pattern
persists across both tasks, although POS task has less variation.

#### Effect of the script

The result difference between the Fidel and Romanized versions highlights the
script’s effect on model performance. Models pre-trained on data predominantly
in the Latin script struggle significantly with tasks in the Fidel script, as
evidenced by the drastically lower F1 scores for most models on the Fidel
script compared to the Latin script. This suggests a strong bias towards the
script used during the training phase, with models performing better on
scripts they have seen before. This is in line with the result reported by
Purkayastha et al. (2023); Muller et al. (2021), which shows the romanization
of unknown script boosts transfer performance in multilinqaul model. But we
also observe this in all of the monolingual models we test. The Roberta model
seems to be an exception, showing a good performance before Romanization,
though Romanization also boosts its performance.

#### Language relatedness

Here, we observe no difference between the two models (m-bert and bert-
Arabic), pre-trained in Arabic, and the four English-based models. English
bert-base and Arabic bert can be directly compared since they are trained with
the same training objective, with the same model size and tokenizer, but with
different languages. Still, they show mixed results, with a slight improvement
for the NER task and decreased performance for the POS task.

#### Model-size and Tokenizers

To explain the possible performance difference between Roberta and the other
models, two possible factors can be considered: model size and tokenizer. When
comparing the models’ size in terms of the number of parameters, roberta-base
is the largest with 125 million parameters. It is followed closely by
CANINE-c, which has 121 million parameters. The bert-base, m-Bert, and bert-
arabic models are of equal size, each with 110 million parameters. The Albert
base model is significantly smaller and has only 12 million parameters.
However, the performance does not consistently follow the model size. The
other possible variable is the tokenizer. Roberta is trained with BPE over raw
bytes instead of Unicode characters. The results show that the BPE
representation enables the model to leverage knowledge to that target language
that benefits the downstream tasks even if the script is not included in the
model’s vocabulary as such.

##  5 Conclusion

In this paper, we explore cross-lingual transfer in less explored but
challenging settings where the target language is not seen during pre-training
and has a unique non-Latin script. Our analysis shows considerable differences
in cross-lingual transfer performance among various models, possibly
attributable to two key factors: the size of the pre-trained model and the
specific tokenizer used during pre-training. The model’s size could impact its
ability to capture and generalize across multiple languages, including
language distant from the pre-trained model’s language with new scripts. In
light of recent works that show the importance of tokenizers for
underresourced languages, we found that the difference in tokenizers causes a
proportional difference in performance. Meanwhile, the tokenizer –responsible
for breaking down text into pieces the model can understand—- also plays a
crucial role in the model’s performance across different languages.

##  6 Limitations

In our analysis, we try to control or consider various factors that could
influence the comparative performance of different models. However, it must be
noted that residual differences in model parameters and the extent of pre-
training data may have contributed to the observed disparities in results.
Furthermore, our study did not involve training a model from scratch with
constant architecture, parameters, data, and domain while varying only the
tokenizer or the model size. This limitation precludes a definitive conclusion
about the isolated effect of the tokenizer on model performance. Hence, we
identify the need for additional research where these elements are carefully
controlled. Such investigations would enable a more robust understanding of
the tokenizer’s role and interaction with other model characteristics in
cross-lingual transfer learning.

## References

  * Adelani et al. (2021) David Ifeoluwa Adelani, Jade Z. Abbott, Graham Neubig, Daniel D’souza, Julia Kreutzer, Constantine Lignos, Chester Palen-Michel, Happy Buzaaba, Shruti Rijhwani, Sebastian Ruder, Stephen Mayhew, Israel Abebe Azime, Shamsuddeen Hassan Muhammad, Chris Chinenye Emezue, Joyce Nakatumba-Nabende, Perez Ogayo, Aremu Anuoluwapo, Catherine Gitau, Derguene Mbaye, Jesujoba O. Alabi, Seid Muhie Yimam, Tajuddeen Gwadabe, Ignatius Ezeani, Rubungo Andre Niyongabo, Jonathan Mukiibi, Verrah Otiende, Iroro Orife, Davis David, Samba Ngom, Tosin P. Adewumi, Paul Rayson, Mofetoluwa Adeyemi, Gerald Muriuki, Emmanuel Anebi, Chiamaka Chukwuneke, Nkiruka Odu, Eric Peter Wairagala, Samuel Oyerinde, Clemencia Siro, Tobius Saul Bateesa, Temilola Oloyede, Yvonne Wambui, Victor Akinode, Deborah Nabagereka, Maurice Katusiime, Ayodele Awokoya, Mouhamadane Mboup, Dibora Gebreyohannes, Henok Tilaye, Kelechi Nwaike, Degaga Wolde, Abdoulaye Faye, Blessing K. Sibanda, Orevaoghene Ahia, Bonaventure F. P. Dossou, Kelechi Ogueji, Thierno Ibrahima Diop, Abdoulaye Diallo, Adewale Akinfaderin, Tendai Marengereke, and Salomey Osei. 2021.  Masakhaner: Named entity recognition for african languages.  In _2nd AfricaNLP Workshop Proceedings, AfricaNLP@EACL 2021, Virtual Event, April 19, 2021_.
  * Antoun et al. (2020) Wissam Antoun, Fady Baly, and Hazem Hajj. 2020.  AraBERT: Transformer-based model for Arabic language understanding.  In _Proceedings of the 4th Workshop on Open-Source Arabic Corpora and Processing Tools, with a Shared Task on Offensive Language Detection_ , pages 9–15, Marseille, France. European Language Resource Association.
  * Artetxe et al. (2019) Mikel Artetxe, Sebastian Ruder, and Dani Yogatama. 2019.  On the cross-lingual transferability of monolingual representations.  _arXiv preprint arXiv:1910.11856_.
  * Cao et al. (2020) Steven Cao, Nikita Kitaev, and Dan Klein. 2020.  Multilingual alignment of contextual word representations.  _CoRR_ , abs/2002.03518.
  * Cañete et al. (2023) José Cañete, Gabriel Chaperon, Rodrigo Fuentes, Jou-Hui Ho, Hojin Kang, and Jorge Pérez. 2023.  Spanish pre-trained bert model and evaluation data.
  * Clark et al. (2022) Jonathan H. Clark, Dan Garrette, Iulia Turc, and John Wieting. 2022.  <scp>canine</scp>: Pre-training an efficient tokenization-free encoder for language representation.  _Transactions of the Association for Computational Linguistics_ , 10:73–91.
  * Conneau et al. (2019) Alexis Conneau, Kartikay Khandelwal, Naman Goyal, Vishrav Chaudhary, Guillaume Wenzek, Francisco Guzmán, Edouard Grave, Myle Ott, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Unsupervised cross-lingual representation learning at scale.  _arXiv preprint arXiv:1911.02116_.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 4171–4186, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Gezmu et al. (2021) Andargachew Mekonnen Gezmu, Binyam Ephrem Seyoum, Michael Gasser, and Andreas Nürnberger. 2021.  Contemporary amharic corpus: Automatically morpho-syntactically tagged amharic corpus.  _CoRR_ , abs/2106.07241.
  * Hammarström (2016) Harald Hammarström. 2016.  Linguistic diversity and language evolution.  _Journal of Language Evolution_ , 1(1):19–29.
  * Joshi et al. (2020) Pratik Joshi, Sebastin Santy, Amar Budhiraja, Kalika Bali, and Monojit Choudhury. 2020.  The state and fate of linguistic diversity and inclusion in the NLP world.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pages 6282–6293, Online. Association for Computational Linguistics.
  * Karthikeyan et al. (2020) K Karthikeyan, Zihan Wang, Stephen Mayhew, and Dan Roth. 2020.  Cross-lingual ability of multilingual bert: An empirical study.
  * Kudo and Richardson (2018) Taku Kudo and John Richardson. 2018.  SentencePiece: A simple and language independent subword tokenizer and detokenizer for neural text processing.  In _Proceedings of the 2018 Conference on Empirical Methods in Natural Language Processing: System Demonstrations_ , pages 66–71, Brussels, Belgium. Association for Computational Linguistics.
  * Lauscher et al. (2020) Anne Lauscher, Vinit Ravishankar, Ivan Vulić, and Goran Glavaš. 2020.  From zero to hero: On the limitations of zero-shot language transfer with multilingual Transformers.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 4483–4499, Online. Association for Computational Linguistics.
  * Liu et al. (2021) Zhuang Liu, Wayne Lin, Ya Shi, and Jun Zhao. 2021.  A robustly optimized bert pre-training approach with post-training.  In _China National Conference on Chinese Computational Linguistics_ , pages 471–484. Springer.
  * Muller et al. (2021) Benjamin Muller, Antonios Anastasopoulos, Benoît Sagot, and Djamé Seddah. 2021.  When being unseen from mBERT is just the beginning: Handling new languages with multilingual language models.  In _Proceedings of the 2021 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , pages 448–462, Online. Association for Computational Linguistics.
  * Pfeiffer et al. (2020) Jonas Pfeiffer, Ivan Vulić, Iryna Gurevych, and Sebastian Ruder. 2020.  MAD-X: An Adapter-Based Framework for Multi-Task Cross-Lingual Transfer.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 7654–7673, Online. Association for Computational Linguistics.
  * Pfeiffer et al. (2021) Jonas Pfeiffer, Ivan Vulić, Iryna Gurevych, and Sebastian Ruder. 2021.  UNKs everywhere: Adapting multilingual language models to new scripts.  In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing_ , pages 10186–10203, Online and Punta Cana, Dominican Republic. Association for Computational Linguistics.
  * Pires et al. (2019) Telmo Pires, Eva Schlinger, and Dan Garrette. 2019.  How multilingual is multilingual BERT? In _Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics_ , pages 4996–5001, Florence, Italy. Association for Computational Linguistics.
  * Purkayastha et al. (2023) Sukannya Purkayastha, Sebastian Ruder, Jonas Pfeiffer, Iryna Gurevych, and Ivan Vulić. 2023.  Romanization-based large-scale adaptation of multilingual language models.
  * Rust et al. (2021) Phillip Rust, Jonas Pfeiffer, Ivan Vulić, Sebastian Ruder, and Iryna Gurevych. 2021.  How good is your tokenizer? on the monolingual performance of multilingual language models.  In _Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , pages 3118–3135, Online. Association for Computational Linguistics.
  * Schuster and Nakajima (2012) Mike Schuster and Kaisuke Nakajima. 2012.  Japanese and korean voice search.  In _2012 IEEE international conference on acoustics, speech and signal processing (ICASSP)_ , pages 5149–5152. IEEE.
  * Sennrich et al. (2016) Rico Sennrich, Barry Haddow, and Alexandra Birch. 2016.  Neural machine translation of rare words with subword units.  In _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 1715–1725, Berlin, Germany. Association for Computational Linguistics.
  * Wang et al. (2020) Zirui Wang, Zachary C Lipton, and Yulia Tsvetkov. 2020.  On negative interference in multilingual models: Findings and a meta-learning treatment.  _arXiv preprint arXiv:2010.03017_.
  * Wu and Dredze (2019) Shijie Wu and Mark Dredze. 2019.  Beto, bentz, becas: The surprising cross-lingual effectiveness of BERT.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , pages 833–844, Hong Kong, China. Association for Computational Linguistics.
  * Wu et al. (2022) Zhengxuan Wu, Isabel Papadimitriou, and Alex Tamkin. 2022.  Oolong: Investigating what makes crosslingual transfer hard with controlled studies.


# The Dynamics of Interacting Swarms 111NRL/MR/6790–18-9782

Carl Kolon  United State Naval Academy
Department of Mathematics
Chauvenet Hall
572C Holloway Road
Annapolis, MD 21402-5002  Ira B. Schwartz  US Naval Research Laboratory
Plasma Physics Division
Nonlinear Systems Dynamics Section
Code 6792
Washington, DC 20375
Contact<EMAIL_ADDRESS>
###### Abstract

Swarms are self-organized dynamical coupled agents which evolve from simple
rules of communication. They are ubiquitous in nature, and becoming more
prominent in defense applications. Here we report on a preliminary study of
swarm collisions for a swarm model in which each agent is self-propelling but
globally communicates with other agents. We generalize previous models by
investigating the interacting dynamics when delay is introduced to the
communicating agents. One of our major findings is that interacting swarms are
far less likely to flock cohesively if they are coupled with delay. In
addition, parameter ranges based on coupling strength, incidence angle of
collision, and delay change dramatically for other swarm interactions which
result in flocking, milling, and scattering.

###### Contents

  1. 1 Introduction
  2. 2 The Basic Swarm Model
  3. 3 Delay-Coupled Swarms
  4. 4 Conclusion

##  1 Introduction

The emergence of complex dynamical behaviors from simple local interaction
rules between pairs of agents in a group is a widespread phenomenon over a
range of application domains, and results in coherent behavior called a swarm.
Many striking examples can be found in biological systems, from the
microscopic (ex., aggregates of bacterial cells or the collective motion of
skin cells in wound healing) [4, 32, 17] to large-scale aggregates of fish,
birds, and even humans [40, 14, 18]. These systems are particularly
interesting to the robotics community because they allow simple individual
agents to achieve complex tasks in ways that are scalable, extensible, and
robust to failures of individual agents. In addition, these aggregate
behaviors are able to form and persist in spite of complicating factors such
as communication delay and restrictions on the number of neighbors each agent
is able to interact with, heterogeneity in agent dynamics, and environmental
noise. These factors define how swarm behaviors are effected by changes due to
internal and external actors, allow the swarm to change its patterns and
function [36, 35, 37, 15], as well as break apart and reform [29].

A number of studies show that even with simple interaction protocols, swarms
of agents are able to converge to organized, coherent behaviors. Existing
literature on the subject provides a wide selection of both agent-based [14,
18, 43, 40] and continuum models [10, 39, 32]. One of the earliest agent-based
models of swarming is Reynolds’s boids [33], which simulates the motion of a
group of flocking birds. The boids follow three simple rules: collision
avoidance, alignment with neighbors, and attraction to neighbors. Since the
publication of Reynolds’s paper, many models based on “zones” of attraction,
repulsion, and/or alignment have been used as a means of realistically
modeling swarming behaviors [26, 38, 44]. Systematic numerical studies of
discrete flocking based on alignment with nearest neighbors we carried out by
Vicsek et al. [42]. Stochastic interactions between agents are modeled in
[31]. In recent years, improved computer vision algorithms have allowed
researchers to record and analyze the motions of individual agents in
biological flocks, and formulating more accurate, empirical models for
collective motion strategies of flocking species including birds and fish [2,
16, 5].

Despite the multitude of available models, how group motion properties emerge
from individual agent behaviors is still an active area of research. For
example, [45] presents a simulation-based analysis of the different kinds of
motion in a fish-schooling model; the authors map phase transitions between
different aggregate behaviors as a function of group size and maximum number
of neighbors that influence the motion of each “fish”. In [18], the authors
use simulation to study transitions in aggregate motions of prey in response
to a predator attack.

Interaction delay is a ubiquitous problem in both naturally-occurring and
artificial systems, including blood cell production and coordinated flight of
bats [22, 3, 28, 12, 11]. Communication delay can cause emergence of new
collective motion patterns and lead to noise-induced switching between
bistable patterns [23, 24, 20]; this, in turn, can lead to instability in
robotic swarming systems [44, 21]. Understanding the effects of delay and
latency is key to understanding many swarm behaviors and functions in natural,
as well as engineered, systems [37]. As another controlling factor of swarm
function, network communication topology maybe used to improve and generalize
swarm function. For example in [15], it was shown how to take known basic
swarm modes, such as flocking, milling and rotating states, and use network
topology to create new hybrid states that form new functions, such as
shielding important swarm controlling agents.

In addition, many models make the mathematically simple but physically
implausible assumption that swarms are globally coupled (that is, each agent
is influenced by the motion of all other agents in the swarm) [30, 8, 7, 18,
41]. Global coupling is easier to analyze and a reasonable assumption in cases
of high-bandwidth communication, with a sufficiently small number of agents.
In contrast, we are interested in the collective motion patterns that emerge
when global communication cannot be achieved. New behaviors can unexpectedly
emerge when the communication structure of a network is altered, as in [46],
where the stability of solutions for compromise dynamics over an Erdös-Renyi
communication network is considered. However, in our system, we show
robustness of emergent motion patterns to loss of communication links in
presence of delayed coupling [37].

A third effect, which we do not consider here, is agent heterogeneity. Most
existing work assumes that the members of the swarm are identical. However,
many practical applications involve swarms that are composed of agents with
differing dynamical properties from the onset, or that become different over
time due to malfunction or aging. Swarm heterogeneity leads to interesting new
collective dynamics such as spontaneous segregation of the various populations
within the swarm; it also has the potential to erode swarm cohesion. In
biology, for example, it has been shown that sorting behavior of different
cell types during the development of an organism can be achieved simply by
introducing heterogeneity in inter-cell adhesion properties [34, 13]. In
robotic systems, allowing for heterogeneity in dynamical behaviors of swarm
agents gives greater flexibility in system design, and is therefore desirable
not only from a theoretical but also from a practical point of view [15].

In this preliminary work we are interested in swarms that are interacting.
That is, the entire swarm is split into two parts initially, and then the two
components approach each other as flocks. We then are interested in
understanding the results based on physical parameters, such as
attraction/repulsion length scales, delayed communication, and coupling
strength. The model for the swarm we use is based on the the employed in [9],
which describe a mathematically swarm model using the Morse potential.
Recently, the authors in [1] studied the effects of swarm collisions using the
Morse potential as a function of incident angle of the interacting flocks.
However, in any robotic swarm there must be some delay between sensing other
agents, and control actuation. We investigated the result of applying delay to
the coupling terms of a swarm model. As a result, every agent behaves as
though the other agents are at the positions they were at $\tau$ time units
before. We found that the flocking state is unlikely following collisions of
delay-coupled swarms, even for fairly small values of $\tau$. The implication
is that it should be possible to have one swarm modify the intent of another
by either capturing the swarm, or redirecting the mean direction of the flock.

##  2 The Basic Swarm Model

The swarm model we used for this preliminary study follows the model in [9].
It is a system of $N$ agents in $n$-dimensional space with position vectors
$x_{i}\in\mathbb{R}^{n}$, acting under the following equation of motion:

$\ddot{x}_{i}=(\alpha-\beta\lvert\dot{x}_{i}\rvert^{2})\dot{x}_{i}-\frac{\lambda}{N}\sum_{j=1,i\neq
j}^{N}\nabla_{x_{i}}U(x_{i},x_{j}),$ (1)

where $\alpha,\beta,\lambda$ are constants,
$U:\mathbb{R}^{2n}\rightarrow\mathbb{R}$ is a potential function of the two
agents’ position, and $\nabla_{x_{i}}$ represents taking the gradient with
respect to $x_{i}$. Many potential functions can be used, as in [27], but a
useful choice is the scaled Morse potential:

$U(x_{i},x_{j})=C\exp(-\lvert x_{i}-x_{j}\rvert/l)-\exp(-\lvert
x_{i}-x_{j}\rvert).$ (2)

In particular, the parameter $l$ can be used to set a ratio of repulsive and
attractive length scales.

Following [1], we fix the values of $\alpha$, $\beta$, $C$, and $l$. We also
constrain our swarm to the plane. We choose $\alpha=1$, $\beta=5$,
$C=\frac{10}{9}$, $l=\frac{3}{4}$. This means that, at long ranges,
$\nabla_{x_{i}}U(x_{i},x_{j})$ will be attractive, and at short ranges it will
be repulsive (see Fig. 1). Furthermore, the parameters choice places the swarm
in the region which is defined as ”catastrophic” [9], meaning that increasing
the number of agents decreases the space between agents, so that the total
size of the swarm approaches a limit in space.

(a) Morse potential as a function of distance.

(b) Morse potential of a variable point on a plane and a point at the origin.

Figure 1: The Morse potential with $C=\frac{10}{9}$ and $l=\frac{3}{4}$.

When $\nabla_{x_{i}}U(x_{i},x_{j})\equiv 0$ for all $i,j\leq N$ with $i\neq
j$, the dynamics is based solely on self-propulsion, and the resulting
configuration of agents is called a flock [19]. That is, all of the agents
when started in parallel, will translate with the same asymptotic velocity and
direction. In the attraction-repulsion model, a flock undergoes translating
motion at a constant speed $\sqrt{\alpha/\beta}$, and can withstand
sufficiently small perturbations [6]. However, due to translation symmetry,
there will always be one eigenvalue in the linearization about the flocking
state equal to zero, which means the stability is neutral at best. Hence small
perturbations may destabilize the flocking dynamics [25].

Since we are interested in flock collisions, two separate flocks are
initialized and pointed towards each other with some incident angle, $\theta$.
We are interested in the behavior after collision, and how it varies for
varying $\lambda$, $\theta$, and $\tau$.

Figure 2: The incident angle of the collision, $\theta$, which we vary as a
parameter.

Flocking is not the only behavior that swarms display. They also can perform
milling behavior, where the agents rotate around a stationary center of mass.
To differentiate between the two behaviors, we calculate the polarization of
the swarm as the directed sum of velocities. Specifically, for a system of $N$
agents with position vectors $x_{i}$, the polarization $P$ of the flock is
defined as:

$P(x_{1},x_{2},\ldots,x_{N})=\left\lvert\frac{\sum_{i=1}^{N}\dot{x}_{i}}{\sum_{i=1}^{N}\lvert\dot{x}_{i}\rvert}\right\rvert$

Two limiting behaviors can be immediately seen to be reflected in the
polarization. If the velocities are coherent, we expect a value of $P\approx
1$. If the velocities point in opposite directions, we expect $P\approx 0$.
Therefore, polarization is a good measure of whether the system displays
flocking ($P\approx 1$), or milling ($P\approx 0$). Additionally, two
colliding swarms can scatter, in which they do not form a coherent group which
either flocks or mills. In this case, we expect $P$ to be neither close to 0
nor close to 1. See Fig. 3 for examples of the three types of swarm
interactions.

(a) Flocking.

(b) Milling.

(c) Scattering.

Figure 3: Flocking, milling, and scattering behavior in a Morse swarm.

##  3 Delay-Coupled Swarms

As mentioned above, since robotic systems possess latency in both
communication and control actuation, we generalized the model to include a
time delay $\tau$ in the interaction term, creating the following delay
differential equation from equation (1):

$\ddot{x}_{i}(t)=(\alpha-\beta\lvert\dot{x}_{i}(t)\rvert^{2})\dot{x}_{i}(t)-\frac{\lambda}{N}\sum_{j-1,i\neq
j}^{N}\nabla_{x_{i}}U(x_{i}(t),x_{j}(t-\tau)),$ (3)

with $U$ defined in Eq. (2). Since this is a delay differential equation, it
requires an initial history function of position and velocity defined over
$[-\tau,0]$, for which we used translating motion of the two flocks along
linear trajectories. We simulated two flocks of 25 agents, each swarm on a
collision path with the other. We varied the values of $\lambda$, the coupling
strength, and $\theta$, the incident angle of the collision.

As a benchmark, we first simulated the results by performing the experiment
without delay;i.e., $\tau=0$, corresponding to a system with instantaneous
interactions. We simulated 400 collisions and recorded the polarization after
transients were removed, about 100 time units. We generated the plot in Fig. 4
showing how polarization can be used to quantify the resulting dynamics after
the swarms collide.

FlockingScatteringMillingMetastability Figure 4: Polarization as a function of
$\lambda$ and $\theta$. Note the distinct regions of flocking, milling,
scattering, and metastability.

From Fig. 4, we immediately observe several large scale features in parameter
space. For low values of $\lambda$, scattering is likely unless the swarm
motion is almost parallel in its initial state. For intermediate values of
$\lambda$, flocking is likely to be observed, while for high values of
$\lambda$, milling dominates the parameter region. A region of metastability
exists between the flocking and milling states. The black lines are used to
guide the eye to show the boundaries separating the final swarm states.

To see how the milling comes into play for the delay-coupled swarm, we start
the flocks approaching each other at zero incidence angle. Figure 5
illustrates the progression from the initial flocks to the final milling
swarm. It shows the potential of how one might use a defensive swarm to stop
another by converting its flocking state, which is a traveling center of mass,
to a milling state with a center of mass that is stationary.

Figure 5: Time snapshots of swarm capture to a milling state. Delay here is
$\tau=0.1$, and other parameters were chosen in the milling regime of Fig.

In exploring the effects of delay on the parameter basins of swarm final
states, we computed the polarization for three distinct values of delay:
$\tau=$ 0, 0.1, and 0.2. The results of the average polarization plots using
the above procedure are in Fig. 6.

(a) $\tau=0$

(b) $\tau=0.1$

(c) $\tau=0.2$

Figure 6: Polarization as a function of $\lambda$ and $\theta$ for 3 distinct
values of $\tau$.

We can note several things about the introduction of delay in panels of Fig.
6. Reference from the $\tau=0$ case, one see a gradually shrinking of the
flocking region, with an increase in milling parameters. Even at the
relatively small delay of $\tau=0.2$, flocking behavior nearly vanishes.
Instead, scattering behavior becomes more likely for low $\lambda$, and
milling behavior becomes more likely for high $\lambda$.

Figure 7: A plot of polarization as a function of $\lambda$ and $\tau$ for
fixed $\theta$. Note that flocking behavior all but disappears above
$\tau\approx 0.1$.

To check the overall dependence of swarm final state on both coupling strength
and delay, we fixed the incidence angle $\theta$ at $\frac{\pi}{4}$, and
varied $\tau$ finely, to observe polarization as a function of $\tau$ and
$\lambda$. The results are shown in Fig. 7. This confirms that as delay
increases steadily, milling behavior largely replaces flocking behavior.

##  4 Conclusion

The preliminary results of this work reveal certain features about how
interacting swarms upon collision modify their behavior. One major finding is
that the existence of a swarming flocking state depends strongly on the
delayed interaction between agents. If delays can be induced into the
communication network of a flocking swarm, the results show that even small
delay disrupts the motion of the center of mass, thus disabling the flocking
state.

In the case of globally coupled interacting swarms which are equal in number,
the resulting swarm states may be in one of three states: Flocking, Milling,
and Scattered. However, the results from this study depend heavily on the
initial states, such as velocity and incidence angle, in addition to delay
magnitude. Nonetheless, the final states of colliding swarms may result in
modifying the center of mass velocity, change the navigation course of the
swarm from its original path, and capture another swarm. Th potential of using
one swarm to modify another also depends on how it is sensing other agents,
which could be through lidar, vision, or from a third party, such as a mother
ship equipped with radar.

The sensitivity of the flocking state also has some other consequences for
modeling. For example, it shows that self-organizing real flocking biological
systems which may be modeled by swarms may have more complicated underlying
dynamics and adaptive controls than are not included in the deterministic
equations of motion. There may be some other factors not modeled here which
biological flocks may use to compensate for actuation delays, noise, and
communication latency. In particular, since the linear stability of the
flocking state is neutral, there must be other individual agent controls to
make the flocking more robust to perturbations.

The preliminary results here also show that globally coupled swarms with
delayed communication cannot be accurately modeled by an instantaneous model.
Similar to work shown in [37], delays cause new states to emerge that are not
observed in systems with instantaneous communication. Moreover, it is also
possible to create new when delay is included in the coupling terms along with
changes in network communication topology which is not global, as shown in
[15]. In the case of a given Morse potential here, simulations accurately
modeling delay-coupled swarms must take the delay into account to get accurate
results.

Finally, since only the Morse potential was used here to model interacting
swarm dynamics, other communication and coupling schemes need to be studied to
see what kinds of general statements may be made, and which conclusion pertain
only to specific swarm models. The end results of such research should prove
to open new areas of interacting and combating swarms in the future.

## References

  * [1] D. Armbruster, S. Martin, and A. Thatcher, Elastic and inelastic collisions of swarms, Physica D: Nonlinear Phenomena, 344 (2017), pp. 45 – 57\.
  * [2] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, V. Lecomte, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, and V. Zdravkovic, Interaction ruling animal collective behavior depends on topological rather than metric distance: evidence from a field study., Proceedings of the National Academy of Sciences of the United States of America, 105 (2008), pp. 1232–1237.
  * [3] S. Bernard, J. Bélair, and M. C. Mackey, Bifurcations in a white-blood-cell production model, Comptes Rendus Biologies, 327 (2004), pp. 201–210.
  * [4] E. O. Budrene and H. C. Berg, Dynamics of formation of symmetrical patterns by chemotactic bacteria, Nature, 376 (1995), pp. 49–53.
  * [5] D. S. Calovi, U. Lopez, S. Ngo, C. Sire, H. Chaté, and G. Theraulaz, Swarming, schooling, milling: Phase diagram of a data-driven fish school model, New Journal of Physics, 16 (2014).
  * [6] J. Carrillo, Y. Huang, and S. Martin, Nonlinear stability of flock solutions in second-order swarming models, Nonlinear Analysis: Real World Applications, 17 (2014), pp. 332–343.
  * [7] Y. Chen and T. Kolokolnikov, A minimal model of predator-swarm interactions, Journal of the Royal Society, Interface, 11 (2014), p. 20131208.
  * [8] Z. Chen, H. Liao, and T. Chu, Clustering in multi-agent swarms via medium-range interaction, EPL (Europhysics Letters), 96 (2011), p. 40015.
  * [9] M. R. D’Orsogna, Y. L. Chuang, A. L. Bertozzi, and L. S. Chayes, Self-propelled particles with soft-core interactions: Patterns, stability, and collapse, Phys. Rev. Lett., 96 (2006), p. 104302.
  * [10] L. Edelstein-Keshet, D. Grunbaum, and J. Watmough, Do travelling band solutions describe cohesive swarms? An investigation for migratory locusts, Journal of Mathematical Biology, 36 (1998), pp. 515–549.
  * [11] E. Forgoston and I. B. Schwartz, Delay-induced instabilities in self-propelling swarms, Physical Review E, 77 (2008), p. 035203.
  * [12] L. Giuggioli, T. J. McKetterick, and M. Holderied, Delayed Response and Biosonar Perception Explain Movement Coordination in Trawling Bats, PLOS Computational Biology, 11 (2015), p. e1004089.
  * [13] F. Graner, Simulation of the differential adhesion driven rearrangement of biological cells, Physical Review E, 47 (1993), pp. 2128–2154.
  * [14] D. Helbing and P. Molnar, Social force model for pedestrian dynamics, Physical Review E, 51 (1995), pp. 4282–4286.
  * [15] J. Hindes, K. Szwaykowska, and I. B. Schwartz, Hybrid dynamics in delay-coupled swarms with ”mothership” networks, PHYSICAL REVIEW E, 94 (2016), p. 032306.
  * [16] Y. Katz, K. Tunstrom, C. C. Ioannou, C. Huepe, and I. D. Couzin, Inferring the structure and dynamics of interactions in schooling fish, Proceedings of the National Academy of Sciences, 108 (2011), pp. 18720–18725.
  * [17] R. M. Lee, D. H. Kelley, K. N. Nordstrom, N. T. Ouellette, and W. Losert, Quantifying stretching and rearrangement in epithelial sheet migration, New Journal of Physics, 15 (2013).
  * [18] S.-H. Lee, Predator’s attack-induced phase-like transition in prey flock, Physics Letters A, 357 (2006), pp. 270–274.
  * [19] H. Levine, W.-J. Rappel, and I. Cohen, Self-organization in systems of self-propelled particles, Physical Review E, 63 (2000).
  * [20] B. S. Lindley, L. Mier-y-Teran Romero, and I. B. Schwartz, Noise Induced Pattern Switching in Randomly Distributed Delayed Swarms, in 2013 American Control Conference, Ieee, June 2013, pp. 4587–4591.
  * [21] Y. Liu, K. M. Passino, and M. Polycarpou, Stability analysis of one-dimensional asynchronous swarms, IEEE Transactions on Automatic Control, 48 (2003), pp. 1848–1854.
  * [22] A. Martin and S. Ruan, Predator-prey models with delay and prey harvesting, Journal of Mathematical Biology, 43 (2001), pp. 247–267.
  * [23] L. Mier-y-Teran Romero, E. Forgoston, and I. B. Schwartz, Noise, Bifurcations, and Modeling of Interacting Particle Systems, in Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems, Jan. 2011, pp. 3905–3910.
  * [24] , Coherent Pattern Prediction in Swarms of Delay-Coupled Agents, IEEE Transactions on Robotics, 28 (2012), pp. 1034–1044.
  * [25] L. Mier-y Teran-Romero, E. Forgoston, and I. B. Schwartz, Coherent pattern prediction in swarms of delay-coupled agents, IEEE Transactions on Robotics, 28 (2012), pp. 1034–1044.
  * [26] J. M. Miller, A. Kolpas, J. P. J. Neto, and L. F. Rossi, A continuum three-zone model for swarms, Bulletin of mathematical biology, 74 (2012), pp. 536–61.
  * [27] A. Mogilner, L. Edelstein-Keshet, L. Bent, and A. Spiros, Mutual interactions, potentials, and individual distance in a social aggregation, Journal of mathematical biology, 47 (2003), pp. 353–89.
  * [28] N. A. M. Monk, Oscillatory Expression of Hes1, p53, and NF-kB Driven by Transcriptional Time Delays, Current Biology, 13 (2003), pp. 1409–1413.
  * [29] D. S. Morgan and I. B. Schwartz, Dynamic coordinated control laws in multiple agent models, Physics Letters A, 340 (2005), pp. 121–131.
  * [30] S. Motsch and E. Tadmor, A New Model for Self-organized Dynamics and Its Flocking Behavior, Journal of Statistical Physics, 144 (2011), pp. 923–947.
  * [31] C. Nilsen, J. Paige, O. Warner, B. Mayhew, R. Sutley, M. Lam, A. J. Bernoff, and C. M. Topaz, Social aggregation in pea aphids: experiment and random walk modeling, PloS ONE, 8 (2013), p. e83343.
  * [32] A. A. Polezhaev, R. A. Pashkov, A. I. Lobanov, and I. B. Petrov, Spatial patterns formed by chemotactic bacteria Escherichia coli, The International Journal of Developmental Biology, 50 (2006), pp. 309–314.
  * [33] C. W. Reynolds, Flocks, herds and schools: A distributed behavioral model, ACM SIGGRAPH Computer Graphics, 21 (1987), pp. 25–34.
  * [34] M. S. Steinberg, Reconstruction of Tissues by Dissociated Cells, Science, 141 (1963), pp. 401–408.
  * [35] K. Szwaykowska, L. Mier-y Teran Romero, , and I. B. Schwartz, Collective motions of heterogeneous swarms, IEEE TRANSACTIONS ON AUTOMATION SCIENCE AND ENGINEERING, 12 (2015), pp. 810–818.
  * [36] K. Szwaykowska, L. Mier-y Teran-Romero, and I. B. Schwartz, Patterned dynamics of delay-coupled swarms with random communication graphs, in Decision and Control (CDC), 2015 IEEE 54th Annual Conference on, IEEE, 2015, pp. 6496–6501.
  * [37] K. Szwaykowska, I. B. Schwartz, L. M.-y.-T. Romero, C. R. Heckman, D. Mox, and M. A. Hsieh, Collective motion patterns of swarms with delay coupling: Theory and experiment, Physical Review E, 93 (2016), p. 032307.
  * [38] I. Tarras, N. Moussa, M. Mazroui, Y. Boughaleb, and A. Hajjaji, Collective Behavior of Interacting Particles: Radius-Dependent Phase Transition, Modern Physics Letters B, 27 (2013), p. 1350028.
  * [39] C. M. Topaz and A. L. Bertozzi, Swarming Patterns in a Two-Dimensional Kinematic Model for Biological Groups, SIAM Journal on Applied Mathematics, 65 (2004), pp. 152–174.
  * [40] K. r. Tunstrø m, Y. Katz, C. C. Ioannou, C. Huepe, M. J. Lutz, and I. D. Couzin, Collective states, multistability and transitional behavior in schooling fish, PLoS computational biology, 9 (2013), p. e1002915.
  * [41] F. Vecil, P. Lafitte, and J. Rosado Linares, A numerical study of attraction/repulsion collective behavior models: 3D particle analyses and 1D kinetic simulations, Physica D, 260 (2013), pp. 127–144.
  * [42] T. Vicsek, A. Czirok, E. Ben-Jacob, I. Cohen, and O. Shochet, Novel Type of Phase Transition in a System of Self-Driven Particles, Physical Review Letters, 75 (1995), pp. 1226–1229.
  * [43] , Novel type of phase transition in a system of self-driven particles, 2006.
  * [44] C. Virágh, G. Vásárhelyi, N. Tarcai, T. Szörényi, G. Somorjai, T. Nepusz, and T. Vicsek, Flocking algorithm for autonomous flying robots, Bioinspiration & biomimetics, 9 (2014), p. 025012.
  * [45] S. V. Viscido, J. K. Parrish, and D. Grünbaum, The effect of population size and number of influential neighbors on the emergent properties of fish schools, Ecological Modelling, 183 (2005), pp. 347–363.
  * [46] J. H. von Brecht, T. Kolokolnikov, A. L. Bertozzi, and H. Sun, Swarming on Random Graphs, Journal of Statistical Physics, 151 (2013), pp. 150–173.


# Hierarchical Bayesian Regression for Multi-Location Sales Transaction
Forecasting

John Mark Agosta Azure Data
Microsoft
Mountain View, CA, USA  Mario Inchiosa  Azure Data
Microsoft
Mountain View, CA, USA

###### Abstract

The features in many prediction models naturally take the form of a hierarchy.
The lower levels represent individuals or events. These units group naturally
into locations and intervals or other aggregates, often at multiple levels.
Levels of groupings may intersect and join, much as relational database tables
do. Besides representing the structure of the data, predictive features in
hierarchical models can be assigned to their proper levels. Such models lend
themselves to hierarchical Bayes solution methods that “share” results of
inference between groups by generalizing over the case of individual models
for each group versus one model that aggregates all groups into one.

In this paper we show our work-in-progress applying a hierarchical Bayesian
model to forecast purchases throughout the day at store franchises, with
groupings over locations and days of the week. We demonstrate using the stan
package on individual sales transaction data collected over the course of a
year. We show how this solves the dilemma of having limited data and hence
modest accuracy for each day and location, while being able to scale to a
large number of locations with improved accuracy.

##  1 Introduction

The field of hierarchical statistical modeling blossomed once simulation
methods such as Markov Chain Monte Carlo (MCMC) made it feasible to solve
models without analytic solutions (Gelman et al. [2013]). Approximation
methods for multi-level models existed in statistics, but subsequently the
field adopted these simulation methods and their Bayesian versions (Gelman and
Hill [2007]). Similarly in AI these techniques appeared in the form of
Probabilistic Graphical Models (PGMs) (Koller and Friedman [2009]). These
graphical DAGs have spawned a bewildering variety of synonymous terms, going
back to the original “Bayesian Belief Networks,” (Pearl [2009] with current
interest in causal reasoning introducing yet another term: “Causal Graphical
Models” Concurrently hierarchical modelling is well represented in the
Economics literature (Geffner et al. [2022]). The similarity in data
structures leads to some curiosity of cross fertilization with Statistical
Relational Learning (Getoor and Tasker [2007]).

It is our belief that these methods are widely applicable and deserve wider
recognition. Most data do not come in flat tables, but come in relational
tables representing different levels of aggregation. The goal of this paper is
to demonstrate the value of modeling hierarchy explicitly in an actual
application.

###  1.1 “Point of Sale” Demand Forecasting

In this paper we present a sales demand model that applies over the fleet of
fast-food store locations, to be used to inform real-time customer demand.
Fast-food items have a short shelf-life; There is no carry-over from one
prediction period to the next. Each location has control of how many products
it prepares in anticipation of customer demand. Demand is uncertain, and each
store needs to anticipate how many items to set aside in anticipation of the
demand in the current period. There is a trade-off: If an item is not sold
punctually— within less than an hour—it is discarded. When too many are made,
at the end of the period items are discarded. When not enough are made,
customers’ demand is not satisfied, creating an “opportunity cost” of lost
revenue.

###  1.2 The Current Cloud-IT infrastructure

With the current IT infrastructure, demand is calculated centrally by a
separate model for each franchise location that generates an estimate of its
demand for the next period. Estimates are transmitted to the “Point of Sales”
(POS) terminals at each location by means of the distributed “Cloud”
application that manages all data interaction between individual store
locations and the firm.

The demand model makes up a small part of this recently implemented
centralized system. After the system went into operation and a few months of
data had been collected, an initial demand prediction model was built using
proprietary automated time-series software. The predictions either tended
toward the median value of zero, or failed to converge and our data science
group was brought in to formulate a replacement model. A primary design
concern is the computational cost of scaling. The intent is to deploy
forecasts to thousands of stores in near-real-time. Our approach is presented
here, where in Section2 we detail the problem and the available data; then in
Section3 explain the hierarchical model we developed. Section4 presents our
results, where we consider these in the context of the original problem and
future approaches in Section5.3.

##  2 Description of the Problem

Forecasting the demand distribution—–the rate that customers arrive and the
quantity of products they request promises not only to reduce cost due to
minimizing waste, but also to increase revenue from anticipating times of
unmet demand. Optimizing production to meet uncertain demand is characteristic
of many production and inventory problems. When as in this case, there is no
inventory hold-over from one period to the next, this specifically is known as
the “news-vendor” problem, (Porteus [2002]) a reference to a newsboy or girl
who must decide how many papers to buy for the day based on their estimate of
how many they will sell.

Sales transactions data consists of a sequence of events in continuous time—a
_stochastic renewal process._ Each event consists of a random number of items
ordered at a random time. The events appear as a dis-continuous function of
time. A sense of the arrival-events data and how the arrival rate varies over
the course of a day is shown in Figure 1. Such event processes do not
constitute a time-series in the conventional sense, confounding conventional
time-series forecasting methods that inevitably fail when applied to such
data. The best one can hope for is to approximate the underlying arrival rate
as a smoothly varying function as the time-series of interest.

Given the uncertainty in the arrival process even with a known, fixed arrival
rate, there is a finite limit in the savings possible by any continuous
function approximation to arrival events. In the simplest case one can imagine
arrivals described by a Poisson distribution whose variance equals the arrival
rate. In actuality arrivals tend to be “over-dispersed”, making variance even
larger. Large variance implies that any prediction based on the arrival rate
is likely to be far from the actual value. For the same reason conventional
time-series error estimates, such as “MAE” are of the order of the arrival
rate parameter, and are of little use for evaluation.

To emphasize the point, as the prediction interval for perishable inventory
gets smaller, the inherent increase in variability in arrivals leads to
unavoidable losses. To illustrate, if the rate is such that one customer is
expected to arrive in a period there are roughly equal probabilities of no
customers arriving, the one customer actually arriving, or more than one
customer arriving. Thus at any instant, even with the best prediction,
inventory will more likely than not be either under or over the instantaneous
predicted amount.

So given the intrinsic uncertainty of a stochastic renewal process that best
describes demand for items, the best one can hope for is to predict an
instantaneous demand rate. We make this concession, apply a rudimentary
regression model at the daily forecast level and focus on modelling the
variability at the level of store locations and days; hence the attractiveness
of a hierarchical approach. The hierarchical model is more interpretable, and
since the upper level data is more concise by a couple orders of magnitude,
the model can scale to thousands of stores.

Short of being able to model the larger in-store production optimization
problem, we leave it up to the store manager to decide on the actual
production policy as informed by our prediction of the instantaneous rate of
demand.

Figure 1: Sales events as time and item count at one location, for one day.
The red curve is a $2^{nd}$-order fit to the log of the item count. Figure 2:
Binned sales events over 15 minute intervals for the same one location and
day. The red curve is a $2^{nd}$-order fit to the log of the binned item
count.

###  2.1 Forming the Hierarchy

At the lowest level of the hierarchy are the individual sales transaction
events consisting of 1,271,639 transactions at 49 stores over the course of
about half a year, for which we have data from 243 days of data, about 150
which are usable. We group the daily samples by location and day_of_week, and
fit a parameterized arrival rate curve to that location-day to create the
upper level data. At the upper level of the hierarchy, the number of samples
is reduced significantly, to
$|locations|\times|days|\times|local\\_model\\_coefficients|$. There are 8649
“location-day” records, one for each day and location, each described by 3
coefficients, a reduction to 2% of the original data size! The schema for the
upper level data is show in Table3. The relations among the
groupings—location, day-of-week and transaction are analogous to those in a
relational data schema, with location and day_of_week as common keys.

We form two upper-level groupings; by location and by day-of-week. Each
grouping could have co-variates; for instance, location could be influenced by
regional economic conditions, and day_of_week by holidays. These extensions to
the model fit naturally into our formulation and arguably could substantially
improve predictions, but here we focus on the basic hierarchical design, and
inclusion of upper-level co-variates is not considered.

###  2.2 Data Sources and Preparation

We were provided by the client with a snapshot of the aforementioned 1,271,639
POS daily sales transaction records over the past year. Each sales transaction
event records the time and number of items purchased by a customer with one-
minute resolution. The sales record fields are shown in Table 1.

This transaction event data is binned on 15-minute intervals to generate
counts on a fixed time grid reducing the transactions to 357,224 records.
Intervals when no sales occur are preserved and recorded as zero counts.
Binning is necessary to fit an estimate of the instantaneous transaction rate
not biased by the timing of transaction events.

To create the upper hierarchy levels the binned data for each day and location
were summarized by a curve fit by 3 coefficients. Figure 2 is an example of a
curve for one case. These cases make up the coefficients data set used for the
hierarchical model, with fields described in Table 3. We assign records in the
coefficients data set randomly with probability $1/2$ to train and test sets.

Inspection of the data over the year’s interval showed a slight increasing
demand trend plus a substantial interval of missing data for the middle
quarter of the year. Our initial analysis, using the early data for train and
later for test revealed that the data set was not stationary. This approach
confounded any test results on the data split and was abandoned.

Characteristic of hierarchical models, our premise is that the variation
within the two groups, location and day_of_week explain the variation in the
coefficients data. This boxplot, Figure3, shows the variation of the first
coefficient $c0$, between and within days of the week. A similar plot for
stores is shown in Figure4. As one might expect, the finer-grained store
grouping shows less within-group variation and more between-group variation
than among days of the week.

Figure 3: These barplots show the variation in the $c0$ coefficient for each
day of the week. The values for each day varies widely, with some weekly
variation peaking on Saturdays. Figure 4: Sales events as time and item count
at one location, for one day. The red curve is a $2^{nd}$-order fit to the log
of the item count. Table 1: Fields in the sales transaction dataset.

Field name | type | Description
---|---|---
LocationNumber | int | Vendor identifier
SalesDayName | string | Day of the week
DailyMinutesOpen | int | Daily operating minutes
DateTimePlaced | datetime | Transaction timestamp
SalesAsMinutes | double | Minutes since opening
Quantity | int | Number of items purchased

Table 2: Fields in the hierarchy dataset.

Field name | type | Description
---|---|---
LocationNumber | int | Store identifier
Day | int | Calendar Day
SalesDayName | string | Day of the week
Coefficient0 | double | Constant
Coefficient1 | double | Slope
Coefficient2 | double | Curvature

##  3 Sales Demand Model

Given the stochastic nature of minute-to-minute demand, most of the value in
prediction comes from how this rate depends on “upper level” factors. As
mentioned at the upper levels the data is grouped by location and day of the
week. Our stochastic renewal process assumes that variation from minute to
minute is not predictable from preceding events. We approximate the
instantaneous demand rate function for each day, and model it by three
parameters:

  * $c0$

: the average daily demand,

  * $c1$

: the daily trend,

  * $c2$

: the “peakedness” of the rate.

We arrive at these values by regressing the centered daily sales arrivals
using a $2^{nd}$ order polynomial—equivalent to an orthonormal basis of the
first 3 Lagrange polynomials—one has 3 independent values that describe each
day and store, of which there are, as mentioned, 8649 instances. $c0$
coefficients tend to small positive values, as Figures3 and 4 show. $c1$ and
$c2$ coefficients are several orders of magnitude smaller, clustering around
zero. This suggests that the daily transaction rate is relatively constant
with minor tendencies to either rise or fall during the day. Each regression
might encompass a few hundred points, but some instances are sparse and
include just a few points. The ability for hierarchical Bayesian methods to
“share” information among units from upper levels improves predictions in
these cases.

The regression step is computationally inexpensive and, in a production
system, would be done at each location, offloading much of the central
computing load to make it possible to scale the model to 1000s of locations.
Certainly more sophisticated local models can be applied, but we use this
approximation that we pre-calculate of the lowest level of the hierarchy to
focus attention on the core problem.

###  3.1 The Graphical Model

The network diagram, Figure 5 shows the graphical model using plate notation.
Each upper level is represented by a separate overlapping plate. The shaded
inner plate is solved separately by the regression step and fed into the upper
level model. This approximation can be removed in the future, but it is
adequate to demonstrate the benefits of the hierarchy. In the statistical
literature this design is called a “two-way random effects model.”111See
Gelman and Hill [2007], p.245. Despite using it, the author argues against use
of the the term “random effects”, since the Bayesian formulation has no need
for the random effects – fixed effects distinction.

The model equation, with the distribution assumptions is described by these
four equations. The superscripts, $D$ and $J$ indicate the day and location
grouping levels. Any parameter without an explicit prior (e.g. $\mu$) is
assumed to have a uniform prior.

$\displaystyle y_{dj}$ $\displaystyle=z^{(D)}_{d}+z^{(J)}_{j}+\mu$ (1)
$\displaystyle z^{(D)}_{d}$
$\displaystyle\sim\mathcal{N}(\tau_{d},\sigma^{(D)})$ (2) $\displaystyle
z^{(J)}_{j}$ $\displaystyle\sim\mathcal{N}(\tau_{j},\sigma^{(J)})$ (3)
$\displaystyle y$ $\displaystyle\sim\mathcal{N}(y_{dj},\sigma)$ (4)

We run this model separately for each class of coefficients,
$(c0_{dj},c1_{dj},c2_{dj})$ by substituting them for the $y$ variable, to
create three versions of the model.

The inputs to the model are the $y$, indexed by location $j$ and day-of-week
$d$. The model estimates the distribution of the contribution of each day-of-
week $z^{(D)}$, and location $z^{(J)}$. Four global parameter distributions
are output for $\sigma^{(D)}$, $\sigma^{(J)}$, $\sigma$, and $\mu$. The MCMC
simulation of the model outputs marginal posteriors over each of the $4+7+49$
parameters conditioned on the observed data $y$.

The stan file for these equations in shown in Figure 9 in the Appendix.

Figure 5: The graphical model network for a model with two upper levels, one
for $J$ locations, the other for $D$ days of the week. Each node represents a
distribution in the generative model. Nodes on a plate are replicated by the
plate index. Overlapping plates indicate replication by the product of the
indexes. Subscripts for nodes on a plate are implied. Nodes outside any plate
are global parameters. Sales transaction data $x$ is modelled separately as 3
component parameters, implying 3 separate models for $y$.

###  3.2 Training

The stan code to run the model is shown in Figure 6, showing the sampling data
read as inputs, and the parameters (par) as outputs. An MCMC training run on
the train dataset of 4302 samples completes in less than 1 minute per chain on
a current laptop. We observed less-than-linear increases in processing time
with larger samples, when running on the combined train and test dataset.

Runtime for each model was less than 4 minutes. To assure convergence for $c1$
and $c2$ values the $y$ data was rescaled to have mean near 1.

⬇

model <- stan_model("regression_coef.stan")

rc_fit <- sampling(model,

data=c("N","D", "J", "d_id", "j_id", "y"), iter=4000, chains=4)

params <- summary(rc_fit,

pars=c("mu", "s_d", "s_j", "s_elipson", "z_d", "z_j"))$summary

Figure 6: rstan code to run the model.

Stan returns several diagnostics including tests for sampling convergence, for
which our model has no issues. Among other diagnostics is a log-likelihood
measure useful for comparisons among runs:

    # mean   se_mean       sd
    # lp__   1410.121     0.1888

Another sanity check is to see the estimates of variance, as one would in
conventional analysis of variance. We see that the sum of variances closely
approximates the $y$ variance.

$\displaystyle 0.445$
$\displaystyle=\sqrt{\sigma^{{(D)}^{2}}+\sigma^{{(J)}^{2}}+\sigma^{2}}$
compared to $\displaystyle 0.419$ $\displaystyle=\sigma_{y}$

The fraction of explained variation is

$R^{2}=1-\sigma^{2}/\sigma_{y}^{2}=0.638$

##  4 Evaluation

Despite the MCMC results returning the full sample marginals for all estimated
parameters, we use a conventional accuracy measure for prediction of the
coefficients. _Bias_ is simply the sum of the difference in the means of
predicted to actual. _RMSE_ is the square root of the average squared
deviation of the differences between predicted and actual. The model predicts
the three coefficient values for each of the upper-level units: the 49
locations and 7 days of the week. These can be compared to a baseline
prediction, of just using the average coefficient value for the hold-out set
over those units. Those averages are exactly the means of the boxplots shown
in Figures4 and 3 This is a challenging test; the average hold-out values are
themselves an accurate predictor.

Ideally one would want to evaluate improvements at the transactions level and
not only at the level of the coefficients fit to the daily transactions. As
argued, the variation of events around the instantaneous demand rate is so
large it would obscure any improvements to the coefficients that determine the
instantaneous demand rate.

First as a sanity check we compute _bias_ and _RMSE_ error for “test on
train”, for $c0$, by predicting from the values that were trained on. These
errors are negligible, as expected:

    Test on ’train’
       $bias
       -0.0003

       $rmse
       0.00876

Full evaluation tests were run for the hierarchical model, for both the
location predictions and day-of-week predictions for the 3 separate
components, and compared to baseline averages. The smaller of the errors
between the baseline average prediction and the hierarchical prediction is
shown in bold, in Table3.

Table 3: RMSE Errors.

Coefficient | Group | Average | Hierarchy
---|---|---|---
c0 | Location | 0.0398 | 0.0381
c0 | Day-Of-Week | 0.0277 | 0.252
c1 | Location | 2.860e-4 | 2.838e-4
c1 | Day-Of-Week | 6.61e-5 | 6.46e-5
c2 | Location | 1.012e-06 | 0.980e-6
c2 | Day-Of-Week | 2.858e-7 | 2.614e-7

Except for the Day-of-week prediction for $c0$ coefficients the hierarchical
model incrementally outperforms the baseline. Observing that using the average
value as a predictor obtains an RMSE error of just a few percent, this is an
impressive improvement by the hierarchical simulation model and recommends its
use in production.

Figure 7: Predicted coefficient0 versus location means from the actual test
set as a predictor. Each point is one location. The red line shows $x=y$.
Figure 8: Predicted coefficient1 versus location means from the actual test
set.

##  5 Discussion

###  5.1 Alternative approaches

Popular today are numerous deep-learning based time-series forecasting tools
such asKuznetsov and Mariet [2018], which considers the whether multiple
series should be learnt in a common model or separately. Although this
approach may approximate the ability to generalize over multiple (e.g. the 324
time-of-day location) series, it has none of the advantages of the method
presented. We believe that such a method invests vast effort in learning from
an stochastic event series, where simple smoothing methods suffice, at the
expense of analyzing the effects of groupings and the contextual variables on
the more parsimonious upper-level data where interesting effects occur. It
would scale to the thousands of locations at only comparatively large
computational cost.

An “off-the-shelf” option for hierarchical modelling is to use a current
spatial-temporal forecasting tool instead, such as FOST. FOST combines
traditional and ML-based temporal modeling with spatial modeling via graph
neural networks. FOST requires two data inputs: time series data and “spatial”
data in the form of a list of directed graph edges and their weights. Edges
can represent truly spatial data, such as distances or travel times.
Additionally, edges can represent relationships such as correlations between
the items represented by graph nodes. If the nodes represent products for
sale, then edges can represent cannibalization or halo effects. A halo effect
is when reducing the price of one product increases sales of a related
product—a cannibalization effect is the opposite.

As for spatial effects, it depends on how close together stores are. Thus we
could create a graph of locations, with edges weighted by distance or travel
time. Alternatively, we could create a graph with edges weighted by the
correlation of sales between pairs of stores.

Effects among locations expressed on a graph structure is a level of model
complexity not captured by a hierarchy. Relations in a hierarchy are mediated
by presence of common parents, not by pair-wise relations. The choice of model
depends speculatively on the relevance of relations that would be expressed in
a graph, the feasibility of predicting from event process data, and
scalability to the number locations.

###  5.2 Discussion

The apparent non-stationarity of the current dataset has been managed by
ignoring the temporal separation of test and train splits. Arguably this is of
minor consequence because the groupings by day-of-week and location averaged
over the year do not relate to each other temporally.

We have many ideas how to extend the current model that at this point just
demonstrates the feasibility of this approach. Extensions of the Graphical
Model could bring in temporal predictors to see if they do have any effect.
One could test if the conditional independence of day-of-week and location is
valid, but a priori there is little reason to think it is not. Most
interestingly one could introduce contextual variables to the upper levels for
economic, demographic, weather, holiday, spatial and other effects that are
presumed determinants of demand.

Simulation methods such as MCMC do not always converge in fixed time, so a
more reliable inference approximation method would be preferable. Using
Expectation Propagation in place of MCMC, as for instance implemented by
_infer.net_ Wang and Wand [2011] is worth trying.

By only forecasting recorded transactions we overlook demand by customers
turned away because no items were available for sale. These events are not
captured. This biases estimates downward. An outstanding challenge in
forecasting is to predict unmet demand, so one would know that a sale would be
missed because items were not available. One possible way to address this is
to model the censoring process. We assume a distribution for transactions, say
Poisson, or a more general renewal distribution, then “extrapolate” the
distribution to estimate extreme values not seen in the data.

Alternatively if one had sales transactions and production amounts data then
the the full demand distribution could be estimated. Then one could tell when
production limits sales (e.g. when a store is sold-out). So by analyzing sales
data when stores are not sold out the full demand curve could be fit to
estimate unmet sales and remove the source of forecast bias.

An obvious oversight in this model formation is that lack of a “supply” model
to complement the “demand” model. A comprehensive model would recommend the
optimal production schedule for perishable items. For that we’d need a
corresponding production model to complement this demand model. The value of
improvements due to the model cannot be evaluated without also observing and
modeling real-time production in which case one would also need historical
data on actual production.

###  5.3 Conclusion

We’ve demonstrated promising results on an initial framework for a
hierarchical model of sales transactions over groupings of locations and days.
As in many applied modelling problems, there are natural groupings of the
samples that can be expressed hierarchically, avoiding the need to pre-process
the data to de-normalize it such as done by “one-hot” encoding. Such problems
translate directly into Bayesian Graphical Models for which there are now
mature solution tools. The advantages of formulating hierarchical as Graphical
Models are several. 1. Groups with limited data benefit from inference from
similar neighboring groups, increasing the effective sample size of groups.
The degree of interaction is a consequence of the sample size within groups.
2. Co-variates can be applied directly to the relevant upper levels. 3. Models
may be federated at lower levels for distributed processing. 4. The full value
of Bayesian methods in terms of explainability and the derivation of true
predictive distributions can be exploited.

###### Acknowledgements.

Without the express help of my data science colleagues, and numerous sales,
marketing, and engineering staff both in-house and at the client, none of this
would be possible.

## References

  * Geffner et al. [2022] Tomas Geffner, Javier Antoran, Adam Foster, Wenbo Gong, Chao Ma, Emre Kiciman, Amit Sharma, Angus Lamb, Martin Kukla, Nick Pawlowski, Miltiadis Allamanis, and Cheng Zhang.  Deep end-to-end causal inference, 2022.  URL https://arxiv.org/abs/2202.02195.
  * Gelman and Hill [2007] Andrew Gelman and Jennifer Hill.  _Data Analysis Using Regression and Multilevel / Hierarchical Models_.  Cambridge University Press, 2007.
  * Gelman et al. [2013] Andrew Gelman, John B. Carlin, Hal S. stern, David B. Dusson, Aki Vehtari, and Donald B. Rubin.  _Bayesian Data Analysis, Third Edition_.  CRC Press, 2013.
  * Getoor and Tasker [2007] Lise Getoor and Ben Tasker.  _Introduction to Statistical Learning Theory_.  MIT Press, 2007.
  * Koller and Friedman [2009] Daphne Koller and Nir Friedman.  _Probabilistic Graphical Models_.  MIT Press, 2009.
  * Kuznetsov and Mariet [2018] Vitaly Kuznetsov and Zelda Mariet.  Foundations of sequence-to-sequence modeling for time series, 2018.  URL https://arxiv.org/abs/1805.03714.
  * Pearl [2009] Judea Pearl.  _Causality: Models, Reasoning and Inference_.  Cambridge University Press, USA, 2nd edition, 2009.  ISBN 052189560X.
  * Porteus [2002] E. L. Porteus.  _Foundations of Stochastic Inventory Theory_.  Stanford University Press, 2002.
  * Wang and Wand [2011] S. S. J. Wang and M. P. Wand.  Using infer.net for statistical analyses.  _The American Statistician_ , 65(2):115–126, 2011\.  10.1198/tast.2011.10169.  URL https://doi.org/10.1198/tast.2011.10169.

## Appendix

⬇

// This Stan program defines a simple model, with a two separate

// groups

// – J locations

// – D days of the week.

data {

int<lower=0> N; // row count

int<lower=0> D; // The number of days of the week

int<lower=0> J; // The number of locations

int<lower=1, upper=D> d_id[N]; // map rows into days of week

int<lower=1, upper=J> j_id[N]; // map rows into locations.

vector[N] y; // outcome – the predicted value for that dow & location.

}

// The parameters accepted by the model. Our model

// includes a eta and sigma for each effect

parameters {

real mu;

vector[D] eta_d;

vector[J] eta_j;

real<lower=0> s_d;

real<lower=0> s_j;

real<lower=0> s_elipson;

}

transformed parameters{

vector[J] z_j;

vector[D] z_d;

vector[N] yhat;

z_j = s_j * eta_j; // random effects

z_d = s_d * eta_d;

for (i in 1:N) {

yhat[i] <- mu + z_d[d_id[i]] + z_j[j_id[i]];

}

}

// The model to be estimated.

model {

eta_d ~ normal(0,1);

eta_j ~ normal(0,1);

y ~ normal(yhat, s_elipson);

}

Figure 9: File regression_coef.stan of a two-way random effects model