Some losses, such as the loss of reputation, are difficult to estimate owing
to the lack of historical data. However, the proposed framework would still be
applicable if those losses were quantified using alternative approaches. For
example, reputational losses could possibly be estimated using a business
analytics approach based on assumptions about the impact of reputation on
future cash flows or equity value changes. For other losses, of which the
stakes are usually considered too high to analyze from a cost-benefit
perspective, such as impacts on national security, a cost-benefit argument can
be made nevertheless by stating that benefits are so high that any cost is
worthwhile.

Evidence for the practicality and usefulness of the proposed framework can be
found in some already existed cyber risk management products and services in
the industry. Other than the aforementioned tool built by [3], KPMG has also
developed a scenario-based cyber risk modeling and quantification model, which
has a structure similar to the cascade model used in this study (see [32]).
Comparing to the existing solutions in the market, the framework proposed in
this paper could provide a scientific foundation and hence offer more reliable
and convincing results in the longer run.

Several parts of this paper could certainly be further studied. First, we
assume that financial losses due to a vulnerability can be linearly reduced by
implementing a control. It is possible in practice that the control does not
scale down the losses but increases the probability of zero loss (see, such
as, [31]); in this case, the loss tensor is no longer an element-wise
multiplication of the cascade model tensor and the raw loss tensor. But
rather, one needs to modify the impact with a point mass at zero loss. Second,
the case study in this paper is based on the publicly acquirable dataset
developed and maintained by Advisen Ltd., which could not reveal sensitive
information of victim companies like internal configuration for
vulnerabilities, implemented controls, and assets, as well as external
exposure for threats. Therefore, we make assumptions for illustrative purposes
of the case study. In practice, when implementing the proposed framework,
companies could avoid making these assumptions and establish the links among
threats, vulnerabilities, and assets based on their own actual cybersecurity
conditions. Third, we mentioned that the adopters of this framework may need
to frequently revise the risk assessment and capital allocation plan due to
substantial changes in the cyber risk environment. Alternatively, a multi-
period model can potentially be built to account for those changes and how
cyber risk management plans should be updated accordingly. Lastly, given that
the purpose of this study is managing the cyber risk of a single organization,
from the organization’s own perspective, we treat inter-organizational
incidents, such as the spread of viruses, the same way as other incidents
without inter-organizational dependence. Future studies may take a multi-agent
approach to the cyber risk management for a group of organizations with the
loss correlation among them, such as via common vulnerabilities, taken into
consideration.

## References

  * [1] Abdulmajeed Alahmari and Bob Duncan  “Cybersecurity Risk Management in Small and Medium-Sized Enterprises: A Systematic Review of Recent Evidence”  In _2020 International Conference on Cyber Situational Awareness, Data Analytics and Assessment (CyberSA)_ IEEE, 2020
  * [2] David L. Alderson, Gerald G. Brown and William M. Carlyle  “Assessing and Improving Operational Resilience of Critical Infrastructures and Other Systems”  In _Bridging Data and Decisions_ INFORMS, 2014, pp. 180–215
  * [3] Arx Nimbus, LLC.  “The 2021 Digital Risk Checklist”, 2021  URL: https://www.arxnimbus.com/cyber-security-risk-whitepaper-download
  * [4] Sarah Bensalem, Nicolás Hernández-Santibáñez and Nabil Kazi-Tani  “A Continuous-Time Model of Self-Protection”  In _Finance and Stochastics_ 27.2, 2023, pp. 503–537
  * [5] Sarah Bensalem, Nicolás Hernández Santibáñez and Nabil Kazi-Tani  “Prevention Efforts, Insurance Demand and Price Incentives Under Coherent Risk Measures”  In _Insurance: Mathematics and Economics_ 93, 2020, pp. 369–386
  * [6] Rainer Böhme and Gaurav Kataria  “Models and Measures for Correlation in Cyber-Insurance”  In _5th Annual Workshop on the Economics of Information Security, WEIS 2006_ , 2006  URL: http://weis2006.econinfosec.org/docs/16.pdf
  * [7] Rainer Böhme, Stefan Laube and Markus Riek  “A Fundamental Approach to Cyber Risk Analysis”  In _Variance_ 12.2, 2019, pp. 161–185
  * [8] Huseyin Cavusoglu, Birendra Mishra and Srinivasan Raghunathan  “The Effect of Internet Security Breach Announcements on Market Value: Capital Market Reactions for Breached Firms and Internet Security Developers”  In _International Journal of Electronic Commerce_ 9.1 Informa UK Limited, 2004, pp. 70–104
  * [9] James Cebula, Mary Popeck and Lisa Young  “A Taxonomy of Operational Cyber Security Risks Version 2”, 2014
  * [10] Xiaowei Chen, Wing Fung Chong, Runhuan Feng and Linfeng Zhang  “Pandemic Risk Management: Resources Contingency Planning and Allocation”  In _Insurance: Mathematics and Economics_ 101B Elsevier BV, 2021, pp. 359–383
  * [11] Wing Fung Chong, Runhuan Feng and Longhao Jin  “Holistic Principle for Risk Aggregation and Capital Allocation”  In _Annals of Operations Research_ , 2021
  * [12] Wing Fung Chong, Daniel Linders, Zhiyu Quan and Linfeng Zhang  “Incident-Specific Cyber Insurance”, 2023  arXiv:2308.00921 [q-fin.RM]
  * [13] CIS  “CIS Controls”, 2021  Center for Internet Security, Inc.  URL: https://learn.cisecurity.org/cis-controls-download
  * [14] Matteo Crosignani, Marco Macchiavelli and André F. Silva  “Pirates Without Borders: The Propagation of Cyberattacks Through Firms’ Supply Chains”  In _Journal of Financial Economics_ 147.2, 2023, pp. 432–448
  * [15] Jan Dhaene, Andreas Tsanakas, Emiliano A. Valdez and Steven Vanduffel  “Optimal Capital Allocation Principles”  In _Journal of Risk and Insurance_ 79.1, 2012, pp. 1–28
  * [16] Daniel DiMase, Zachary A. Collier, Kenneth Heffner and Igor Linkov  “Systems Engineering Framework for Cyber Physical Security and Resilience”  In _Environment Systems and Decisions_ 35.2 Springer ScienceBusiness Media LLC, 2015, pp. 291–300
  * [17] Christian Eckert, Nadine Gatzert and Madeline Schubert  “Analyzing spillover effects from data breaches to the US (cyber) insurance industry”  In _The European Journal of Finance_ 0.0, 2022, pp. 1–24
  * [18] Isaac Ehrlich and Gary S. Becker  “Market Insurance, Self-Insurance, and Self-Protection”  In _Journal of Political Economy_ 80.4 University of Chicago Press, 1972, pp. 623–648
  * [19] Thomas M. Eisenbach, Anna Kovner and Michael Junho Lee  “Cyber risk and the U.S. financial system: A pre-mortem analysis”  In _Journal of Financial Economics_ 145.3, 2022, pp. 802–826
  * [20] Martin Eling and Kwangmin Jung  “Copula Approaches for Modeling Cross-sectional Dependence of Data Breach Losses”  In _Insurance: Mathematics and Economics_ 82 Elsevier BV, 2018, pp. 167–180
  * [21] Martin Eling and Jan Wirfs  “What Are the Actual Costs of Cyber Risk Events?”  In _European Journal of Operational Research_ 272.3 Elsevier BV, 2019, pp. 1109–1119
  * [22] Mario Frustaci, Pasquale Pace, Gianluca Aloi and Giancarlo Fortino  “Evaluating Critical Security Issues of the IoT World: Present and Future Challenges”  In _IEEE Internet of Things Journal_ 5.4 Institute of ElectricalElectronics Engineers (IEEE), 2018, pp. 2483–2495
  * [23] Edward Furman and Ričardas Zitikis  “Weighted Risk Capital Allocations”  In _Insurance: Mathematics and Economics_ 43.2, 2008, pp. 263–269
  * [24] Lawrence A. Gordon and Martin P. Loeb  “The Economics of Information Security Investment”  In _ACM Transactions on Information and System Security_ 5.4, 2002, pp. 438–457
  * [25] Lawrence A. Gordon, Martin P. Loeb, William Lucyshyn and Lei Zhou  “Increasing Cybersecurity Investments in Private Sector Firms”  In _Journal of Cybersecurity_ 1.1, 2015, pp. 3–17
  * [26] Lawrence A. Gordon, Martin P. Loeb and Lei Zhou  “Investing in Cybersecurity: Insights from the Gordon-Loeb Model”  In _Journal of Information Security_ 07.2 Scientific Research Publishing, Inc,, 2016, pp. 49–59
  * [27] Hartford Steam Boiler  “HSB CyberOneTM Coverage”, 2017  URL: https://www.munichre.com/content/dam/munichre/contentlounge/website-pieces/documents/CyberOne_HSB.pdf/_jcr_content/renditions/original./CyberOne_HSB.pdf
  * [28] Hemantha SB Herath and Tejaswini C Herath  “Copula-Based Actuarial Model for Pricing Cyber-Insurance Policies”  In _Insurance Markets and Companies_ 2.1 Business Perspectives Ltd., 2011
  * [29] Annette Hofmann and Richard Peter  “Self-Insurance, Self-Protection, and Saving: On Consumption Smoothing and Risk Management”  In _Journal of Risk and Insurance_ 83.3 Wiley, 2015, pp. 719–734
  * [30] Alain Janvry and Elisabeth Sadoulet  “Optimal Share Contracts with Moral Hazard on Effort and in Output Reporting: Managing the Double Laffer Curve Effect”  In _Oxford Economic Papers_ 59.2 Oxford University Press (OUP), 2006, pp. 253–274
  * [31] Shinichi Kamiya, Jun-Koo Kang, Jungmin Kim, Andreas Milidonis and René M. Stulz  “Risk management, firm reputation, and the impact of successful cyberattacks on target firms”  In _Journal of Financial Economics_ 139.3, 2021, pp. 719–749
  * [32] KPMG  “Cyber Risk Modelling and Quantification”, 2020  URL: https://assets.kpmg.com/content/dam/kpmg/uk/pdf/2020/11/cyber-risk-modelling-and-quantification.pdf
  * [33] Zhaoxi Liu, Wei Wei, Lingfeng Wang, Chee-Wooi Ten and Yeonwoo Rho  “An Actuarial Framework for Power System Reliability Considering Cybersecurity Threats”  In _IEEE Transactions on Power Systems_ , 2020
  * [34] Alessandro Mazzoccoli and Maurizio Naldi  “Robustness of Optimal Investment Decisions in Mixed Insurance/Investment Cyber Risk Management”  In _Risk Analysis_ 40.3, 2020, pp. 550–564
  * [35] Mario Menegatti  “Optimal prevention and prudence in a two-period model”  In _Mathematical Social Sciences_ 58.3 Elsevier BV, 2009, pp. 393–397
  * [36] NIST  “Framework for Improving Critical Infrastructure Cybersecurity”, 2018  National Institute of StandardsTechnology  URL: https://nvlpubs.nist.gov/nistpubs/CSWP/NIST.CSWP.04162018.pdf
  * [37] Chen Peng, Maochao Xu, Shouhuai Xu and Taizhong Hu  “Modeling Multivariate Cybersecurity Risks”  In _Journal of Applied Statistics_ 45.15, 2018, pp. 2718–2740
  * [38] Keyun Ruan  “The Point of Diminishing Return on Cyber Risk Investment”  In _Digital Asset Valuation and Cyber Risk Management_ Elsevier, 2019, pp. 99–115
  * [39] SANS Institute  “Mapping the Critical Security Controls (CSC) v4.0 to the Verizon VERIS 2013 Threats - Executive Summary”, 2013  URL: https://www.sans.org/media/critical-security-controls/critical_security_controls_v4.0.pdf
  * [40] Fabrizio Smeraldi and Pasquale Malacaria  “How to Spend It: Optimal Investment for Cyber Security”  In _Proceedings of the 1st International Workshop on Agents and CyberSecurity_ ACM Press, 2014
  * [41] Gernot Tragler, Jonathan P. Caulkins and Gustav Feichtinger  “Optimal Dynamic Allocation of Treatment and Enforcement in Illicit Drug Control”  In _Operations Research_ 49.3 Institute for Operations Researchthe Management Sciences (INFORMS), 2001, pp. 352–362
  * [42] University of Oxford and Axis  “White Paper: Calculating Residual Cyber-risk”, 2020  University of OxfordAxis  URL: https://www.cs.ox.ac.uk/projects/ACVAR/White%20Paper%20-%20Calculating%20residual%20cyber-risk.pdf
  * [43] Verizon Security Research & Cyber Intelligence Center  “Vocabulary for Event Recording and Incident Sharing (VERIS)”, 2021  URL: https://github.com/vz-risk/veris
  * [44] Aliza Vigderman and Gabe Turner  “The Best Business Antivirus Protection Software in 2023”, 2023  URL: https://www.security.org/antivirus/best/business/
  * [45] Wenye Wang and Zhuo Lu  “Cyber Security in the Smart Grid: Survey and Challenges”  In _Computer Networks_ 57.5, 2013, pp. 1344–1371
  * [46] Jude Wanniski  “Taxes, Revenues, and the Laffer Curve”  In _The Public Interest_ 50, 1978, pp. 3–16
  * [47] Maochao Xu, Kristin M. Schweitzer, Raymond M. Bateman and Shouhuai Xu  “Modeling and Predicting Cyber Hacking Breaches”  In _IEEE Transactions on Information Forensics and Security_ 13.11 Institute of ElectricalElectronics Engineers (IEEE), 2018, pp. 2856–2871
  * [48] Gabriela Zeller and Matthias Scherer  “Risk Mitigation Services in Cyber Insurance: Optimal Contract Design and Price Structure”  In _The Geneva Papers on Risk and Insurance - Issues and Practice_ 48.2, 2023, pp. 502–547

##  Appendix A Mapping Relationship between Control and Asset

Control Name | Assets
---|---
Devices | Applications | Users | Network | Data
1\. Inventory and Control of Hardware Assets | 1 |  |  |  |
2\. Inventory and Control of Software Assets |  | 1 |  |  |
3\. Continuous Vulnerability Management |  | 1 | 1 |  |
4\. Controlled Use of Administrative Privileges |  |  | 1 |  |
5\. Secure Configuration for Hardware and Software on Mobile Devices, Laptops, Workstations and Servers |  | 1 |  |  |
6\. Maintenance, Monitoring and Analysis of Audit Logs |  |  |  | 1 |
7\. Email and Web Browser Protections |  | 1 |  | 1 |
8\. Malware Defenses | 1 |  |  | 1 |
9\. Limitation and Control of Network Ports, Protocols and Services | 1 |  |  |  |
10\. Data Recovery Capabilities |  |  |  |  | 1
11\. Secure Configuration for Network Devices, such as Firewalls, Routers and Switches |  |  |  | 1 |
12\. Boundary Defense | 1 |  | 1 | 1 |
13\. Data Protection |  |  |  |  | 1
14\. Controlled Access Based on the Need to Know |  |  |  | 1 | 1
15\. Wireless Access Control | 1 |  |  | 1 |
16\. Account Monitoring and Control |  |  | 1 |  |
17\. Implement a Security Awareness and Training Program |  |  | 1 |  |
18\. Application Software Security |  | 1 |  |  |
19\. Incident Response and Management | 1 | 1 |  | 1 | 1
20\. Penetration Tests and Red Team Exercises |  |  |  |  |

Table 9: Controls and assets defined in the CIS Controls. Value 1 represents
that the control mitigates the vulnerability in the corresponding asset. A
blank cell means there is no relationship between the control and the asset.

##  Appendix B Mapping Relationship between Threat and Control

Threats | CIS Controls |  |  |
---|---|---|---|---
1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20
Tampering (alter physical form or function) | 1 |  | 1 |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  |
Backdoor (enable remote access) |  | 1 | 1 | 1 | 1 |  |  |  |  |  |  | 1 | 1 |  |  |  |  |  |  |
Use of stolen authentication credentials |  |  |  |  |  |  |  |  |  |  |  |  |  | 1 | 1 |  |  |  |  |
Export data to another site or system |  |  |  |  |  |  |  |  |  | 1 |  |  | 1 |  |  |  | 1 |  |  |
Use of Backdoor or C2 channel |  | 1 | 1 | 1 | 1 |  |  |  |  |  |  | 1 | 1 |  |  |  |  |  |  |
Phishing (or any type of *ishing) |  | 1 | 1 | 1 | 1 |  |  |  |  |  |  | 1 | 1 |  |  |  |  |  |  |
Command and control (C2) |  | 1 | 1 | 1 | 1 |  |  |  |  |  |  | 1 | 1 |  |  |  |  |  |  |
Downloader (pull updates or other malware) |  | 1 | 1 | 1 | 1 |  |  |  |  |  |  | 1 |  |  |  |  |  |  |  |
Brute force or password guessing attacks |  |  |  |  |  |  |  |  |  |  |  |  |  | 1 | 1 |  |  |  |  |
Spyware, keylogger or form-grabber (capture user input or activity) |  | 1 | 1 | 1 | 1 |  |  |  |  |  |  | 1 |  |  |  |  |  |  |  |
Capture data stored on system disk |  |  |  |  |  |  |  |  |  |  |  |  |  | 1 | 1 |  | 1 |  |  |
System or network utilities (e.g., PsTools, Netcat) |  | 1 | 1 | 1 | 1 |  |  |  |  |  |  | 1 |  |  |  |  |  |  |  |
Abuse of system access privileges |  |  |  |  |  |  |  |  |  |  |  | 1 |  | 1 | 1 | 1 |  |  |  |
Ram scraper or memory parser (capture data from volatile memory) |  | 1 | 1 | 1 | 1 |  |  |  |  |  |  | 1 |  |  |  |  |  |  |  |
Use of unapproved hardware or devices | 1 |  |  |  |  |  |  |  |  |  |  |  |  |  |  |  | 1 |  |  |
SQL injection |  |  |  |  |  | 1 |  |  |  |  |  |  |  |  |  |  |  |  |  |
Embezzlement, skimming, and related fraud |  |  |  |  |  |  |  |  |  |  |  |  |  | 1 | 1 |  |  |  |  |
Theft (taking assets without permission) | 1 |  |  |  |  |  |  |  |  |  |  |  |  | 1 |  |  | 1 |  |  |
Bribery or solicitation |  |  |  |  |  |  |  |  |  |  |  | 1 |  | 1 | 1 |  |  |  |  |
Disable or interfere with security controls |  |  | 1 |  |  |  |  |  |  | 1 |  |  |  |  |  |  |  |  |  |
Password dumper (extract credential hashes) |  | 1 | 1 | 1 | 1 |  |  |  |  |  |  | 1 |  |  |  |  |  |  |  |
Misconfiguration |  |  | 1 |  |  |  |  |  |  | 1 |  |  |  |  |  |  |  |  |  |
Programming error (flaws or bugs in custom code) |  |  |  |  |  | 1 |  |  |  |  |  |  |  |  |  |  |  |  |  |
Misdelivery (direct or deliver to wrong recipient) |  |  |  |  |  | 1 |  |  |  |  |  |  |  |  | 1 |  | 1 |  |  |
Loss or misplacement | 1 |  | 1 |  |  |  |  |  |  |  |  |  |  |  |  |  | 1 |  |  |
Table 10: The mapping between threats and CIS Controls from [39].

##  Appendix C Algorithm for Finding Set $\mathfrak{I}$

Input: $\omega_{ik}$ and
$\overline{K}_{ik}\left(\mathbf{M},{\mathbf{I}}\right)$, for $i=1,2,\dots,l$,
$k=1,2,\dots,n$; $b\left(\mathbf{M},{\mathbf{I}}\right)$

Output: $\mathfrak{I}$

1

2$\ell\leftarrow l\times n-1$

3${\mathfrak{I}}^{0}\leftarrow$ list of indices of sorted
$\omega_{ik}\overline{K}_{ik}\left(\mathbf{M},{\mathbf{I}}\right)$, such that
$\omega_{{\mathfrak{I}}^{0}[0]}\overline{K}_{{\mathfrak{I}}^{0}[0]}\left(\mathbf{M},{\mathbf{I}}\right)\leq\omega_{{\mathfrak{I}}^{0}[1]}\overline{K}_{{\mathfrak{I}}^{0}[1]}\left(\mathbf{M},{\mathbf{I}}\right)\leq\dots\leq\omega_{{\mathfrak{I}}^{0}[\ell]}\overline{K}_{{\mathfrak{I}}^{0}[\ell]}\left(\mathbf{M},{\mathbf{I}}\right)$

4for  _$\iota\leftarrow 0$ to $\ell$_ do

5 ${\mathfrak{I}}\leftarrow{\mathfrak{I}}^{0}[\iota,\dots,\ell]$

6 if
_$b\left(\mathbf{M}{,\mathbf{I}}\right)+\sum_{\tau=0}^{\ell-\iota}\left(\frac{\omega_{{\mathfrak{I}}[0]}}{\omega_{{\mathfrak{I}}[\tau]}}\overline{K}_{{\mathfrak{I}}[0]}\left(\mathbf{M},{\mathbf{I}}\right)-\overline{K}_{{\mathfrak{I}}[\tau]}\left(\mathbf{M},{\mathbf{I}}\right)\right)\geq
0$_ then

7 break

8 end if

9

10 end for

return _${\mathfrak{I}}$_

Algorithm 1 Identification Procedure of Set $\mathfrak{I}$


XX XX 8 Month Publication Title

Department: Reproducible Research

# Managing Software Provenance to Enhance Reproducibility in Computational
Research

Akash Dhruv  Anshu Dubey  Argonne National Laboratory, Lemont, IL, 60439, USA

(2023)

###### Abstract

Scientific processes rely on software as an important tool for data
acquisition, analysis, and discovery. Over the years sustainable software
development practices have made progress in being considered as an integral
component of research. However, management of computation-based scientific
studies is often left to individual researchers who design their computational
experiments based on personal preferences and the nature of the study. We
believe that the quality, efficiency, and reproducibility of computation-based
scientific research can be improved by explicitly creating an execution
environment that allows researchers to provide a clear record of traceability.
This is particularly relevant to complex computational studies in high-
performance computing (HPC) environments. In this article, we review the
documentation required to maintain a comprehensive record of HPC computational
experiments for reproducibility. We also provide an overview of tools and
practices that we have developed to perform such studies around Flash-X, a
multi-physics scientific software.

††journal: Computing in Science & Engineering

## INTRODUCTION

Experimental and observational scientists follow a rigorous process of
recording their work. For many, receiving training to do this diligently is
considered a high priority for several excellent reasons. Unless the exact
conditions of the experiment are recorded it cannot be repeated. A complete
specification of assumptions made during the experiment, and the recording of
the methodology followed for interpreting and analyzing the outcome is
required for reproducibility – a hallmark of the scientific process.
Additional factors may include the expense of running experiments, and the
impact placement of probes might have on the recorded data. In observational
sciences, similar or greater constraints may result from the rarity of events
– the observer must be ready to record the event when it occurs. To make the
most of the event the observer would typically wish to record as many details
as possible.

This scientific discipline is largely lacking in computational science. It has
traditionally been assumed in computational science circles that careful
recording of experiments may not be necessary because running the software
again with identical input will produce identical results. This is a false
assumption because software is rarely stationary. Every instance of using it
to obtain new scientific results tweaks it in some ways, thereby changing some
of its behavior. In high-performance computing (HPC) environments, the focus
of this article, changes in the system software stack and hardware may also
cause changes in results generated by the software. As the scale and
complexity of the software and the problems it is trying to solve grow, so
does the expense of running computational experiments. Some high-profile
scientific failures in various disciplines and scrutiny of scientific output
during COVID-19 have put a spotlight on the rigor of running and recording
computational experiments. Although workflow management tools have been
developed to make it easier to run simulations, and archive and analyze the
data, there is no community-wide adoption of rigor in requiring provenance of
the code and data, and the log of the experiments conducted.

Laboratory notebooks, or scientific notebooks (referred to as lab notebooks
from here on) have long been used as a crucial tool in experimental and
observational sciences by researchers, engineers, and students to record and
document their work, experiments, observations, and findings in a structured
and organized manner. Their purpose is to maintain a detailed and accurate
record of all activities related to a research project or experiment.
Maintaining a well-kept lab notebook is not only considered a scientific best
practice but also a requirement in many research and academic settings. It
helps ensure transparency, accountability, and the reliability of scientific
work, which is essential for advancing knowledge and innovation.

Individual researchers often use a directory-based organization for their
computational experiments which are designed based on their personal
preferences, the nature of the studies, and the need to balance time and
computing resources. This directory-based design is helpful in nesting
experiments that cover different parametric spaces. Leveraging this design and
explicitly creating an execution environment around software and its
dependencies can improve the quality, efficiency, and reproducibility of their
studies. Such an environment, which is analogous to a laboratory space, can
help focus researchers to prioritize scientific rigor and develop tools and
practices to manage the generation, protection, storage, and analysis of data.
The construction of this virtual laboratory environment should include the
integration of a lab notebook that can provide a record of software
configurations and research decisions.

The question then is, what do lab notebooks look like for an HPC computational
experiment, and how should they be managed? In this paper, we list features
that lab notebooks may require for such computational investigations, and
describe one exemplar of a solution developed for conducting experiments with
Flash-X [1], a multi-physics multi-component software that can be used for
simulations in several science domains.

## LABORATORY NOTEBOOKS FOR COMPUTATIONAL EXPERIMENTS

Since computational work often involves coding, simulations, and data analysis
using software tools, the lab notebook is best maintained in electronic
format, and could include a combination of repositories, spreadsheets, text,
and markdown formats. For reasons discussed
previously111https://www.youtube.com/watch?v=fWpI4S_dvhc, the use of some form
of lab notebook in computational science has several challenges. However,
fundamental sets of activities that must be recorded are fairly simple to
enumerate.

  * •

Title and Purpose: Title or project name and a brief description of its
purpose or objectives.

  * •

Code Repository Links: Links or references to the code repositories and/or
input files used in the experiment.

  * •

Software and Hardware: Specifications of hardware used and system software
stack, libraries, and tools used in the experiment, along with their versions.

  * •

Modifications: Systematic recording of modifications to the software and
hardware during the experiment

  * •

Experiment Design: Description of the algorithms, data sets, parameters, and
any assumptions made for the experiment along with the reasoning behind the
selection of specific runs to be made. A log of steps taken to prepare for the
experiment can be very helpful for future experiments.

  * •

Data Sources: Documentation about any external data used, along with
explanations of how and where it’s stored. Any data preprocessing steps, such
as cleaning, normalization, or transformation.

  * •

Data Storage: Archival storage of produced and collected data with attached
metadata to be usable by other researchers.

  * •

Experimental Runs: Log of each experimental run, including the input
parameters, the date and time of execution, and the resulting output or data.
If version control is in use commit references for any significant updates
should be kept.

  * •

Results and Analysis: Presentation of the results of a computational
experiment, including tables, graphs, and statistical analyses along with a
description of how the results were interpreted and what conclusions were
drawn.

  * •

References: Citations and references to relevant literature, software
documentation, or external resources that influenced the computational
experiment.

Several tools and software applications can help maintain well-organized and
effective lab notebooks for computational experiments. These tools are
designed to streamline documentation, code management, data analysis, and
collaboration. For example, Jupyter Notebooks222https://jupyter.org/ are
interactive, web-based environments for creating and sharing documents that
combine live code (Python, R, etc.), equations, visualizations, and narrative
text. They are already widely used in data science and computational research.
Similarly R Markdown333https://rmarkdown.rstudio.com/ is an authoring format
that integrates code, results, and narrative text into a single document. It’s
commonly used with the R programming language but can be adapted for other
languages as well. Python-based notebooks such as Google
Colab444https://colab.google/ are similar to Jupyter Notebooks but are focused
on Python. They allow documentation of Python code alongside explanations and
visualizations. Several other open-source and commercial solutions exist with
a variety of features. Any of these tools along with a set of well-defined
recording practices can form the basis of an execution environment that
promotes reproducibility of computational experiments. However, none of them
suffice for complex multiphysics HPC computations which generate a huge amount
of data, and require substantial post-processing and analysis.

## THE FLASH-X SOLUTION

Flash-X is community-developed software that is undergoing several modes of
development simultaneously. It is a new version of a long-existing community
code FLASH [2] that has been re-architected to be compatible with
heterogeneous hardware platforms. Several new physics capabilities and an
entirely new method of integration have also been added to the code. As a
consequence, a fairly common occurrence is where ongoing capability and
performance improvement requirements collide with the needs of a domain
science study. The situation is exacerbated when all participants involved in
a study are not familiar with the inner workings and constraints of the code.
In general two different types of experiments are regularly conducted with the
code. One set of experiments measures the performance of different components
of the code, while the other set pertains to domain science investigations.
The requirements and constraints of these experiments differ from one another,
though some of the required meta-information for complete specification is
identical for all experiments. Some of this required meta-information such as
repository version, software stack version, and all the configuration
parameters are recorded in a log file that every execution instance of Flash-X
generates.

The performance experiments typically involve scaling studies with the
implication that the same application is run with different configurations of
degree of parallelism, hardware components in use, and possibly different
implementations of some of the code components. The exploration parameters in
such studies tend to be related to the infrastructural components of the code,
and the analysis is performed on either the performance summary section of the
log file or, if a performance tool is used, then the data generated by the
tool. The output of the simulation itself is not relevant to the study. The
domain science investigations explore the parameter space of the physics
involved. These experiments typically have little variation in the degree of
parallelism or hardware in use. The output of the simulation is an important
artifact here and can be quite large.

An example of a workflow that might be encountered by a developer who is also
a user of the code for doing scientific investigations is shown in Figure 1a.
Such developers may perform regular testing of the physics and infrastructure
components, keep track of the performance changes, and conduct scientific
experiments at the same time. The collection of lab notebooks [3] described in
Figure 1b provides a good example of what may be done in such circumstances.
These notebooks are all seeded from a general notebook [4] being maintained by
the developer in question. In this instance, the developer is simultaneously
conducting the following experiments,

  * •

Flow-Boiling-3D: Lab notebook for production runs of three-dimensional
multiphase flow boiling simulations.

  * •

ImBound-Mapping-Optimization: Lab notebook to investigate and optimize mapping
of Lagrangian particles on block-structured AMR grids for the purpose of
developing an immersed boundary method for fluid-structure interaction
problems.

  * •

AMReX-Bittree-Performance: Lab notebook to improve scaling of Flash-X
applications that use AMReX in octree mode.

  * •

Outflow-Forcing-BubbleML: Reproducibility capsule for research articles on
using Flash-X simulations to develop scientific machine learning models for
thermal science applications.

Figure 1: (a) A common scenario in computational sciences where a developer
has to conduct individual experiments and collaborate with other researchers
who may not be experts in using their software. (b) Example use case of lab
notebooks for Flash-X development. The developer maintains a personal
repository that serves as a seed for setting up experiments for physics and
performance related studies that can be shared with collaborators.

Figure 2 depicts the directory tree that might be set up to record this set of
experiments. Here, the first subdirectory, software records all the software
configuration information that is not directly available from the log file. In
this example there is an external dependency AMReX [5], a library that
provides adaptive mesh refinement capabilities. The build specifications of
external libraries are not known to Flash-X configurator, and must be
explicitly recorded to be able to exactly reproduce the execution environment.
The second subdirectory simulation describes the application instance that is
used in the experiment. If an experiment uses more than one application
instance it is expected that a separate subdirectory would be created for
each. Because Flash-X has a high degree of composability, and therefore a huge
collection of tests, it is also important to list the tests that were used to
ensure ongoing correctness of the code throughout the duration of the
experiment. Organizing the tests within this structure also enables quick
running of the test-suite as described later.

.1 Experiment. .2 software/. .3 amrex/. .4 setupAMReX.sh. .4 Jobfile. .3
flashx/. .4 setupFlashX.sh. .4 Jobfile. .2 simulation/. .3 FlowBoiling/. .4
flashBuild.sh. .4 flashRun.sh. .4 Jobfile. .4 jobnode.archive/. .5 <mm-dd-
yyyy>/. .2 tests/. .3 Tests.suite. .3 runTests.sh. .3 Jobfile. .2
environment.sh. .2 Jobfile.

Figure 2: Directory tree for a Flash-X experiment.

The process of seeding and managing these notebooks is implemented using a
lightweight command line tool, Jobrunner [6], that enforces a directory-based
inheritance for application configurations encoded in shell scripts to manage
computational experiments. The directory tree is designed to split different
shell commands that set up the software stack and simulation in an optimal
way. For instance, a combination of environment.sh and setupAMReX.sh is used
to configure AMReX most suitably for the experiment. Here, environment.sh
documents and sets variables that store compiler information and installation
location of different libraries, and setupAMReX.sh stores information specific
to building AMReX with the necessary options. The Jobfiles that exist at every
level of the tree encode information that Jobrunner can use to stitch together
these bash files to perform different tasks.

⬇

# Location:

# Experiment

job:

setup:

- environment.sh

submit:

- environment.sh

\par# Location:

# Experiment/software/amrex

job:

setup:

- setupAMReX.sh

\par# Location:

# Experiment/simulation/FlowBoiling

job:

setup:

- flashSetup.sh

submit:

- flashSubmit.sh

archive:

- "*hdf5*"

- "*.log"

Figure 3: Contents of Jobfiles at different locations along the directory
tree.

Figure 3 shows contents of the Jobfiles at different nodes along the directory
tree described in Figure 2. During the execution of Jobrunner commands, files
assigned to the respective commands are picked up to execute tasks in the
target node of the directory tree. For instance, the commands described in
Figure 4 build AMReX, clone and configure an application instance of Flash-X
for multiphase pool boiling problem, compile the assembled code, and then
execute the created binary to obtain results. At the same time, local testing
of the code is performed using Flash-X’s custom testing
framework555https://github.com/Flash-X/Flash-X-Test which uses tests encoded
in Tests.suite along with the environment.sh and runTests.sh scripts to build
and execute tests. Note that environment.sh, located at the root of the
project directory is used for each Jobrunner command, providing consistency
between tests and experiments.

Jobrunner hides error-prone individual steps of the experiments while allowing
organized and explicit documentation of configuration options that can be
easily modified by editing the shell scripts. The directory tree in Figure 2
can be easily modified and redesigned based on the requirements of the
experiments to include markdown notes and analysis files. For full
functionality and documentation for Jobrunner see A. Dhruv [6].

⬇

# Setting up software and dependencies

jobrunner setup software/amrex

jobrunner setup software/flashx

\par# Setting up and running experiments

jobrunner setup simulation/FlowBoiling

jobrunner submit simulation/FlowBoiling

\par# Archiving results to jobnode.archive

jobrunner archive \ simulation/FlowBoiling

\par# Running Flash-X test suite

jobrunner submit tests

Figure 4: Jobrunner commands for setting up dependencies, running tests and
experiments, and archive data. These commands are executed from the root of
the directory-tree

Data archiving is the last remaining concern for reproducibility. This is
implemented using Jobrunner’s archive command which picks up file patterns
listed in Jobfiles (see Figure 3) and moves them to a jobnode.archive/<mm-dd-
yyyy> directory under the target node of an experiment (see
Experiment/simulation/FlowBoiling in Figure 2). The directories containing the
data are eventually moved to a cloud-based archival service along with a clone
of the Github source repository to preserve the tree structure. Note that the
raw data is itself not included in the repository because that can be quite
large. It is organized in such a way that the archive can be unpacked below
the directory structure maintained in the repository to exactly reconstruct
all the artifacts of the experiment. See Lab-Notebooks/Outflow-Forcing-
BubbleML [3] for an example lab notebook for recent publications.

## CONCLUSION

An increasing emphasis on reproducibility and greater scrutiny of
computational science results is slowly changing the perception of what
constitutes a good computational experiment. The concept of maintaining
laboratory notebooks in computational sciences has been gaining popularity to
enforce structure and rigor in scientific studies. This is a welcome change,
and deserves encouragement. However, unlike more traditional scientific
disciplines, computational scientists have challenges in how frequently their
execution environments may change, and the kind of impact they may have on the
continuity of their experiments. Additionally, they often participate in teams
that may be geographically diverse, therefore they need distributed digital
mechanisms to record their work. Tools such as Jupyter notebook and Github
repositories can help, but need additional care to fully capture the
provenance of a computational experiment. We have presented an approach that
addresses many of the challenges faced by software that is being developed
while also being used for production. It is our hope that tools and ideas
presented here serve as a motivation for other scientists to design and
organize their experiments. Integration of well-organized laboratory notebooks
with reproducibility and data capsules can improve the quality of scientific
artifacts, and enhance productivity of collaborative research.

## ACKNOWLEDGMENTS

The authors would like to acknowledge Jared O’Neal’s work in promoting the
culture of using lab notebooks for Flash-X experiments, and for
conceptualizing the format that has been adopted by the team and described in
this article.

This work was partially supported by the Exascale Computing Project
(17-SC-20-SC), a collaborative effort of the US Department of Energy Office of
Science and the National Nuclear Security Administration, and the Laboratory
Directed Research and Development Program supported by Argonne. The submitted
manuscript was created by UChicago Argonne, LLC, operator of Argonne National
Laboratory (“Argonne”). Argonne, a U.S. Department of Energy Office of Science
laboratory, is operated under Contract No. DE-AC02-06CH11357. The U.S.
Government retains for itself, and others acting on its behalf, a paid-up
nonexclusive, irrevocable worldwide license in said article to reproduce,
prepare derivative works, distribute copies to the public, and perform
publicly and display publicly, by or on behalf of the Government. The
Department of Energy will provide public access to these results of federally
sponsored research in accordance with the DOE Public Access Plan.
http://energy.gov/downloads/doe-public-access-plan.

## References

  * [1] A. Dubey, K. Weide, J. O’Neal, A. Dhruv, S. Couch, J. A. Harris, T. Klosterman, R. Jain, J. Rudi, B. Messer, M. Pajkos, J. Carlson, R. Chu, M. Wahib, S. Chawdhary, P. M. Ricker, D. Lee, K. Antypas, K. M. Riley, C. Daley, M. Ganapathy, F. X. Timmes, D. M. Townsley, M. Vanella, J. Bachan, P. M. Rich, S. Kumar, E. Endeve, W. R. Hix, A. Mezzacappa, and T. Papatheodore.  Flash-X: A multiphysics simulation software instrument.  SoftwareX, 19:101168, 2022.
  * [2] A. Dubey, K. Antypas, M. K. Ganapathy, L. B. Reid, K. Riley, D. Sheeler, A. Siegel, and K. Weide.  Extensible component-based architecture for flash, a massively parallel, multiphysics simulation code.  Parallel Computing, 35(10):512–522, 2009.
  * [3] A. Dhruv.  Reproducibility capsules and execution environments for computational experiments.  https://github.com/Lab-Notebooks, 2023.
  * [4] A. Dhruv.  Execution environment for Flash-X software development.  https://github.com/akashdhruv/Multiphysics-Simulations, 2023.
  * [5] W. Zhang, A. Almgren, V. Beckner, J. Bell, J. Blaschke, C. Chan, M. Day, B. Friesen, K. Gott, D. Graves, M. Katz, A. Myers, T. Nguyen, A. Nonaka, M. Rosso, S. Williams, and M. Zingale.  AMReX: a framework for block-structured adaptive mesh refinement.  Journal of Open Source Software, 4(37):1370, May 2019.
  * [6] A. Dhruv.  Jobrunner: A command line tool to manage execution environments for computational studies.  Zenodo, 10.5281/zenodo.8322955, September 2023.  https://github.com/akashdhruv/Jobrunner

Akash Dhruv, is a Postdoctoral Fellow in the Mathematics and Computer Science
Division at Argonne National Laboratory where he performs research in DevOps
and performance engineering for scientific computing applications, along with
building computational pipelines to integrate simulation with machine learning
workflows. He received his Ph.D. in mechanical and aerospace engineering from
George Washington University, Washington D.C., and a B.Tech in mechanical
engineering from National Institute of Technology, Surat.

Anshu Dubey, is a Senior Computational Scientist in the Mathematics and
Computer Science Division at Argonne National Laboratory and a Senior
Scientist (CASE) at the University of Chicago. She was previously at Lawrence
Berkeley National Laboratory, and before that she was the associate director
of the Flash Center for Computational Science at the University of Chicago and
also led the CS/Applications group, who develop, maintain, and distribute the
FLASH code. She received her Ph.D. in computer science from Old Dominion
University and a B.Tech. in electrical engineering from the Indian Institute
of Technology, New Delhi.


0cm section[0.5cm] []0.5cm [] subsection[1.25cm] []0.75cm []
subsubsection[2.25cm] []1.0cm []

# Extending normality: A case of unit distribution generated from the moments
of the standard normal distribution

Miguel S. Concha-Aracena<EMAIL_ADDRESS>Departamento de
Matemática, Universidad de Atacama, Copiapó, Chile Leonardo Barrios-Blanco
<EMAIL_ADDRESS>Departamento de Matemática, Universidad
de Atacama, Copiapó, Chile David Elal-Olivero<EMAIL_ADDRESS>Departamento
de Matemática, Universidad de Atacama, Copiapó, Chile Paulo Henrique Silva
<EMAIL_ADDRESS>Department of Statistics, Federal University of Bahia,
Salvador, Brazil Diego Carvalho do Nascimento<EMAIL_ADDRESS>Departamento de Matemática, Universidad de Atacama, Copiapó, Chile

###### Abstract

This article presents an important theorem, which shows that from the moments
of the standard normal distribution one can generate density functions
originating a family of models. Additionally, we discussed that different
random variable domains are achieved with transformations. For instance, we
adopted the moment of order two, from the proposed theorem, and transformed
it, which allowed us to exemplify this class as unit distribution. We named it
as Alpha-Unit (AU) distribution, which contains a single positive parameter
$\alpha$ ($\text{AU}(\alpha)\in[0,1]$). We presented its properties and showed
two estimation methods for the $\alpha$ parameter, the maximum likelihood
estimator (MLE) and uniformly minimum-variance unbiased estimator (UMVUE)
methods. In order to analyze the statistical consistency of the estimators, a
Monte Carlo simulation study was carried out, where the robustness was
demonstrated. As real-world application, we adopted two sets of unit data, the
first regarding the dynamics of Chilean inflation in the post-military period,
and the other regarding the daily maximum relative humidity of the air in the
Atacama Desert. In both cases shown, the AU model is competitive, whenever the
data present a range greater than 0.4 and extremely heavy asymmetric tail. We
compared our model against other commonly used unit models, such as the beta,
Kumaraswamy, logit-normal, simplex, unit-half-normal, and unit-Lindley
distributions.
Keywords: Asymmetry accommodation; rates and proportions; single-parameter
distribution; unit distribution; water monitoring;

##  Introduction

Statistical methodology plays an important role in quantitative methods, given
the hypothesis testing and inferential procedures. Nonetheless, the comparison
across features is given based on a generated function estimated from the data
information. Most often, mild suppositions are taken compromising the
generalization of the results.

Under the perspective of statistical generalization (inferential method), some
challenges are found for bounded distribution estimation. For instance, the
confidence interval, which is often adopted from the maximum likelihood
estimation approach and asymptotic supposition, is also assumed. Specially,
interval estimation can be seen off the parameter space domain.

One exemplification is the case where bounded information data are observed,
nonetheless, normality is commonly assumed to be true. This is the case of
proportion/rate data, which are double bounded in the lower limit equal to 0
and upper limit equal to 1. Relative humidity is an example of this scenario
where every decision-making should be $\in[0,1]$ [fonseca2021water,
bayer2018beta], or commonly rates used in the field of finance, economics and
demography, to list a few.

In the case of rates and proportions processes, as well as other processes
whose variable of interest assumes values in the range $(0,1)$, there is a
well-represented class of models, the unit distributions family, which deals
with this type of double-bounded data. Among the existing unit distributions,
we can cite the power distribution, beta distribution [ferrari2004beta],
Kumaraswamy distribution [kumaraswamy1980generalized], unit-logistic
distribution [tadikamalla1982systems], simplex distribution
[barndorff1991some], unit-Weibull distribution [mazucheli2018unitw,
mazucheli2020unitw], unit-Lindley distribution [mazucheli2019ul], unit-half
normal distribution [bakouch2021flexible], unit log-log distribution
[korkmaz2021unitll], modified Kumaraswamy and reflected modified Kumaraswamy
distributions [sagrillo2021modified], unit-Teissier distribution
[krishna2022unitt], unit extended Weibull families of distributions
[guerra2021unit], unit folded normal distribution [korkmaz2022unit], unit-Chen
distribution [korkmaz2022unitc], and Marshall-Olkin reduced Kies distribution
[afify2022new].

Despite the applicability of the unit distributions in double-bounded
variables, another important fact is that the interval estimation for the
parameter may also be limited in a domain (like positive real number). In this
manner, we also presented an inferential alternative through the delta method.

This work starts by presenting an important theorem that transforms from a
modification of the standard normal distribution into a class of density
distributions that can be seen as unit. Then, as an exemplification, a case of
second moment was chosen to illustrate the usefulness of this class of
probabilistic models. This class of distributions shows to be competitive for
high-frequency data with range greater than 0.4, important to real-world
applications, whereas classical unit distribution fails [santana2022beta].
Additionally, two different data sets were selected to illustrate the
adjustment of the proposed model. The first is related to Chilean inflation
(ultimate post-military era), and the second is from the dryest area of the
planet (excluding the north and south poles).

###  Motivation

The normal distribution is very important in the history of statistics, where
numerous modifications to this distribution have been proposed in the
literature [stahl2006evolution, limpert2011problems]. An interesting fact
related to the normal distribution is that its even moments can be used to
generate new distributions, as is the case that we will show below, through a
definition and a result embodied in a theorem that accounts for the
characterization of these new distributions.

Definition 1. A random variable $B$ is said to be distributed according to a
Bimodal Normal (BN) distribution of order $k$, that is, $B\sim\text{BN}(k)$,
discussed in [elal2010alpha], if its probability density function (PDF) is
given by

$f(b\mid k)=\frac{1}{c}b^{2k}\phi(b),\qquad b\in\mathbb{R},$ (1.1)

where $\phi(\cdot)$ is the PDF of a standard normal distribution,
$c=\prod_{j=1}^{k}(2j-1)$ and $k=\\{1,2,3,\ldots\\}$.

This class of distributions is always bimodal, where the observed modes move
away when the order $k$ increases (see Figure 1).

Figure 1: Density function of the BN distribution by varying the parameter $k$
(displayed on the top of each chart).

It is interesting to mention that transformations derived from the
$\text{BN}(k)$ distribution may lead to other domains of interest, e.g., the
unit domain. For example, let $B\sim\text{BN}(k)$, then by adding a scale
parameter $\alpha$, the transformation $\alpha|B|\in\mathbb{R}^{+}$, and then
the transformation $e^{-\alpha|B|}\in[0,1]$. Therefore, the stochastic
characterization of a $\text{BN}(k)$ distribution can be obtained according to
the following theorem.

Theorem 1. Let $W_{1}$ and $W_{2}$ be independent random variables, where
$W_{1}$ is such that $\mathbb{P}(W_{1}=1)=\mathbb{P}(W_{1}=-1)=1/2$ and
$W_{2}\sim\chi^{2}_{2k+1}$. Then,

$W_{1}\sqrt{W_{2}}\sim\text{BN}(k).$ (1.2)

So, this theorem is mainly motivated by the result that shows that if
$X\sim\text{BN}(k)$, then $X^{2}\sim\chi^{2}_{2k+1}$. The demonstration is
presented in Appendix A.

##  Distribution of the Second Moment of the Unit-Normal Distribution

In this section, we will discuss a new unit distribution, named Alpha-Unit,
which presents a single parameter, $\alpha$. Whereas it will be presented its
stochastic representations (probability density and cumulative distribution
functions), moments, characteristic function, and how to generate random
numbers from it.

By taking the general theorem presented, and considering $k=1$, that is,
considering the second moment of the standard normal distribution and its
transform, a new unit distribution called Alpha-Unit will be illustrated.
However, as $k$ increases, the concentration of the distribution intensifies.

###  Properties and Characterization

Definition 2. (Alpha-Unit distribution). A random variable $X$ follows an
Alpha-Unit (AU) distribution with parameter $\alpha>0$, that is,
$X\sim\text{AU}(\alpha)$, if its PDF is given by

$f_{X}(x\mid\alpha)=\frac{2}{x\alpha}\left(\frac{\ln(x)}{\alpha}\right)^{2}\phi\left(\frac{\ln(x)}{\alpha}\right),\qquad
0<x\leq 1.$ (2.1)

Remark 1. If $X\sim\text{AU}(\alpha)$, then its PDF is unimodal.

Demonstration. The maxima of the AU distribution are studied, for which the
criterion of the first derivative is first considered:

$\frac{df_{X}(x\mid\alpha)}{dx}=\frac{2}{x\alpha^{2}}\frac{\ln(x)}{\alpha}\phi\left(\frac{\ln(x)}{\alpha}\right)\left[\frac{2}{x}-\frac{\ln(x)}{x}-\frac{[\ln(x)]^{2}}{\alpha}\frac{1}{x\alpha}\right]=0.$

Solving algebraically for $x$, we obtain:

$x=\begin{cases}e^{\left(\frac{-\alpha^{2}+\sqrt{\alpha^{4}+8\alpha^{2}}}{2}\right)}&\mbox{(i)
}\\\ \\\
e^{\left(-\frac{\alpha^{2}+\sqrt{\alpha^{4}+8\alpha^{2}}}{2}\right)}&\mbox{(ii)
}\end{cases}.$

To see if either or both expressions are solutions, it must be true that

$e^{-h}\ \ \&\ \ h>0.$

By working algebraically, it can be seen that this is only true for (i),
therefore, the AU distribution is unimodal. $\blacksquare$

Proposition 1. If $X\sim\text{AU}(\alpha)$, then its $r$-th order moment is
given by

$\mathbb{E}[X^{r}]=2e^{\left(\frac{r^{2}\alpha^{2}}{2}\right)}\left[\left(1+r^{2}\alpha^{2}\right)\left(1-\Phi(r\alpha)\right)-r\alpha\phi(r\alpha)\right].$
(2.2)

Demonstration.

$\mathbb{E}[X^{r}]=\int_{0}^{1}x^{r}\frac{2}{x\alpha}\left(\frac{\ln(x)}{\alpha}\right)^{2}\phi\left(\frac{\ln(x)}{\alpha}\right)dx.$

By making the change of variables:

$\begin{cases}u=\frac{1}{\alpha}\ln(x)\quad\Rightarrow\quad e^{u\alpha}=x\\\
\\\ du=\frac{1}{\alpha x}dx\quad\Rightarrow\quad\alpha
e^{u\alpha}du=dx\end{cases},$

then substituting into the previous equation and developing algebraically, we
have:

$\mathbb{E}[X^{r}]=2e^{\frac{\alpha^{2}r^{2}}{2}}\int_{-\infty}^{0}u^{2}\frac{1}{\sqrt{2\pi}}e^{-\frac{(u-\alpha
r)^{2}}{2}}du.$

Then, by making another change of variables: $h=u-\alpha r$, $dh=du$; and
replacing these expressions in the previous equation, we have:

$\displaystyle\mathbb{E}[X^{r}]$
$\displaystyle=2e^{\frac{\alpha^{2}r^{2}}{2}}\int_{-\infty}^{-\alpha
r}\left(h+\alpha r\right)^{2}\frac{1}{\sqrt{2\pi}}e^{-\frac{h^{2}}{2}}dh$
$\displaystyle=2e^{\frac{\alpha^{2}r^{2}}{2}}\int_{-\infty}^{-\alpha
r}\left(h^{2}+2h\alpha r+\alpha^{2}r^{2}\right)\phi(h)dh$
$\displaystyle=2e^{\frac{\alpha^{2}r^{2}}{2}}\left(\int_{-\infty}^{-\alpha
r}h^{2}\phi(h)dh+2\alpha r\int_{-\infty}^{-\alpha
r}h\phi(h)dh+\alpha^{2}r^{2}\int_{-\infty}^{-\alpha r}\phi(h)dh\right).$

Solving the integrals, we get:

$\mathbb{E}[X^{r}]=2e^{\frac{\alpha^{2}r^{2}}{2}}\left[\alpha r\phi(\alpha
r)+\left(1-\Phi(\alpha r)\right)-2\alpha r\phi(\alpha
r)+\alpha^{2}r^{2}\left(1-\Phi(\alpha r)\right)\right].$

Then, solving algebraically, we arrive at Proposition 1. $\blacksquare$

From Proposition 1, we obtain the mean and variance of the $\text{AU}(\alpha)$
model as follows:

$\displaystyle\mathbb{E}[X]$
$\displaystyle=2e^{\frac{\alpha^{2}}{2}}\left[(1+\alpha^{2})\left(1-\Phi(\alpha)\right)-\alpha\phi(\alpha)\right],$
$\displaystyle\mathbb{V}\text{ar}[X]$
$\displaystyle=\mathbb{E}[X^{2}]-\left(\mathbb{E}[X]\right)^{2}$
$\displaystyle=2e^{2\alpha^{2}}\left[(1+4\alpha^{2})\left(1-\Phi(2\alpha)\right)-2\alpha\phi(2\alpha)\right]-4e^{\alpha^{2}}\left[(1+\alpha^{2})\left(1-\Phi(\alpha)\right)-\alpha\phi(\alpha)\right]^{2},$

where $\Phi(\cdot)$ is the cumulative distribution function (CDF) of a
standard normal distribution.

Remark 2. As illustration, Figure 2 shows the generated asymmetry and kurtosis
based on the chosen $\alpha$ parameter of the AU distribution.

Figure 2: Density function of the AU distribution by varying the parameter
$\alpha$ (displayed on the top of each chart). Whereas
$B\sim\text{BN}(1)\rightarrow B^{2}\sim\chi^{2}_{3}$, then the AU model was
generated from $X=e^{-\alpha|B|}$.

Proposition 2. If $X\sim\text{AU}(\alpha)$, then its CDF is given by

$F_{X}(x\mid\alpha)=2\Phi\left(\frac{\ln(x)}{\alpha}\right)-2\left(\frac{\ln(x)}{\alpha}\right)\phi\left(\frac{\ln(x)}{\alpha}\right).$
(2.3)

Demonstration. By definition of the CDF, we have:

$F_{X}(x\mid\alpha)=\int_{0}^{x}\frac{2}{t\alpha}\left(\frac{\ln(t)}{\alpha}\right)^{2}\phi\left(\frac{\ln(t)}{\alpha}\right)dt.$

By making the change of variables:

$\begin{cases}u=\frac{\ln(t)}{\alpha}\quad\Rightarrow\quad e^{u\alpha}=t\\\
\\\ du=\frac{1}{\alpha t}dt\quad\Rightarrow\quad\alpha
e^{u\alpha}du=dt\end{cases},$

then substituting into the previous equation and reducing expressions
algebraically, we get:

$F_{X}(x\mid\alpha)=2\int_{-\infty}^{\frac{\ln(x)}{\alpha}}u^{2}\phi(u)du.$

Calculating the integral, we have:

$\displaystyle F_{X}(x\mid\alpha)$
$\displaystyle=2\left[-u\phi(u)\Big{|}_{-\infty}^{\ln(x)/\alpha}+\int_{-\infty}^{\ln(x)/\alpha}\phi(u)du\right]$
$\displaystyle=2\left[-\left(\frac{\ln(x)}{\alpha}\right)\phi\left(\frac{\ln(x)}{\alpha}\right)+\Phi\left(\frac{\ln(x)}{\alpha}\right)\right].$

Multiplying and commuting, we arrive at the expression of Proposition 2.
$\blacksquare$

Additionally, if $X$ denotes the monitored variable, then the PDF of $X$ is
given by (2.1). Also, consider that the probability of false alarm (known as
type I error) is $\pi$. Thus, we have:

$\mathbb{P}(X<\text{LCL}\mid\alpha)=\mathbb{P}(X>\text{UCL}\mid\alpha)=\pi/2,$

where $\alpha$ is the in-control process parameter (that is, the parameter
that controls the quality characteristic based on the in-control state), LCL
and UCL are the lower and upper control chart limits, respectively. Given the
CDF $F_{X}(x\mid\alpha)$, then the quantile function of $X$ is defined by
$Q(p\mid\alpha)=F_{X}^{-1}(p\mid\alpha)$, and can be obtained by setting to 0
and solving (numerically) for $x$ the following equation:

$\Phi\left(\frac{\ln(x)}{\alpha}\right)-\left(\frac{\ln(x)}{\alpha}\right)\phi\left(\frac{\ln(x)}{\alpha}\right)-\frac{p}{2},\qquad\text{for}\,\,\,0<p<1.$
(2.4)

Following [bayer2018CIE], the control limits and centerline (CL) of the
proposed AU control chart are given by

$\text{LCL}=Q\left(\pi/2\mid\alpha\right),\qquad\text{CL}=\mathbb{E}[X\mid\alpha],\\\
\qquad\text{UCL}=Q\left(1-\pi/2\mid\alpha\right),$ (2.5)

where $Q(.)$ is the quantile function of the $\text{AU}(\alpha)$ distribution.

Proposition 3. If $X\sim\text{AU}(\alpha)$, then its moment-generating
function (MGF) is given by

$\psi_{X}(t\mid\alpha)=2\sum_{k=0}^{\infty}\frac{t^{k}}{k!}e^{\left(\frac{k^{2}\alpha^{2}}{2}\right)}\left[\left(1+k^{2}\alpha^{2}\right)\left(1-\Phi(k\alpha)\right)-k\alpha\phi(k\alpha)\right].$
(2.6)

Demonstration. By definition of the MGF, we have:

$\psi_{X}(t\mid\alpha)=\mathbb{E}\left[e^{tx}\right]=\int_{0}^{1}e^{tx}\frac{2}{x\alpha}\left(\frac{\ln(x)}{\alpha}\right)^{2}\phi\left(\frac{\ln(x)}{\alpha}\right)dx.$

By making the following change of variables:

$\begin{cases}u=\frac{\ln(x)}{\alpha}\quad\Rightarrow\quad e^{u\alpha}=x\\\
\\\ du=\frac{1}{\alpha x}dx\quad\Rightarrow\quad\alpha
e^{u\alpha}du=dx\end{cases},$

then substituting and simplifying into the previous equation, we get:

$\displaystyle\psi_{X}(t\mid\alpha)$
$\displaystyle=2\int_{-\infty}^{0}e^{\left(te^{u\alpha}\right)}u^{2}\phi(u)du$
$\displaystyle=2\int_{-\infty}^{0}\sum_{k=0}^{\infty}\frac{t^{k}e^{u\alpha
k}}{k!}u^{2}\phi(u)du.$

Working algebraically, we have:

$\psi_{X}(t\mid\alpha)=2\sum_{k=0}^{\infty}\frac{t^{k}}{k!}e^{\left(\frac{\alpha^{2}k^{2}}{2}\right)}\int_{-\infty}^{0}u^{2}\frac{1}{\sqrt{2\pi}}e^{\left(\frac{-(u-\alpha
k)^{2}}{2}\right)}du.$

By making the following change of variables: $h=u-\alpha k$, then substituting
into the previous equation, we have:

$\psi_{X}(t\mid\alpha)=2\sum_{k=0}^{\infty}\frac{t^{k}}{k!}e^{\left(\frac{\alpha^{2}k^{2}}{2}\right)}\int_{-\infty}^{-\alpha
k}(h+\alpha k)^{2}\phi(h)dh.$

Then, solving the integral and adjusting algebraically, we arrive at the
expression of Proposition 3. $\blacksquare$

Algorithm 1 Random number generation from the $\text{AU}(\alpha)$ model.

Step 1. Generate a random number $X_{1}\sim\chi^{2}_{3}$.

Step 2. Generate a random number $U\sim\text{Uniform}(0,1)$. If $U\leq 1/2$,
set $V=\sqrt{X_{1}}$; otherwise, $V=-\sqrt{X_{1}}$.

Step 3. Based on the numbers obtained, generate the variable $Y=\alpha|V|$,
where $\alpha$ is a (positive) scale parameter and $|V|$ is a Bimodal Half-
Normal (BHN).

Step 4. Conclude with the number generated by Step 2 as a negative power of
base $e$, that is, $X=e^{-Y}=e^{-\alpha|V|}\in[0,1]$.

Step 5. Repeat Steps 1-4 $n$ times to obtain a random sample of size $n$ from
the $\text{AU}(\alpha)$ model.

The pseudo-code describes the important steps for the generation of random
numbers from the AU distribution. Further proofs are attached under Appendix
B.

##  Inference

In this section, we will discuss the parameter estimation adopting the
uniformly minimum-variance unbiased estimator (UMVUE) and maximum likelihood
estimator (MLE) approaches. First, it will be demonstrated that the UMVUE is
obtained straightforward since the proposed AU distribution is part of the
exponential family. Later, the MLE will be also discussed, in which it will
help to estimate not only with the point estimation of the $\alpha$ parameter,
but also with the interval estimation. We enrolled the reasoning considering
the asymptotic convergence in distribution of the parameter estimator, as well
as we adapted a transformation which ensures that the interval of the
parameter will always be on its domain (the delta method). The delta
transformation procedure will enable the correct inferences and the standard
error calculation associated with the parameter estimate. Later on, we will
present a simulation study to illustrate these theoretical results.

###  UMVUE through the Exponential Family

Many of the distributions used in statistics belong to the exponential family,
thereby implying a considerable advantage over other models that do not belong
to this family. Such an advantage is declared significantly when it comes to
calculating the statistic $T(\bm{X})$ of a random sample
$\bm{X}=\left(X_{1},X_{2},\ldots,X_{n}\right)^{\top}$. We will show below that
the proposed $\text{AU}(\alpha)$ model belongs to this family.

A random variable $X$ is said to belong to a one-parameter exponential family
if its associated PDF $f(\cdot\mid\theta)$ can be written in the form

$f(x\mid\theta)=\exp\left\\{c(\theta)T(x)+d(\theta)+S(x)\right\\}.$

Let $X\sim\text{AU}(\alpha)$, then the PDF of $X$ can be written in
exponential form as follows:

$f(x\mid\alpha)=\exp\left\\{-\frac{1}{2\alpha^{2}}[\ln(x)]^{2}-3\ln(\alpha)+\ln\left(\frac{[\ln(x)]^{2}}{x\sqrt{2\pi}}\right)\right\\}.$

Then, $X$ belongs to a one-parameter exponential family if we define

$c(\alpha)=-\frac{1}{2\alpha^{2}},\quad T(x)=[\ln(x)]^{2},\quad
d(\alpha)=-3\ln(\alpha),\quad\text{and}\quad
S(x)=\ln\left(\frac{[\ln(x)]^{2}}{x\sqrt{2\pi}}\right).$

Let $\bm{x}=\left(x_{1},x_{2},\ldots,x_{n}\right)^{\top}$ be an observation
from the random sample $\bm{X}=\left(X_{1},X_{2},\ldots,X_{n}\right)^{\top}$,
with $X_{i}\sim\text{AU}(\alpha)$, for $i=1,2,\ldots,n$, then the joint PDF
presented in exponential form would be given by

$f(\bm{x}\mid\alpha)=\exp\left\\{-\frac{1}{2\alpha^{2}}\sum_{i=1}^{n}[\ln(x_{i})]^{2}-3n\ln(\alpha)+\sum_{i=1}^{n}\ln\left(\frac{[\ln(x_{i})]^{2}}{x_{i}\sqrt{2\pi}}\right)\right\\},$

from which it can be concluded that the statistic
$T(\bm{X})=\sum_{i=1}^{n}[\ln(X_{i})]^{2}$ is sufficient and complete, once
the AU model is part of the exponential family.

Proposition 4. Let $\bm{X}=\left(X_{1},X_{2},\ldots,X_{n}\right)^{\top}$ be a
random sample with $X_{i}\sim\text{AU}(\alpha)$, for $i=1,2,\ldots,n$, and
$T(\bm{X})=\sum_{i=1}^{n}[\ln(X_{i})]^{2}$, then

$W_{n}=\frac{1}{\alpha^{2}}T(\bm{X})\sim\chi_{3n}^{2}.$

Demonstration. If $D=\left[\frac{\ln(X)}{\alpha}\right]^{2}$, then
$D\sim\chi^{2}_{3}$. Thus, $n$ independent and identically distributed samples
of $D$ will have the sum of $n$ $\chi^{2}_{3}$, which will result in a
$\chi^{2}$ distribution with degrees of freedom equal to $3n$, that is,
$\chi^{2}_{3n}$, since

$\displaystyle F_{D}(d)$ $\displaystyle=\mathbb{P}(D\leq
d)=\mathbb{P}\left(\left[\frac{\ln(X)}{\alpha}\right]^{2}\leq
d\right)=\mathbb{P}\left(-\sqrt{d}\leq\frac{\ln(X)}{\alpha}\leq\sqrt{d}\right)$
$\displaystyle=\mathbb{P}\left(-\alpha\sqrt{d}\leq\ln(X)\leq\alpha\sqrt{d}\right)=\mathbb{P}\left(\ln(X)\leq\alpha\sqrt{d}\right)-\mathbb{P}\left(\ln(X)\leq-\alpha\sqrt{d}\right)$
$\displaystyle=1-\mathbb{P}\left(\ln(X)\leq-\alpha\sqrt{d}\right)=1-\mathbb{P}\left(X\leq
e^{-\alpha\sqrt{d}}\right)=1-F_{X}\left(e^{-\alpha\sqrt{d}}\right),$

so,

$\displaystyle f_{D}(d)$ $\displaystyle=\frac{\partial F_{D}(d)}{\partial
d}=f_{X}\left(e^{-\alpha\sqrt{d}}\right)\left(e^{-\alpha\sqrt{d}}\right)\left(\frac{\alpha}{2\sqrt{d}}\right)$
$\displaystyle=\frac{2}{\alpha
e^{-\alpha\sqrt{d}}}\left(\frac{-\alpha\sqrt{d}}{\alpha}\right)^{2}\phi\left(\frac{-\alpha\sqrt{d}}{\alpha}\right)e^{-\alpha\sqrt{d}}\frac{\alpha}{2\sqrt{d}}$
$\displaystyle=\frac{1}{\sqrt{d}}\left(\sqrt{d}\right)^{2}\frac{1}{\sqrt{2\pi}}e^{-\frac{(\sqrt{d})^{2}}{2}}=\frac{1}{\sqrt{2\pi}}d^{1/2}\exp\left(-d/2\right)\
\ \equiv\ \ \chi^{2}_{3}.\quad\blacksquare$

Proposition 5. Let $\bm{X}=\left(X_{1},X_{2},\ldots,X_{n}\right)^{\top}$ be a
random sample with $X_{i}\sim\text{AU}(\alpha)$, for $i=1,2,\ldots,n$, and
$T(\bm{X})=\sum_{i=1}^{n}[\ln(X_{i})]^{2}$, then

$S(\bm{X})=\frac{\Gamma\left(\frac{3n}{2}\right)\sqrt{2}}{\Gamma\left(\frac{3n+1}{2}\right)}\sqrt{T(\bm{X})}$

is an unbiased estimator of $\alpha$.

Demonstration. Remembering that if $X\sim\text{Gamma}(a,b)$ distribution, then
$\mathbb{E}[X^{k}]=\frac{\Gamma(a+b)}{b^{k}\Gamma(a)}$. Once the parameter
$\alpha$ is observed to be squared, it will be necessary to release it to find
an unbiased estimator. So, considering the random variable $W_{n}^{1/2}$ (with
$W_{n}$ as defined in Proposition 4), then

$\mathbb{E}\left[(W_{n})^{1/2}\right]=\frac{\Gamma\left(\frac{3n}{2}+\frac{1}{2}\right)}{2^{1/2}\Gamma\left(\frac{3n}{2}\right)},$

so,

$\displaystyle\mathbb{E}\left[\left(\frac{1}{\alpha^{2}}T(\bm{X})\right)^{1/2}\right]$
$\displaystyle=\frac{\Gamma\left(\frac{3n}{2}+\frac{1}{2}\right)}{2^{1/2}\Gamma\left(\frac{3n}{2}\right)}$
$\displaystyle\mathbb{E}\underbrace{\left[\sqrt{T(\bm{X})}\frac{\Gamma\left(\frac{3n}{2}\right)\sqrt{2}}{\Gamma\left(\frac{3n}{2}+\frac{1}{2}\right)}\right]}_{S(\bm{X})}$
$\displaystyle=\alpha.\quad\blacksquare$

Remark 3. Considering the two previous propositions and resorting to the
Lehmann-Scheffé theorem, we can conclude that $S(\bm{X})$ is UMVUE for
$\alpha$.

###  Estimation using the Maximum Likelihood Method

Let $\bm{x}=\left(x_{1},x_{2},\ldots,x_{n}\right)^{\top}$ be a realization of
the random sample $\bm{X}=\left(X_{1},X_{2},\ldots,X_{n}\right)^{\top}$ taken
from the $\text{AU}(\alpha)$ model. Then, the log-likelihood function is given
by

$\ell(\alpha)=\text{constant}-3n\ln(\alpha)-\Sigma_{i=1}^{n}\ln(x_{i})+2\Sigma_{i=1}^{n}\ln(\ln(x_{i}))-\frac{1}{2\alpha^{2}}\Sigma_{i=1}^{n}[\ln(x_{i})]^{2}.$

The MLE of $\alpha$, i.e., $\widehat{\alpha}$, is found by solving the
following equation:

$\frac{d\ell(\alpha)}{d\alpha}=-\frac{3n}{\alpha}+\frac{1}{\alpha^{3}}\Sigma_{i=1}^{n}[\ln(x_{i})]^{2}=0,$

giving

$\widehat{\alpha}=\left\\{\frac{1}{3n}\sum_{i=1}^{n}[\ln(x_{i})]^{2}\right\\}^{1/2}.$

On the other hand, the second derivative of $\ell(\alpha)$ evaluated at
$\alpha=\widehat{\alpha}$ is negative, concluding that $\widehat{\alpha}$ is
MLE for $\alpha$.

It is known that, under certain regularity conditions,

$\sqrt{n}\left(\widehat{\alpha}-\alpha\right)\xrightarrow{\text{D}}\text{N}\left(0,I^{-1}(\alpha)\right),$

where
$I(\alpha)=-\mathbb{E}\left[\frac{d^{2}\ell(\alpha)}{d\alpha^{2}}\right]=\frac{6n}{\alpha^{2}}$.

A two-sided $100(1-\pi)\%$ confidence interval for $\alpha$ can be calculated
by

$\left[\widehat{\alpha}-z_{1-\pi/2}\,\sqrt{\mathbb{V}\text{ar}\left[\widehat{\alpha}\right]},\
\widehat{\alpha}+z_{1-\pi/2}\,\sqrt{\mathbb{V}\text{ar}\left[\widehat{\alpha}\right]}\right],$
(3.1)

where $z_{q}$ is the $q$-th percentile of the standard normal distribution.
The variance of $\widehat{\alpha}$ can be approximated by the inverse of the
observed Fisher information, as

$\mathbb{V}\text{ar}\left[\widehat{\alpha}\right]=I^{-1}\left(\widehat{\alpha}\right)=\frac{\widehat{\alpha}^{2}}{6n}.$
(3.2)

Since $\alpha$ is a positive value and we cannot guarantee that the lower
limit of the interval (3.1) is positive, we resort to the delta method to
remedy such a situation. For this, we define the function
$g:[0,\infty)\rightarrow\mathbb{R}$ as $g(\alpha)=\ln(\alpha)$, and knowing
that

$\sqrt{n}\left(g(\widehat{\alpha})-g(\alpha)\right)\xrightarrow{\text{D}}\text{N}\left(0,I^{-1}\left(\alpha\right)\left[\frac{dg\left(\alpha\right)}{d\alpha}\right]^{2}\right),$

we can then obtain an approximate two-sided $100(1-\pi)\%$ confidence interval
for $\alpha$ by

$\left[\widehat{\alpha}\exp\left(-\frac{z_{1-\pi/2}}{\sqrt{6n}}\right),\
\frac{\widehat{\alpha}}{\exp\left(-\frac{z_{1-\pi/2}}{\sqrt{6n}}\right)}\right].$
(3.3)

###  Simulation Study

To illustrate the presented inferences for the estimation of the AU
distribution, in this subsection we compared (via simulation study) the MLE
versus the UMVUE. Moreover, we considered the scenarios where the parameter
$\alpha=\\{0.1,0.3,0.5,0.7,1.1,1.5\\}$, considering sample sizes
$n=\\{100,200,500\\}$, through the Monte Carlo method with N=1,000
repetitions. All this procedure took into account the random number generator
for the $\text{AU}(\alpha)$ distribution presented in Algorithm 1. All
analyses made in this work adopted the open-source R software [Rcran].

In order to compare the performance of the proposed estimators (MLE and
UMVUE), since the true parameter value is known, we adopted the performance
metrics bias and mean squared error (MSE), which are defined as follows:

$\text{Bias}(\alpha)=\frac{1}{\text{N}}\sum_{i=1}^{\text{N}}(\widehat{\alpha}_{i}-\alpha)\quad\text{and}\quad\text{MSE}(\alpha)=\frac{1}{\text{N}}\sum_{i=1}^{\text{N}}(\widehat{\alpha}_{i}-\alpha)^{2},$

where $\widehat{\alpha}_{i}$ is the estimate for $\alpha$ in the $i$-th
iteration (point estimation). Additionally, based on the asymptotic results
presented in this work, we also calculated the confidence interval (CI) length
adopting the delta method from Equation (3.3) (interval estimation). That is,
it analyzed the average of all the upper limits of the confidence interval, as
well as the average of all the lower limits, and then calculated their
difference.

Table 1 shows the obtained average estimates of the parameter $\alpha$,
varying the sample size $n$, as well as the corresponding bias, MSE and CI
length (this last only for MLE) results.

Table 1: Average estimates, bias and MSE for the MLE and UMVUE of the single parameter ($\alpha$) of the AU distribution, considering different sample sizes ($n$). |  | MLE |  |  | UMVUE
---|---|---|---|---|---
$n$ | $\alpha$ | Estimate | Bias | MSE | CI Length |  | Estimate | Bias | MSE
100 | 0.1 | 0.0998 | -0.0001 | 1.6930e-05 | 0.0160 |  | 0.0999 | -8.2264e-05 | 1.6165e-05
200 |  | 0.0999 | -9.8758e-05 | 8.7306e-06 | 0.0113 |  | 0.0999 | -5.7124e-05 | 8.7314e-06
500 |  | 0.0999 | -3.3400e-06 | 3.5542e-06 | 0.0071 |  | 0.1000 | 1.3327e-05 | 3.5555e-06
100 | 0.3 | 0.2996 | -0.0003 | 0.0002 | 0.0480 |  | 0.2999 | -8.0656e-05 | 0.0002
200 |  | 0.2997 | -0.0002 | 7.8575e-05 | 0.0339 |  | 0.2998 | -0.0001 | 7.8582e-05
500 |  | 0.2999 | -1.0020e-05 | 3.1987e-05 | 0.0214 |  | 0.3002 | 0.0002 | 3.0979e-05
100 | 0.5 | 0.4994 | -0.0005 | 0.0004 | 0.0800 |  | 0.4999 | -0.0001 | 0.0004
200 |  | 0.4997 | -0.0004 | 0.0002 | 0.0565 |  | 0.4997 | -0.0003 | 0.0002
500 |  | 0.4999 | -1.6700e-05 | 8.8855e-05 | 0.0357 |  | 0.5000 | 6.6637e-05 | 8.8888e-05
100 | 0.7 | 0.6992 | -0.0007 | 0.0008 | 0.1120 |  | 0.6998 | -0.0002 | 0.0008
200 |  | 0.6993 | -0.0006 | 0.0004 | 0.0791 |  | 0.6996 | -0.0004 | 0.0004
500 |  | 0.6999 | -2.3380e-05 | 0.0001 | 0.0501 |  | 0.7000 | 9.3291e-05 | 0.0001
100 | 1.1 | 1.0987 | -0.0012 | 0.0020 | 0.1760 |  | 1.0997 | -0.0003 | 0.0020
200 |  | 1.0989 | -0.0010 | 0.0010 | 0.1244 |  | 1.0994 | -0.0006 | 0.0010
500 |  | 1.0999 | -3.6741e-05 | 0.0004 | 0.0787 |  | 1.1001 | 0.0001 | 0.0004
100 | 1.5 | 1.4983 | -0.0016 | 0.0038 | 0.2400 |  | 1.4996 | -0.0004 | 0.0038
200 |  | 1.4985 | -0.0014 | 0.0019 | 0.1696 |  | 1.4991 | -0.0008 | 0.0019
500 |  | 1.4999 | -5.0101e-05 | 0.0008 | 0.1073 |  | 1.5002 | 0.0002 | 0.0007

The asymptotic convergence of the MLE, towards the robustness, was noticed as
the sample size increases. In addition, both MLE and UMVUE’s bias and MSE are
small and tend to decrease as $n$ gets larger. On the other hand, the CI
length also decreases as the sample size increases.

As a last summary, regarding the robustness of the estimators, it was taken
the difference between the MLE and UMVUE estimates, considering each different
sample size $n$, then the interquartile range (IQR) was calculated per sample
size group. That is,
$\text{IQR}^{(ni)}(\widehat{\alpha_{1}}^{(ni)}_{\text{MLE}}-\widehat{\alpha_{1}}^{(ni)}_{\text{UMVUE}},\ldots,\widehat{\alpha_{j}}^{(ni)}_{\text{MLE}}-\widehat{\alpha_{j}}^{(ni)}_{\text{UMVUE}})$,
where $ni=\\{100,200,500\\}$ and
$\alpha_{j}=\\{\alpha_{1}=0.1,\alpha_{2}=0.3,\ldots,\alpha_{6}=1.5\\}$. For
instance, the IQR for $n=100$ was $0.00053$, when for $n=200$ reduced to
$0.00025$, and $n=500$ resulted in $0.00012$. This shows, in summary, that
every time the sample size is large, the error range gets smaller, regardless
the value of the $\alpha$ parameter.

##  Real-World Exemplifications

In this section, we exemplified two applications adopting the AU distribution
with real-world problems. The first case is with respect to the dynamics of
Chilean inflation in the post-military dictatorship period. And the second
case is regarding the relative humidity of the air in the northern Chilean
city of Copiapó (Atacama region).

The Chilean inflation data are recorded annually, whose values considered
range from 1992 to 2021. These are based on the period after the military
dictatorship of 1973-1990. It was analyzed the dynamic of the inflation data
(in %), which were standardized by min-max transformation, resulting in a unit
response variable (value between 0 and 1). The years 1990 and 1991 were
excluded, since they are considered as a period of transition. Then, total
amount of observations were 30 years (from 1992 until 2021).

On the other hand, the relative air humidity data cover a period from February
2015 to October 2022, with a one-hour recording format (104,415 observations).
Then, this dataset was transformed into daily maximum observation (6,226
observations).

###  Chilean Inflation (Post-Military Era)

Figure 3 presents the dynamics of Chilean inflation in the post-military
dictatorship period, demonstrating stability between the years 1999 and 2008.
The right panel shows the time series for inflation, where time is measured in
years, from year 1 (1992) until year 30 (2021). The left panel shows the
accumulation of the values of the time series, where a predominant trend is
shown around 0.1 of the inflation rate.

Figure 3: Chilean inflation in the period 1992-2021 (post-military era). The
histogram on the left shows a skewness of the data. The dynamic is represented
in the right panel, where a disturbance (outlier) is noticed in the year 2008
(observation #17).

Once the empirical dynamic of these data was analyzed, the most common unit
distributions, presented in the statistical literature, were adjusted. The
upper panel of Figure 4 illustrates the histogram for the inflation data,
where it is compared with different adjusted densities based on the MLE: AU,
beta (BE), Kumaraswamy (KUM), logit-normal (LOGITNO), simplex (SIMPLEX), unit-
half-normal (UHN), and unit-Lindley (ULINDLEY). The lower panel of the same
figure shows the fitted CDFs superimposed to the empirical CDF (ECDF).

Figure 4: Estimated densities superimposed to the histogram (top-chart), and
estimated CDFs superimposed to the ECDF (bottom-chart) (Chilean inflation
data).

In order to quantify the performance of the adjusted models, we analyzed the
Akaike Information Criterion (AIC) [akaike1977], and the Bayesian (or Schwarz)
Information Criterion (BIC) [schwarz1978]. The obtained results (see Table 2)
show the AU model as the best-fitted model to this data set. In addition, it
is possible to infer about the average of the phenomenon, that is, the
expectation of the model AU($\widehat{\alpha}=1.2059$), resulting in
$\mathbb{E}[X_{\text{Inflation}}]=0.1948$. In other words, the average Chilean
inflation, post-military era, is 19.48%.

Table 2: Parameter estimates, AIC and BIC values (Chilean inflation data). S.E. $=$ standard error. Model | Parameter Estimate (S.E.) | AIC | BIC
---|---|---|---
$\text{AU}(\alpha)$ | $\widehat{\alpha}=1.205943\,(0.008079)$ | $-47.89$ | $-46.49$
$\text{BE}(\mu,\sigma)$ | $\widehat{\mu}=0.185857\,(0.000496)$ | $-44.58$ | $-41.78$
| $\widehat{\sigma}=0.314688\,(0.001304)$ |  |
$\text{KUM}(\mu,\sigma)$ | $\widehat{\mu}=1.370127\,(0.045522)$ | $-43.63$ | $-40.83$
| $\widehat{\sigma}=7.968427\,(7.750459)$ |  |
$\text{LOGITNO}(\mu,\sigma)$ | $\widehat{\mu}=0.150323\,(0.000457)$ | $-46.23$ | $-43.43$
| $\widehat{\sigma}=0.916938\,(0.014013)$ |  |
$\text{SIMPLEX}(\mu,\sigma)$ | $\widehat{\mu}=0.182462\,(0.000584)$ | $-43.17$ | $-40.37$
| $\widehat{\sigma}=2.854833\,(0.135834)$ |  |
$\text{UHN}(\sigma)$ | $\widehat{\sigma}=0.413894\,(0.002855)$ | $-33.62$ | $-32.22$
$\text{ULINDLEY}(\mu)$ | $\widehat{\mu}=0.186834\,(0.000575)$ | $-41.99$ | $-40.58$

In the following subsection, we will illustrate the performance of the AU
model adopting a high-frequency data set originated from the relative humidity
from a city located in the Atacama Desert.

###  Water Monitoring in Air Humidity

The hydrological regime of the main rivers of Atacama is characterized by ice
sources: water flows from the peaks following the melting of snowfall,
glaciers, and permafrost located in the upper parts of the Andes range. In the
context of climate change, it is therefore essential to understand the
hydrological cycle of these regions in order to set up a sustainable
management policy. Understanding the hydrological cycle requires the
implementation of tools for forecasting river flows, relative humidity,
groundwater reservoirs, or any other water-related quantity monitoring, which
inevitably needs an in-depth knowledge of the physical phenomena that govern
the entire hydrological cycle and, more precisely, the complex interaction
between atmosphere, climate, landforms, ice, snow and river flows.

Additionally, a unique phenomenon called Camanchaca happens, which is a fog
passing by the Copiapó city, recurrent only between midnight until around 10
a.m. Here, we demonstrate the variation of the relative humidity of Copiapó
city, proposing a methodology that can be efficient, adjustable to these data.
Using the daily maximum relative humidity, we compared six different unit
distributions: AU, BE, KUM, LOGITNO, SIMPLEX, and UHN, as shown in Figure 5.

Figure 5: Estimated densities superimposed to the histogram (top-chart), and
estimated CDFs superimposed to the ECDF (bottom-chart) (relative air humidity
data).

After comparing the commonly used unit models, we demonstrate the advantage of
adjusting the AU model against others (visually). Table 3 confirms the best
fit of the AU model based on information criteria (AIC and BIC), as well as
shows the estimation of the parameter(s) of each model.

Table 3: Parameter estimates, AIC and BIC values (relative air humidity data). Model | Parameter Estimate (S.E.) | AIC | BIC
---|---|---|---
$\text{AU}(\alpha)$ | $\widehat{\alpha}=0.1092\,(3.1902\text{e-}07)$ | $-\text{14,023.49}$ | $-\text{14,016.76}$
$\text{BE}(\mu,\sigma)$ | $\widehat{\mu}=0.8476\,(1.2027\text{e-}06)$ | $-\text{13,927.89}$ | $-\text{13,914.41}$
| $\widehat{\sigma}=0.2410\,(4.1119\text{e-}06)$ |  |
$\text{KUM}(\mu,\sigma)$ | $\widehat{\mu}=9.4004\,(0.0141)$ | $-\text{13,605.90}$ | $-\text{13,592.43}$
| $\widehat{\sigma}=2.3882\,(0.0019)$ |  |
$\text{LOGITNO}(\mu,\sigma)$ | $\widehat{\mu}=0.8693\,(3.1376\text{e-}06)$ | $-\text{7,600.43}$ | $-\text{7,586.95}$
| $\widehat{\sigma}=1.2299\,(1.2148\text{e-}04)$ |  |
$\text{SIMPLEX}(\mu,\sigma)$ | $\widehat{\mu}=0.9735\,(1.2959\text{e-}06)$ | 32,477.13 | 32,490.61
| $\widehat{\sigma}=94.0480\,(0.7103)$ |  |
$\text{UHN}(\sigma)$ | $\widehat{\sigma}=99.9900\,(6.5334\text{e-}07)$ | 5,101,018,733.13 | 5,101,018,739.86

After obtaining the parameter estimate for $\alpha$, the AU model was used to
construct a Statistical Process Control (SPC) chart, by calculating a
tolerance upper-lower bound. Moreover, it was adopted the Highest Density
Interval (HDI), considering a confidence of 99% to monitor the daily maximum
relative humidity records (as showed by Figure 6).

Figure 6: SPC control chart, considering a 99% of tolerance based on the AU
model fitted to the daily maximum relative humidity of Copiapó city, Chile,
from February 1st, 2015 until October 4th, 2022. It is observed that 193 days
(3.1%) presented anomaly values (out-of-control signals). The obtained control
limits were: $\text{LCL}=68.56\%$ and $\text{UCL}=97.73\%$.

The expected daily maximum water relative humidity is 76.48% (based on the
adjusted AU model). The obtained control limits, considering a confidence (or
tolerance) of 99%, were: $\text{LCL}=68.56\%$ and $\text{UCL}=97.73\%$. Thus,
the control chart based on the AU model, or simply, AU control chart, is
another exciting and valuable alternative to some well-known SPC tools, which
enlightens the prediction and opens new doors to discuss extreme events in the
Atacama water particles monitoring through probabilistic reasoning.

##  Conclusions

This work showed the competitiveness of the developed Theorem 1 (equation
(1.2)), which enables a great class of distributions that are all from the
exponential family. As an exemplification, we adopted the special case for
$k=1$, equivalent to the moment of order two of the standard normal
distribution, and after some transformations, developed the Alpha-Unit (AU)
distribution. Whereas we dedicated to the unit range given the importance of
this stochasticity representation.

Unit distributions are useful for values that oscillate between 0 and 1, such
as fractions, proportions and rates, among others, or for a set of values
where there is a minimum or maximum limitation, resorting to standardization
through the min-max transformation. Most distributions of this type come from
transforming a random variable with certain distribution so that it takes
values between 0 and 1, as in the case of unit-Lindley distribution
[mazucheli2019ul], which comes from the Lindley distribution [lindley1958,
lindley1965].

There are numerous works based on (unit) distributions, extending a model and
applying it to several areas [korkmaz2021unitll, guerra2021unit,
korkmaz2022unit]. In this work, we introduced and showed the competitiveness
of the AU distribution, especially for data with a range greater than 0.4, or
which present high asymmetry and low decay. Further studies shall investigate
this hypothesis in a wider amount of data sets (through different sorts of
wide data range). Additionally, implementation in this model adopting
hierarchical estimation and spatio-temporal dependence would be useful for
forecast/predictable problems.

## References

  * [1] [Afify et al.]Afify, Nassar, Kumar Cordeiro2022afify2022new Afify, A. Z., Nassar, M., Kumar, D. Cordeiro, G. M. 2022, ‘A new unit distribution: properties, inference, and applications’, Electronic Journal of Applied Statistical Analysis 15(2), 460–484.
  * [2] Akaike et al.1977akaike1977 Akaike, H. et al. 1977, ‘On entropy maximization principle.’, Application of statistics .
  * [3] [Bakouch et al.]Bakouch, Nik, Asgharzadeh Salinas2021bakouch2021flexible Bakouch, H. S., Nik, A. S., Asgharzadeh, A. Salinas, H. S. 2021, ‘A flexible probability model for proportion data: Unit-half-normal distribution’, Communications in Statistics-Case Studies, Data Analysis and Applications 7(2), 271–288.
  * [4] Barndorff-Nielsen Jørgensen1991barndorff1991some Barndorff-Nielsen, O. E. Jørgensen, B. 1991, ‘Some parametric models on the simplex’, Journal of Multivariate Analysis 39(1), 106–116.
  * [5] Bayer, Cintra Cribari-Neto2018bayer2018beta Bayer, F. M., Cintra, R. J. Cribari-Neto, F. 2018, ‘Beta seasonal autoregressive moving average models’, Journal of Statistical Computation and Simulation 88(15), 2961–2981.
  * [6] Bayer, Tondolo Müller2018bayer2018CIE Bayer, F. M., Tondolo, C. M. Müller, F. M. 2018, ‘Beta regression control chart for monitoring fractions and proportions’, Computers & Industrial Engineering 119, 416–426.
  * [7] Elal-Olivero2010elal2010alpha Elal-Olivero, D. 2010, ‘Alpha-skew-normal distribution’, Proyecciones (Antofagasta) 29(3), 224–240.
  * [8] Ferrari Cribari-Neto2004ferrari2004beta Ferrari, S. Cribari-Neto, F. 2004, ‘Beta regression for modelling rates and proportions’, Journal of applied statistics 31(7), 799–815.
  * [9] [Fonseca et al.]Fonseca, Ferreira, Nascimento, Fiaccone, Ulloa-Correa, García-Piña Louzada2021fonseca2021water Fonseca, A., Ferreira, P. H., Nascimento, D. C. d., Fiaccone, R., Ulloa-Correa, C., García-Piña, A. Louzada, F. 2021, ‘Water particles monitoring in the atacama desert: Spc approach based on proportional data’, Axioms 10(3), 154.
  * [10] [Guerra et al.]Guerra, Peña-Ramírez Bourguignon2021guerra2021unit Guerra, R. R., Peña-Ramírez, F. A. Bourguignon, M. 2021, ‘The unit extended weibull families of distributions and its applications’, Journal of Applied Statistics 48(16), 3174–3192.
  * [11] Korkmaz, Altun, Chesneau Yousof2022korkmaz2022unitc Korkmaz, M. Ç., Altun, E., Chesneau, C. Yousof, H. M. 2022, ‘On the unit-chen distribution with associated quantile regression and applications’, Mathematica Slovaca 72(3), 765–786.
  * [12] Korkmaz, Chesneau Korkmaz2022korkmaz2022unit Korkmaz, M. Ç., Chesneau, C. Korkmaz, Z. S. 2022, ‘The unit folded normal distribution: A new unit probability distribution with the estimation procedures, quantile regression modeling and educational attainment applications’, Journal of Reliability and Statistical Studies pp. 261–298.
  * [13] Korkmaz Korkmaz2021korkmaz2021unitll Korkmaz, M. Ç. Korkmaz, Z. S. 2021, ‘The unit log–log distribution: a new unit distribution with alternative quantile regression modeling and educational measurements applications’, Journal of Applied Statistics pp. 1–20.
  * [14] [Krishna et al.]Krishna, Maya, Chesneau Irshad2022krishna2022unitt Krishna, A., Maya, R., Chesneau, C. Irshad, M. R. 2022, ‘The unit teissier distribution and its applications’, Mathematical and Computational Applications 27(1), 12.
  * [15] Kumaraswamy1980kumaraswamy1980generalized Kumaraswamy, P. 1980, ‘A generalized probability density function for double-bounded random processes’, Journal of Hydrology 46(1-2), 79–88.
  * [16] Limpert Stahel2011limpert2011problems Limpert, E. Stahel, W. A. 2011, ‘Problems with using the normal distribution–and ways to improve quality and efficiency of data analysis’, PloS one 6(7), e21403.
  * [17] Lindley1958lindley1958 Lindley, D. V. 1958, ‘Fiducial distributions and bayes’ theorem’, Journal of the Royal Statistical Society. Series B (Methodological) 20(1), 102–107.
  * [18] Lindley1965lindley1965 Lindley, D. V. 1965, Introduction to Probability and Statistics from a Bayesian Viewpoint, Part II: Inference, Cambridge University Press.
  * [19] [Mazucheli et al.]Mazucheli, Menezes Chakraborty2019mazucheli2019ul Mazucheli, J., Menezes, A. F. B. Chakraborty, S. 2019, ‘On the one parameter unit-lindley distribution and its associated regression model for proportion data’, Journal of Applied Statistics 46(4), 700–714.
  * [20] [Mazucheli et al.]Mazucheli, Menezes Ghitany2018mazucheli2018unitw Mazucheli, J., Menezes, A. F. B. Ghitany, M. E. 2018, ‘The unit-weibull distribution and associated inference’, Journal of Applied Probability and Statistics 13(2), 1–22.
  * [21] [Mazucheli et al.]Mazucheli, Menezes, Fernandes, De Oliveira Ghitany2020mazucheli2020unitw Mazucheli, J., Menezes, A., Fernandes, L., De Oliveira, R. Ghitany, M. 2020, ‘The unit-weibull distribution as an alternative to the kumaraswamy distribution for the modeling of quantiles conditional on covariates’, Journal of Applied Statistics 47(6), 954–974.
  * [22] R Core Team2022Rcran R Core Team 2022, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria.
https://www.R-project.org/

  * [23] [Sagrillo et al.]Sagrillo, Guerra Bayer2021sagrillo2021modified Sagrillo, M., Guerra, R. R. Bayer, F. M. 2021, ‘Modified kumaraswamy distributions for double bounded hydro-environmental data’, Journal of Hydrology 603, 127021.
  * [24] [Santana-e Silva et al.]Santana-e Silva, Cribari-Neto Vasconcellos2022santana2022beta Santana-e Silva, J. J., Cribari-Neto, F. Vasconcellos, K. L. 2022, ‘Beta distribution misspecification tests with application to covid-19 mortality rates in the united states’, Plos One 17(9), e0274781.
  * [25] Schwarz et al.1978schwarz1978 Schwarz, G. et al. 1978, ‘Estimating the dimension of a model’, The annals of statistics 6(2), 461–464.
  * [26] Stahl2006stahl2006evolution Stahl, S. 2006, ‘The evolution of the normal distribution’, Mathematics magazine 79(2), 96–113.
  * [27] Tadikamalla Johnson1982tadikamalla1982systems Tadikamalla, P. R. Johnson, N. L. 1982, ‘Systems of frequency curves generated by transformations of logistic variables’, Biometrika 69(2), 461–465.
  * [28]

## Appendix A

This appendix shows the proof that for a random variable

$X\sim\text{BN}(k)\quad\rightarrow\quad X^{2}\sim\chi^{2}_{2k+1}.$

Then,

$\displaystyle F_{X^{2}}(x)$ $\displaystyle=\mathbb{P}\left(X^{2}\leq
x\right)=\mathbb{P}\left(-\sqrt{x}\leq
X\leq\sqrt{x}\right)=2\mathbb{P}\left(X\leq\sqrt{x}\right)-1=2F_{X}\left(\sqrt{x}\right)-1.$

It follows that

$\displaystyle f_{X^{2}}(x)$
$\displaystyle=2f_{X}\left(\sqrt{x}\right)\frac{1}{2\sqrt{x}}=\frac{1}{c}\left(\sqrt{x}\right)^{2k}\phi(\sqrt{x})\frac{1}{\sqrt{x}}=\frac{1}{\prod^{k}_{j=1}(2j-1)}\left(\sqrt{x}\right)^{2k-1}\frac{1}{\sqrt{2\pi}}e^{-\frac{x}{2}}.$

Knowing that
$\Gamma\left(\frac{2k+1}{2}\right)=\prod^{k}_{j=1}(2j-1)\frac{\sqrt{\pi}}{2^{k}}$,
then

$\displaystyle f_{X^{2}}(x)$
$\displaystyle=\frac{1}{\prod^{k}_{j=1}(2j-1)}x^{\frac{2k-1}{2}}\frac{1}{\sqrt{2\pi}}e^{-\frac{x}{2}}=\frac{\sqrt{\pi}}{2^{k}\Gamma\left(\frac{2k+1}{2}\right)}\frac{x^{\frac{2k-1}{2}}}{2^{1/2}\sqrt{\pi}}e^{-\frac{x}{2}}$
$\displaystyle=\frac{1}{\Gamma\left(\frac{2k+1}{2}\right)2^{\frac{2k+1}{2}}}x^{\frac{2k-1}{2}}e^{-\frac{x}{2}}.$

Therefore, $X^{2}\sim\chi^{2}_{2k+1}$.

Besides that, complementation can be taken into account by saying that,
considering $W_{2}\sim\chi^{2}_{2k+1}$ and $\mathbb{P}(W_{1}=\pm 1)=1/2$, then
$B=W_{1}\sqrt{W_{2}}\sim\text{BN}(k)$:

Let $b\geq 0$, then

$\displaystyle F_{B}(b)$ $\displaystyle=\mathbb{P}(B\leq
b)=\mathbb{P}\left(W_{1}\sqrt{W_{2}}\leq b\right)$
$\displaystyle=\mathbb{P}\left(W_{1}\sqrt{W_{2}}\leq b\mid
W_{1}=1\right)\mathbb{P}\left(W_{1}=1\right)+\mathbb{P}\left(W_{1}\sqrt{W_{2}}\leq
b\mid W_{1}=-1\right)\mathbb{P}\left(W_{1}=-1\right)$
$\displaystyle\overset{\text{ind.}}{=}\mathbb{P}\left((1)\sqrt{W_{2}}\leq
b\right)\frac{1}{2}+\mathbb{P}\left((-1)\sqrt{W_{2}}\leq b\right)\frac{1}{2}.$
$\displaystyle\text{Since }b\geq 0,\ \ \text{then}\ \
\mathbb{P}\left((-1)\sqrt{W_{2}}\leq b\right)=1:$
$\displaystyle=\mathbb{P}\left(\sqrt{W_{2}}\leq
b\right)\frac{1}{2}+\frac{1}{2}=\mathbb{P}\left(|W_{2}|\leq
b^{2}\right)\frac{1}{2}+\frac{1}{2}=\mathbb{P}\left(-b^{2}\leq W_{2}\leq
b^{2}\right)\frac{1}{2}+\frac{1}{2}$
$\displaystyle=\frac{1}{2}\bigg{[}\mathbb{P}\left(X\leq
b^{2}\right)-\underbrace{\mathbb{P}\left(X\leq-b^{2}\right)}_{0}\bigg{]}+\frac{1}{2}=\frac{1}{2}\mathbb{P}\left(X\leq
b^{2}\right)+\frac{1}{2}=\frac{1}{2}F_{X}\left(b^{2}\right)+\frac{1}{2}$
Therefore, $\displaystyle f_{B}(b)$ $\displaystyle=\frac{\partial
F_{B}(b)}{\partial
b}=\frac{1}{2}f_{X}\left(b^{2}\right)2b=bf_{X}\left(b^{2}\right)=b\frac{1}{\Gamma\left(\frac{2k+1}{2}\right)2^{\frac{2k+1}{2}}}\left(b^{2}\right)^{\frac{2k+1}{2}-1}e^{-\frac{b^{2}}{2}}$
$\displaystyle=b\frac{1}{\Gamma\left(\frac{2k+1}{2}\right)2^{\frac{2k+1}{2}}}b^{2k-1}e^{-\frac{b^{2}}{2}}=\frac{1}{\frac{\sqrt{\pi}\prod^{k}_{j=1}(2j-1)}{2^{k}}2^{\frac{2k+1}{2}}}b^{2k}e^{-\frac{b^{2}}{2}}$
$\displaystyle=\frac{1}{\prod^{k}_{j=1}(2j-1)}\frac{b^{2k}}{\sqrt{2\pi}}e^{-\frac{b^{2}}{2}}=\frac{1}{\underbrace{\prod^{k}_{j=1}(2j-1)}_{c}}b^{2k}\phi(b).$

Analogously, it is proved for $b<0$.

## Appendix B

The proposed theorem will be illustrated considering $k=1$ to show the origin
of the random numbers that generate the AU distribution.

Proposition 6. If $X\sim\text{BN}(1)$, then

  1. 1.

$f_{X}(x)=x^{2}\phi(x)$ is a bimodal density function;

  2. 2.

$X^{2}\sim\chi_{3}^{2}$;

  3. 3.

Let $Y$ and $V$ be independent random variables with $Y\sim\chi_{3}^{2}$ and
$V$ is such that $\mathbb{P}(V=-1)=\mathbb{P}(V=1)=1/2$, then

$R=V\sqrt{Y}\sim\text{BN}(1).$

Demonstration.

  1. 1.

If $f_{X}(x)$ is bimodal, it would have two maxima, for which the first and
second derivative criteria would be applied:

$\frac{df_{X}(x)}{dx}=0\quad\Rightarrow$
$\frac{d\left(x^{2}\phi(x)\right)}{dx}=2x\phi(x)+x^{2}[-x\phi(x)]=2x\phi(x)-x^{3}\phi(x)=x\phi(x)(2-x^{2})=0.$

Then, it can be seen that the solutions of the previous equation would be:
$x_{1}=0$, $x_{2}=\sqrt{2}$, $x_{3}=-\sqrt{2}$. Hence, applying the second
derivative criterion:

$\frac{d^{2}f_{X}(x)}{dx^{2}}<0\quad\Rightarrow$
$\frac{d\left(x\phi(x)(2-x^{2})\right)}{dx}=\phi(x)(2-x^{2})+x[-x\phi(x)](2-x^{2})+x\phi(x)(-2x).$

Reducing algebraically, we get:

$\frac{d^{2}f_{X}(x)}{dx^{2}}=\phi(x)\left(2-5x^{2}+x^{4}\right)<0.$

The only solutions that satisfy the previous inequality are: $x_{2}=\sqrt{2}$,
$x_{3}=-\sqrt{2}$. Therefore, there are two maxima and the BN distribution is
bimodal.

  2. 2.

Let $W=X^{2}$ and $w>0$:

$\displaystyle F_{W}(w)$ $\displaystyle=\mathbb{P}(W\leq
w)=\mathbb{P}(X^{2}\leq w)=\mathbb{P}(-\sqrt{w}\leq X\leq\sqrt{w})$
$\displaystyle=\mathbb{P}(X\leq\sqrt{w})-\mathbb{P}(X\leq-\sqrt{w})$
$\displaystyle=F_{X}\left(\sqrt{w}\right)-\left[1-F_{X}\left(\sqrt{w}\right)\right]=2F_{X}\left(\sqrt{w}\right)-1.$

Then, by deriving the previous expression, we obtain:

$f_{W}(w)=2f_{X}\left(\sqrt{w}\right)\frac{1}{2\sqrt{w}}=\frac{1}{\sqrt{w}}f_{X}\left(\sqrt{w}\right)=\frac{1}{\sqrt{w}}\left(\sqrt{w}\right)^{2}\phi\left(\sqrt{w}\right)=\frac{1}{\sqrt{2\pi}}w^{1/2}e^{-w/2}.$

Observing the expression above, we have that $X^{2}\sim\chi_{3}^{2}$.
$\quad\blacksquare$

Definition 3. (Bimodal Half-Normal distribution). Let $Y\sim\text{BN}(1)$, if
$Q=\alpha|Y|$ with $\alpha>0$, then we say that $Q$ is distributed according
to a Bimodal Half-Normal (BHN) distribution with parameter $\alpha$ and we
denote it by $Q\sim\text{BHN}(\alpha)$.

Proposition 7. If $Q\sim\text{BHN}(\alpha)$, then the PDF of $Q$ is given by

$f_{Q}(q\mid\alpha)=\frac{2}{\alpha}\left(\frac{q}{\alpha}\right)^{2}\phi\left(\frac{q}{\alpha}\right),\quad
q>0.\\\ $

Demonstration.

$F_{Q}(q)=\mathbb{P}(Q\leq q)=\mathbb{P}(\alpha|Y|\leq
q)=\mathbb{P}\left(-\frac{q}{\alpha}\leq
Y\leq\frac{q}{\alpha}\right)=2\mathbb{P}\left(Y\leq\frac{q}{\alpha}\right)-1=2F_{Y}\left(\frac{q}{\alpha}\right)-1.$

Then, by deriving the previous expression, we obtain:

$f_{Q}(q)=2f_{Y}\left(\frac{q}{\alpha}\right)\frac{1}{\alpha}=\frac{2}{\alpha}\left(\frac{q}{\alpha}\right)^{2}\phi\left(\frac{q}{\alpha}\right).\quad\blacksquare$

Proposition 8. If $Q\sim\text{BHN}(\alpha)$, then

$X=e^{-Q}\sim\text{AU}(\alpha).\\\ $

Demonstration. Let $X=e^{-Q}$, $0<x\leq 1$, then

$\displaystyle F_{X}(x)$ $\displaystyle=\mathbb{P}(X\leq
x)=\mathbb{P}\left(e^{-Q}\leq
x\right)=\mathbb{P}\left(-Q\leq\ln(x)\right)=\mathbb{P}\left(Q\geq-\ln(x)\right)$
$\displaystyle=1-\mathbb{P}\left(Q\leq-\ln(x)\right)=1-F_{Q}\left(-\ln(x)\right).$

By deriving the previous expression, we have:

$f_{X}(x)=f_{Q}(-\ln(x))\frac{1}{x}=\frac{2}{\alpha}\left(\frac{-\ln(x)}{\alpha}\right)^{2}\phi\left(\frac{-\ln(x)}{\alpha}\right)\frac{1}{x}=\frac{2}{\alpha
x}\left(\frac{\ln(x)}{\alpha}\right)^{2}\phi\left(\frac{\ln(x)}{\alpha}\right).\quad\blacksquare$