11institutetext: Kiwamu Izumi 22institutetext: Institute of Space and
Astronautical Science, Japan Aerospace Exploration Agency, 3-1-1 Yoshinodai,
Chuo, Sagamihara, Kanagawa 252-5210 Japan, 22email<EMAIL_ADDRESS>

Karan Jani 33institutetext: Department Physics and Astronomy, Vanderbilt
University, 2301 Vanderbilt Place, Nashville, TN, 37235, USA , 33email:
<EMAIL_ADDRESS>
# Detection Landscape in the Deci-Hertz Gravitational-Wave Spectrum

Kiwamu Izumi corresponding author Karan Jani

###### Abstract

Direct observations of gravitational waves at frequencies around deci-Hertz
will play a crucial role in fully exploiting the potential of multi-messenger
astronomy. In this chapter, we discuss the detection landscape for the next
several decades of the deci-Hertz gravitational-wave spectrum. We provide an
overview of the experimental frontiers being considered to probe this
challenging regime and the astrophysics and fundamental goals accessible
towards them. This includes interferometric observatories in space with
heliocentric and geocentric satellites, cubesats, lunar-based experiments and
atom intereferometry. A major focus of this chapter is towards the technology
behind DECi-hertz Interferometer Gravitational wave Observatory (DECIGO) and
its scientific pathfinder mission concept B-DECIGO.

##  1 Keywords

Instruments, laser interferometry, space mission, lunar science, cubesats,
atomic clock, gravitational-wave astronomy.

##  2 1\. Introduction

The current generation of ground-based gravitational-wave detectors such as
Laser Interferometer Gravitational-Wave Observatory (LIGO)
TheLIGOScientific:2014jea , Virgo TheVirgo:2014hva  and KAGRA Akutsu:2018axf
are sensitive at frequencies above 10 Hz. These L-shaped detectors with a few
km arm-length are suited to probe gravitational waves from the coalescing
binaries of neutron stars and stellar black holes
$({\lesssim}10^{2}~{}M_{\odot})$ Aasi:2013wya . The next-generation of ground-
based gravitational-wave detectors like Einstein Telescope Punturo:2010zz  and
Cosmic Explorer Reitze:2019iox  are expected to be sensitive as low as 5 Hz.
The arm-length of these detectors will be tens of km and include several
technological upgrades from the current generation. With these next-generation
detectors, the sensitive volume for stellar compact binaries will reach
unprecedented cosmological scale, and open the prospects for direct detections
of lower-range to mid-range intermediate-mass black hole
$({\lesssim}10^{3}~{}M_{\odot})$ Jani:2019ffg . The upcoming space-based
gravitational-wave detector Laser Interferometer Space Antenna (LISA)
Audley:2017drz  will have a few million km long arm length, thus measuring
low-frequency gravitational waves near the milli-Hertz range. This debut
space-based gravitational-wave detector will open a new window into the
astrophysics of super-massive black hole binaries of
${\sim}10^{6}-10^{7}~{}M_{\odot}$. Furthermore, the ongoing network of pulsar
timing arrays IPTA):2013lea  are probing the nano-Hertz regime. This can
potentially measure gravitational waves from super-massive black hole binaries
as heavy as ${\sim}10^{9}~{}M_{\odot}$ Mingarelli:2019pyd .

A particularly challenging regime to measure gravitational waves for all the
detectors mentioned above is near the deci-Hertz frequencies
$({\sim}0.1~{}\mathrm{Hz})$. For a terrestrial laser interferometric setup, it
is difficult to go below 1 Hz barrier due to earth’s seismic noise. On the
other hand, space missions like LISA cannot go higher than 0.1 Hz due to laser
shot noise. The universe is filled with rich astrophysical sources in and
around deci-Hertz $(10^{-2}~{}\mathrm{Hz}\sim 1~{}\mathrm{Hz})$ (see
Sedda:2019uro ; Mandel:2017pzd  for a recent review). The chief among them are
the elusive intermediate-mass black holes (IMBHs:
$10^{2}-10^{5}~{}M_{\odot}$), which are crucial clues in understanding the
first population of stars in the universe, as well as the seed formation of
the super-massive black holes found at the centers of galaxies (see
Greene_IMBH ; Bellovary:2019nib  for a recent review). There are a few
promising candidates for IMBHs across its broad mass-range from
electromagnetic observations, though the first confirmed evidence is
attributed to the remnant of binary black hole merger GW190521, observed by
the LIGO and Virgo detectors Abbott_2020_190521 .

Figure 1: Gravitational-wave frequencies corresponding to the evolution of
some of the prominent astrophysical binaries in the deci-Hertz spectrum. The
inspiral of binary neutron star (blue line) and stellar binary black holes
(maroon line) will be observed in the deci-Hertz for years. For intermediate-
mass black holes, IMBHs, (gray lines), depending on their mass-ranges, the
deci-Hertz access would measure their inspiral, merger and ringdown.

For a coalescing binary, the instantaneous gravitational-wave frequency
$f_{\mathrm{GW}}$ at a time $\tau$ before the merger scales inversely with the
chirp mass $\mathcal{M}$ Peters:1964zz . In Fig. 1, we show the evolution of
various astrophysical binaries around the deci-Hz range. Binaries with IMBHs
would generally undergo all the three phases of the coalescence - inspiral,
merger and ringdown - within in the deci-Hertz band. A binary with an upper-
range $({\sim}10^{5}~{}M_{\odot})$ and mid-range $({\sim}10^{4}~{}M_{\odot})$
IMBHs would undergo inspiral for a few minutes to hours in the deci-Hertz
band, followed by merger and ringdown. Binaries with lower-range
$({\sim}10^{3}~{}M_{\odot})$ IMBHs would inspiral in deci-Hertz spectrum for a
few days, though their merger and ringdown may be better accessible with next-
generation ground-based detectors. The expected binary population of IMBHs
remain fairly unknown and the constraints on their merger rates differ
dramatically from equal-mass systems to intermediate-mass ratios
$(m_{1}/m_{2}\sim 100)$. The deci-Hertz spectrum opens a rare observational
window into testing a variety of formation scenarios of IMBHs to cosmological
distances.

Within the frequency band of the current generation detectors such as LIGO,
Virgo and KAGRA $(f_{\mathrm{GW}}=[10,~{}1000]~{}\mathrm{Hz})$, we can only
measure about a second of evolution of binaries of stellar black holes
$(\mathcal{M}\sim 10~{}M_{\odot})$. In the deci-Hertz band, these binary black
holes stay for months (maroon line). The binary neutron star sources would
evolve within deci-Hertz band up to a few years (blue line). This extended
observation is important to measure the orbital eccentricity and spin
orientation of the binaries. Both of these measurements are key to distinguish
among binary formation channels Stephan:2016kwj ; Rodriguez:2016vmx . A deci-
Hertz range detector also acts as an early-alert system. All the binaries that
will be seen by next-generation ground-based detectors such as Einstein
Telescope and Cosmic Explorer, would be seen at the lower frequency many
months ahead by a deci-Hertz detector. This opens an interesting possibility
of multiband gravitational-wave astronomy, wherein the same astrophysical
source is studied in different frequency bands Jani:2019ffg ; Sesana:2016ljz ;
Cutler:2019krq ; Isoyama:2018rjb . An equivalent analogy of such in astronomy
is that of multi-wavelength electromagnetic observations across radio, X-ray
and other bands. The precise knowledge of sky-localization and distance from
deci-Hertz observation will permit high-latency electromagnetic followups and
also stronger constraints on the Hubble constant with “dark sirens” GLOC .

The deci-Hertz spectrum allows a unique peak into the progenitors of Type Ia
supernovae (such as double white dwarfs Toonen:2012jj ), whose observations
will help for better calibration of standard candles GLOC ; Lunar_Harms . The
low-frequency gravitational waves by the asymmetric ejection of neutrino
during core-collapse supernovae will also be in the deci-Hertz regime
Vartanyan:2020nmt . In context of fundamental physics, the deci-Hertz sources
open a wide range of tests for general relativity, including any presence of
scalar fields and dipole radiation. The multiband observations between deci-
Hertz with LIGO-like ground detectors and LISA-like space missions will allow
an independent check on the pre- and post-merger properties of black holes,
thus aiding the tests of the no-hair theorem. For probing Beyond Standard
Model physics, deci-Hertz will particularly be effective for probing bosonic
fields around intermediate-mass black holes merger, sub-solar dark matter
candidates and stochastic gravitational-wave background from electroweak
symmetry breaking (see Sedda:2019uro  for a recent review).

In this chapter, we will discuss the various experimental frontiers that are
being considered to access this challenging yet rich deci-Hertz gravitational-
wave spectrum. We will briefly discuss the astrophysics and fundamental
physics goals that can be advanced with each of these frontiers, and comment
on their feasibility by the next decade. A major part of this chapter focuses
on the Japanese space-mission, DECIGO Kawamura:2011zz , and its planned
technological development.

##  3 Experimental Frontiers

For an ideal deci-Hertz gravitational-wave detector, we expect its peak
sensitivity around ${\sim}0.1$ Hz to be somewhere between that of LISA near
the milli-Hz range (characteristic strain, $h=\delta{L}/L\sim 10^{-21}$) or
that of LIGO around 100 Hz ($h\sim 10^{-22}$). While that certainly is not a
strict criteria, it sets a target sensitivity for detecting expected
gravitational-wave sources. Fig. 2 highlights the four broad categories of
experimental frontiers that are being considered to achieve the desired
sensitivity in and around the deci-Hertz regime. This includes a space-based
mission with (i) heliocentric or (ii) geocentric orbits, (iii) an experiment
on the lunar surface, and (iv) atomic interferometry. Under each category, we
list some of the prominent detector concepts suggested in the literature,
along with their notable advantages and challenges.

Figure 2: The four broad categories of experimental frontiers being considered
for measuring gravitational-wave frequencies in the deci-Hertz regime. For
each category, we list the prominent detector designs and their salient
features.

We plot the sensitivity curves for some of the proposed deci-Hertz detectors
in Fig. 3. The detector sensitivity is measured as the dimensionless
characteristic strain $h=\delta{L}/L$. The line styles refers to the four
categories of experimental frontiers: heliocentric (straight lines),
geocentric (dotted-dashed lines), lunar-based (dotted line) and atom
interferometry (dashed line). In Fig. 4, we show the horizon distance, i.e.
the maximum detection radius, of each of these detectors towards coalescing
binaries. Notice that other than the cubesat mission SAGE, all the proposed
deci-Hertz detector can reach up to redshift $z\sim 100$ for a certain mass-
range of binaries.

A majority of the proposed efforts are geared towards a large-scale
heliocentric space mission similar to that of LISA. With three satellites,
each separated by 2.5 million kms arm-length forming the interferometer, the
LISA mission can, in-principle, access gravitational waves with frequencies as
low as ${\sim}0.1$ mHz and as high as ${\sim}0.1$ Hz. At this higher frequency
end, LISA’s sensitivity is expected to be at most $h\sim 10^{-20}$. A post-
LISA mission with a similar 3-satellite configuration but an overall higher
sensitivity could achieve the target for deci-Hertz observations. Few of such
design ideas include the Advanced Millihertz Gravitational-wave Observatory
(AMIGO) Baibhav:2019rsa , which was proposed to ESA’s Voyage 2050 call, and
the Chinese-based mission Taiji Hu:2017mde .

Figure 3: Sensitivity of the proposed gravitational-wave experiments that can
measure frequencies near the deci-Hz spectrum (0.1 Hz). The thick black line
refers to the most sensitive heliocentric space-mission in this region,
DECIGO. The maroon dotted-dashed line refers to its pathfinder space-mission
in geocentric orbit, B-DECIGO. Other heliocentric space missions shown are
ALIA (sky blue curve) and DeciHertz Observatory (blue curve). The lunar-based
detector, GLOC, is shown with dotted grey lines. More unconventional
approaches such as cubesat mission SAGE (orange dotted-dashed curve) and a
detector based on atomic interferometer (dashed green lines) are also
highlighted.

To achieve the most of the deci-Hertz science, a detector should have high-
sensitivity across a broad frequency coverage of $0.01-1$ Hz. For space-
missions with LISA-like concepts, it will require shorter arm-length and
technologies to reduce the acceleration noise. Two of such prominent design
ideas are the Advanced Laser Interferometer Antenna (ALIA) Bender:2013nsa ;
Baker:2019pnp  and DeciHertz Observatories (DO) Sedda:2019uro , proposed to
ESA’s Voyage 2050 call as well. Their sensitivities and detection radii are
shown in Fig. 3 and Fig. 4 (sky-blue and dark-blue curves), respectively. Both
these missions could detect IMBH binaries across its three orders of mass-
range $(10^{2}-10^{5}~{}M_{\odot})$ to an unprecedented cosmological distance
(redshift $z\sim 100$). Furthermore, they would also detect stellar
$({\sim}10~{}M_{\odot})$ and supermassive $({\sim}10^{6}~{}M_{\odot})$ binary
black holes binaries to $z\lesssim 10$, thus complementing with the science
goals from LISA and LIGO.

More recently, there has been modest proposals for relatively low-cast space-
missions in geocentric orbits that can access deci-Hertz region, albeit not
with the same sensitivity requirement. This includes a triangular
interferometeric constellation with cubesats, called SagnAc interferometer for
Gravitational wavE (SAGE) Lacour:2018nws . The sensitivity and detection
radius for SAGE is shown in Fig. 3 and Fig. 4 (orange dashed-dotted curve).
The detection sensitivity of SAGE allows a survey of upper-range IMBH binaries
$({\sim}10^{5}~{}M_{\odot})$ within redshift $z\lesssim 1$. Other proposals of
such geocentric mission in the literature includes geosynchronous Laser
Interferometer Space Antenna (gLISA) Tinto:2011nr  and the Geostationary
Antenna for Disturbance-Free Laser Interferometry (GADFLI) McWilliams:2011ra .

Figure 4: The Detection radius measured in the units of cosmological redshift
$z$ for the proposed deci-Hertz gravitational-wave experiments. The colors
refer to the detector sensitivities shown in Fig. 3. The y-axis is equivalent
to the horizon distance for an optimally placed, equal-mass, non-spinning
binary. We compute the horizon for SNR=8. In the case of space-based detectors
(all detectors except GLOC), we put a detector lifetime of 5 years. The x-axis
refers to the total mass of the binary as measured in its source frame, i.e.
with the redshift correction $1/(1+z)$. We compute redshift assuming the
Planck 2018 cosmological parameters.

With the advent of the return to Moon, it is now finally possible to think of
a gravitational-wave detector on the surface of the Moon that would take
advantage of the natural conditions. A few recent proposals that are being
studied to access deci-Hertz regime from the Moon are the Gravitational-wave
Lunar Observatory for Cosmology (GLOC) GLOC , Lunar Gravitational-Wave Antenna
(LGWA) Lunar_Harms  and Lunar Seismic Gravitational wave Antenna (LSGA) LSGA .
We discuss these proposals in more details in the later sections. A promising
detector technology is also emerging from atom-based interferometry, which
could access gravitational-waves in 0.01-1 Hz of frequency range. Some of the
prominent design proposals are Mid-band Atomic Gravitational Wave
Interferometric Sensor (MAGIS) Graham:2017pmn  and the Atomic Experiment for
Dark Matter and Gravity Exploration in Space (AEDGE) Bertoldi:2019tck . The
sensitivity and detection radius for a concept atom interferometer (adapted
from Sedda:2019uro , with parameters based on Kolkowitz:2016wyg ) is shown in
Fig. 3 and Fig. 4 (green dashed-curve). For more details regarding the science
goals and design proposals for atom-interferometer, see the chapter “Quantum
sensors with matter waves for GW observation” in this edition by Bertoldi,
Bouyer and Canuel.

The most ideal deci-Hertz detector, which has been extensively studied in the
literature, is the proposed Japanese space-mission DECi-hertz Interferometer
Gravitational wave Observatory (DECIGO). The sensitivity and detection radius
of DECIGO is shown in Fig. 3 and Fig. 4 (black curve). The detector posses a
rare advantage of accessing gravitational-waves to cosmological distances from
stellar to IMBHs to supermassive black holes. The mission design and the
associated technologies for DECIGO are discussed in the next section.

##  4 DECIGO

DECi-hertz Interferometer Gravitational wave Observatory (DECIGO) is the
future Japanese space gravitational wave antenna concept observing the
frequency band from $\sim 10^{-1}$ to 10 Hz. DECIGO was originally proposed in
2001 Seto:2001qf  to fill the frequency gap between the mHz frequency band
covered by LISA and the audio band by the terrestrial interferometers. The
conceptual mission design calls for four identical clusters of observatories
deployed in heliocentric orbits, as sketched in figure 5. Each of the clusters
consists of three drag-free spacecraft, 1000 km apart from each other in
almost equilateral triangular formation, forming six laser links Kawamura_2006
. Each spacecraft houses two floating test masses serving as proof masses (see
Figure 6). The target sensitivity of a single cluster is set to $4\times
10^{-24}\,\textrm{Hz}^{-1/2}$ at around 1 Hz Kawamura:2020pcg .

Figure 5: Orbit of DECIGO. Four clusters of observatories are distributed in
the heliocentric orbit with two of them placed at the same position.

A unique design choice in DECIGO is the utilization of the Fabry-Perot optical
cavities consisting of a pair of two test masses apart by 1000 km to improve
the sensitivity without demanding a significant laser power. However, this
means that the optical lengths of the Fabry-Perot cavities must be kept at
resonances at a precision level of $10^{-9}$ m or less. Therefore, achieving a
resonance of a Fabry-Perot cavity in orbit has been identified as one of the
key challenges. The characteristic design parameters are summarized in table
1.

Design Parameters | DECIGO | B-DECIGO
---|---|---
Laser wavelength | 515 nm | 515 nm
Laser power | 10 W | 1 W
Arm Length | 1000 km | 100 km
Finesse of the Fabry-Perot cavities | 10 | 100
Diameter of mirror (test mass) | 1 m | 30 cm
Mirror mass | 100 kg | 30 kg
Number of observatory clusters | 4 | 1
Table 1: Design parameters for DECIGO and B-DECIGO
doi:10.1142/S0218271818450013 .

One of the main scientific objectives of DECIGO is to conduct searches for
primordial gravitational wave backgrounds such as those amplified during the
inflation period Grishchuk_1993 ; PhysRevD.79.103501 . Besides, DECIGO will be
capable of delivering new insight into astrophysics and fundamental physics
through the observations of various compact star binary systems Seto:2001qf ;
10.1143/PTP.123.1069  and other sources PhysRevD.81.104043 ;
PhysRevD.75.061302 ; Kakizaki:2015wua ; PhysRevLett.124.041804 . In addition,
heterogeneous observations of merger events with other gravitational wave
detectors in different frequency bands should improve the performance of the
astrophysical parameter estimation.

Figure 6: Conceptual design of the interferometer for a single DECIGO cluster.

###  4.1 B-DECIGO and technology developments

In order to demonstrate the technologies relevant to DECIGO and simultaneously
conduct astrophysically important observations in the deci-Hz band, the
pathfinder mission concept called B-DECIGO was recently proposed
10.1093/ptep/ptw127 . The sensitivity and detection radius of B-DECIGO are
shown in Fig. 3 and Fig. 4 (maroon dashed-dotted line). One of the major
design drivers for B-DECIGO was the observatory performance being able to
observe the stellar mass binary black holes up to a redshift of $z\sim 30$.
The main scientific objective is to study the origin of the stellar mass black
holes by mapping out their mass spectrum and event rate as a function of
redshift. B-DECIGO is also capable of contributing to multi-messenger
astronomy by issuing the forecasts for merger events with a sky localization
error of $\sim 0.3\,\textrm{deg}^{2}$ and a merger time accuracy of $\sim 1$
sec for a binary black hole with the masses similar to GW150914 at $z=0.1$, a
day before the merger.

Based on the scientific motivations as well as the consideration for making
the developments less challenging, B-DECIGO is designed to have a shorter arm
length of 100 km and reduced mass of 30 kg for the test masses with only a
single observatory cluster deployed. Geocentric orbits with a high altitude of
2000 km has been considered. The design parameters of B-DECIGO are also listed
in table 1. The target sensitivity is set to $2\times
10^{-23}\,\textrm{Hz}^{-1/2}$ at around 1 Hz. The launch date is currently
foreseen to be some time in 2030’s.

The recent activities for the technology developments relevant to both DECIGO
and B-DECIGO encompass a number of key functionalities and aspects. These
include the laser source and its stabilization systems, the relative angular
sensor for spacecraft 10.1117/12.2536073 , the interferometer control schemes
nagano2020demonstration ; izumi2020backlinked , the orbital design and further
improvement in the sensitivity YAMADA2020126626 .

##  5 Lunar-Based Experiments

The Moon offers a generous environment for gravitational-wave observations.
There is no atmosphere, thus providing a high-quality vacuum just above the
lunar surface. There are no ocean tides or significant ground motion to impact
the differential displacement of the test masses. There are no large-scale
human activities interfering with the detector. Therefore, in principle, one
one can construct a large-scale interferometer like that of LIGO on lunar-
soil. The conceptual idea has been long around Bender_Lunar , though the halt
of crewed missions to Moon since 70s dissipated any further interest from the
community. It is only recently that initiatives such as NASA Artemis,
Commercial Crew and ESA’s European Large Logistics Lander project that there
is a strong possibility of returning to the Moon this decade and building a
permanent base. One of the science priorities for NASA Artemis is in utilizing
the unique environment of Moon to study the universe.

Unlike space-based deci-Hertz detectors, a lunar-based detector can be in
observational mode for decades. This is particularly important in following up
rarer astrophysical processes in the deci-Hertz spectrum. Hardware failures or
upgrading detector technologies is also possible on the Moon, thus providing a
better value for long-term investment. Below we discuss two distinct
experimental techniques for lunar-based gravitational-wave detection, which in
particular opens the access to deci-Hertz gravitational-wave frequencies.

The first technique is adapted for Gravitational-Wave Lunar Observatory
(GLOC), a proposed interferometer on the Moon GLOC . The conceptual design for
the GLOC is expected to be similar to that of the next-generation earth-based
detectors Cosmic Explorer (CE) and Einstein Telescope (ET). The arm-length of
the interferometer can be set to few tens of kms and the L-shaped LIGIO-like
interferometer could be replaced by a triangular geometry such as that of ET
(see Fig. 7 \- left panel). The end station could be covered in a temperature
controlled dorm to host the required optics. Ideal sites for GLOC could be
inside a crater or within one of lava caves. The tentative sensitivity curve
and detection radius of GLOC is shown in Fig. 3 and Fig. 4 (gray dotted line).
The lowest frequency we expect for GLOC is around ${\sim}0.3$ Hz. Below that
frequency, the sensitivity of GLOC is limited by the seismic and thermal
noises. However, compared to other deci-Hertz detectors, GLOC would remain
sensitive at higher frequencies, thus enhancing the overall signal-to-noise
ratio.

Figure 7: Conceptual designs of gravitational-wave detector on the lunar
surface. Left: Gravitational-wave Lunar Observatory (GLOC) GLOC . Right: Lunar
Gravitational-Wave Antenna (LGWA) Lunar_Harms . Figures have been reproduced
by permission of the corresponding authors.

For stellar mass binaries, GLOC could detect them practically across the
observable universe. The geocentric orbit of Moon helps tremendously in
constraining the sky-localization to few arc-seconds in GLOC, thus opening new
dark sirens for measuring Hubble constant. GLOC could measure sub-solar dark
matter candidates upto redshift $z\sim 10$. For IMBH binaries, GLOC could
survey lower- to mid-range $({\sim}10^{3}-10^{4}~{}M_{\odot})$ upto redshift
$z\sim 10$, while can find the upper-range $({\sim}10^{5}~{}M_{\odot})$ within
redshift $z\lesssim 1$. The rare access to ${\sim}1$ Hz in GLOC further opens
avenues to test progenitor mechanisms to Type Ia, thus providing new
calibrations for the standard candles.

A second technique has been shown in the Lunar Gravitational-Wave Antenna
(LGWA) Lunar_Harms  and the Lunar Seismic and Gravitational-Wave Antenna
(LSGA) LSGA . An extension of the resonant-bar detector idea by Weber
PhysRev.117.306 , in these concept designs the entire Moon acts as a large
spherical detector receiving the incoming gravitational waves. The measurement
of the signal is then essentially the vibrational eigenmodes. For LGWA, the
proposal is to place an array of highend seismometers in km long array on the
Moon (see Fig. 7-right panel). They will monitor the normal modes of the Moon
in the frequency spectrum of 1 mHz — 1 Hz. The predicted sensitivity of LGWA
near the deci-Hertz regime would be better than that of LISA, thus serving an
important science driver. For LSGA, the proposal is to deploy 10 km optical
cables in an L-shaped geometry. The laser light passing through the cables
would permit the measurement of the displacement strain from incoming
gravitational waves. Both these designs will be able to achieve the science
goals of deci-Hertz gravitational-wave astronomy discussed earlier in the
chapter.

##  6 Summary

Continuing experimental effort have been seen over the world for designing
gravitational-wave detectors sensitive at the frequencies around deci-Hertz.
Observations in the deci-Hertz band will provide us with unique opportunities
to enhance the ability of the detector network and to study the dark and
relativistic aspects of our Universe. Among the large-scale space-missions,
DECIGO, its scientific pathfinder mission concept B-DECIGO, ALIA and DO are
expected to open the new gravitational window in the deci-Hertz band up to
cosmological distances. For modest geocentric proposals, SAGE will be able to
survey the deci-Hertz spectrum in the local universe. The lunar-based
experiments GLOC, LGWA and LSGA offers a new avenue to measure frequencies
near deci-Hertz, and substantially advance the scope of gravitational-wave
astronomy.

##  7 Cross-Reference

Bertoldi A., Bouyer P. and Canuel, B., “Quantum sensors with matter waves for
GW observation”

## References

  * (1) Aasi J et al. (LIGO Scientific) 2015 Class. Quant. Grav. 32 074001 (Preprint 1411.4547)
  * (2) Acernese F et al. (Virgo) 2015 Class. Quant. Grav. 32 024001 (Preprint 1408.3978)
  * (3) Akutsu T et al. (KAGRA) 2019 Nature Astron. 3 35–40 (Preprint 1811.08079)
  * (4) Abbott B P et al. (KAGRA, LIGO Scientific, Virgo) 2018 Living Rev. Rel. 21 3 (Preprint 1304.0670)
  * (5) Punturo M et al. 2010 Class. Quant. Grav. 27 194002
  * (6) Reitze D et al. 2019 Bull. Am. Astron. Soc. 51 35 (Preprint 1907.04833)
  * (7) Jani K, Shoemaker D and Cutler C 2019 Nature Astron. 4 260–265 (Preprint 1908.04985)
  * (8) Amaro-Seoane P et al. (LISA) 2017 (Preprint 1702.00786)
  * (9) Manchester R N 2013 Class. Quant. Grav. 30 224010 (Preprint 1309.7392)
  * (10) Mingarelli C M F 2019 (Preprint 1901.06785)
  * (11) Sedda M A et al. 2020 Class. Quant. Grav. 37 215011 (Preprint 1908.11375)
  * (12) Mandel I, Sesana A and Vecchio A 2018 Class. Quant. Grav. 35 054004 (Preprint 1710.11187)
  * (13) Greene J E, Strader J and Ho L C 2020 Annual Review of Astronomy & Astrophysics 58 257–312 (Preprint 1911.09678)
  * (14) Bellovary J, Brooks A, Colpi M, Eracleous M, Holley-Bockelmann K, Hornschemeier A, Mayer L, Natarajan P, Slutsky J and Tremmel M 2019 (Preprint %****␣main.tex␣Line␣325␣****1903.08144)
  * (15) Abbott R, Abbott T, Abraham S, Acernese F, Ackley K, Adams C, Adhikari R, Adya V, Affeldt C, Agathos M and et al 2020 Physical Review Letters 125 ISSN 1079-7114 URL http://dx.doi.org/10.1103/PhysRevLett.125.101102
  * (16) Peters P C 1964 Phys. Rev. 136 B1224–B1232
  * (17) Stephan A P, Naoz S, Ghez A M, Witzel G, Sitarski B N, Do T and Kocsis B 2016 Mon. Not. Roy. Astron. Soc. 460 3494–3504 (Preprint 1603.02709)
  * (18) Rodriguez C L, Zevin M, Pankow C, Kalogera V and Rasio F A 2016 Astrophys. J. 832 L2 (Preprint 1609.05916)
  * (19) Sesana A 2016 Phys. Rev. Lett. 116 231102 (Preprint 1602.06951)
  * (20) Cutler C et al. 2019 Bull. Am. Astron. Soc. 51 109 (Preprint 1903.04069)
  * (21) Isoyama S, Nakano H and Nakamura T 2018 PTEP 2018 073E01 (Preprint 1802.06977)
  * (22) Jani K and Loeb A 2020 arXiv e-prints arXiv:2007.08550 (Preprint 2007.08550)
  * (23) Toonen S, Nelemans G and Portegies Zwart S 2012 Astron. Astrophys. 546 A70 (Preprint 1208.6446)
  * (24) Harms J et al. 2021 The Astrophysical Journal 910 1 (Preprint 2010.13726)
  * (25) Vartanyan D and Burrows A 2020 (Preprint 2007.07261)
  * (26) Kawamura S et al. 2011 Class. Quant. Grav. 28 094011
  * (27) Baibhav V et al. 2019 (Preprint 1908.11390)
  * (28) Hu W R and Wu Y L 2017 Natl. Sci. Rev. 4 685–686
  * (29) Bender P L, Begelman M C and Gair J R 2013 Class. Quant. Grav. 30 165017
  * (30) Mueller G, Baker J et al. 2019 Bull. Am. Astron. Soc. 51 243 (Preprint 1907.11305)
  * (31) Lacour S et al. 2019 Class. Quant. Grav. 36 195005 (Preprint 1811.04743)
  * (32) Tinto M, de Araujo J C N, Aguiar O D and da Silva Alves M E 2011 (Preprint 1111.2576)
  * (33) McWilliams S T 2011 (Preprint 1111.3708)
  * (34) Katsanevas S et al. 2020 Ideas for exploring the Moon with a large European lander (ESA) URL https://ideas.esa.int/servlet/hype/IMT?documentTableId=45087607031744010&userAction=Browse&templateName=&documentId=a315450fae481074411ef65e4c5b7746
  * (35) Graham P W, Hogan J M, Kasevich M A, Rajendran S and Romani R W (MAGIS) 2017 (Preprint 1711.02225)
  * (36) El-Neaj Y A et al. (AEDGE) 2020 EPJ Quant. Technol. 7 6 (Preprint 1908.00802)
  * (37) Kolkowitz S, Pikovski I, Langellier N, Lukin M D, Walsworth R L and Ye J 2016 Phys. Rev. D94 124043 (Preprint 1606.01859)
  * (38) Seto N, Kawamura S and Nakamura T 2001 Phys. Rev. Lett. 87 221103 (Preprint astro-ph/0108011)
  * (39) Kawamura S, Nakamura T, Ando M, Seto N, Tsubono K et al. 2006 Classical and Quantum Gravity 23 S125–S131 URL https://doi.org/10.1088%2F0264-9381%2F23%2F8%2Fs17
  * (40) Kawamura S et al. 2020 (Preprint 2006.13545)
  * (41) Kawamura S, Nakamura T, Ando M, Seto N, Akutsu T et al. 2019 International Journal of Modern Physics D 28 1845001 (Preprint https://doi.org/10.1142/S0218271818450013) URL https://doi.org/10.1142/S0218271818450013
  * (42) Grishchuk L P 1993 Classical and Quantum Gravity 10 2449–2477 URL https://doi.org/10.1088%2F0264-9381%2F10%2F12%2F006
  * (43) Kuroyanagi S, Chiba T and Sugiyama N 2009 Phys. Rev. D 79(10) 103501 URL https://link.aps.org/doi/10.1103/PhysRevD.79.103501
  * (44) Yagi K and Tanaka T 2010 Progress of Theoretical Physics 123 1069–1078 ISSN 0033-068X (Preprint https://academic.oup.com/ptp/article-pdf/123/6/1069/5373653/123-6-1069.pdf) URL https://doi.org/10.1143/PTP.123.1069
  * (45) Nishizawa A, Taruya A and Kawamura S 2010 Phys. Rev. D 81(10) 104043 URL https://link.aps.org/doi/10.1103/PhysRevD.81.104043
  * (46) Seto N 2007 Phys. Rev. D 75(6) 061302 URL https://link.aps.org/doi/10.1103/PhysRevD.75.061302
  * (47) Kakizaki M, Kanemura S and Matsui T 2015 Phys. Rev. D92 115007 (Preprint 1509.08394)
  * (48) Dror J A, Hiramatsu T, Kohri K, Murayama H and White G 2020 Phys. Rev. Lett. 124(4) 041804 URL https://link.aps.org/doi/10.1103/PhysRevLett.124.041804
  * (49) Nakamura T, Ando M, Kinugawa T, Nakano H, Eda K, Sato S, Musha M, Akutsu T, Tanaka T, Seto N, Kanda N and Itoh Y 2016 Progress of Theoretical and Experimental Physics 2016 ISSN 2050-3911 093E01 (Preprint https://academic.oup.com/ptep/article-pdf/2016/9/093E01/9621962/ptw127.pdf) URL https://doi.org/10.1093/ptep/ptw127
  * (50) Musha M, Akami K and Suemasa A 2019 Novel satellite positioning system using acousto-optic deflector International Conference on Space Optics — ICSO 2018 vol 11180 ed Sodnik Z, Karafolas N and Cugny B International Society for Optics and Photonics (SPIE) pp 1589 – 1596 URL https://doi.org/10.1117/12.2536073
  * (51) Nagano K, Takeda H, Michimura Y, Uchiyama T and Ando M 2020 Demonstration of a dual-pass differential fabry-perot interferometer for future interferometric space gravitational wave antennas (Preprint 2008.12462)
  * (52) Izumi K and Fujimoto M K 2020 A back-linked fabry-perot interferometer for space-borne gravitational wave observations (Preprint 2011.05483)
  * (53) Yamada R, Enomoto Y, Nishizawa A, Nagano K, Kuroyanagi S, Kokeyama K, Komori K, Michimura Y, Naito T, Watanabe I, Morimoto T, Ando M, Furusawa A and Kawamura S 2020 Physics Letters A 384 126626 ISSN 0375-9601 URL http://www.sciencedirect.com/science/article/pii/S037596012030493X
  * (54) Stebbins R T and Bender P L 1990 AIP Conference Proceedings 202 188–204 (Preprint https://aip.scitation.org/doi/pdf/10.1063/1.39103) URL https://aip.scitation.org/doi/abs/10.1063/1.39103
  * (55) Weber J 1960 Phys. Rev. 117(1) 306–313 URL https://link.aps.org/doi/10.1103/PhysRev.117.306


longtable

# Bayesian non-parametric specification of bathtub shaped hazard rate
functions

Richard Arnold1, Stefanka Chukova1 and Yu Hayakawa2

(14 May 2023)

1School of Mathematics and Statistics, Victoria University of Wellington, New
Zealand,
<EMAIL_ADDRESS><EMAIL_ADDRESS>
2School of International Liberal Studies, Waseda University, Tokyo, Japan,
<EMAIL_ADDRESS>
#### Abstract

Hazard rate functions of natural and manufactured systems often show a bathtub
shaped failure rate. A high early rate of failures is followed by an extended
period of useful working life where failures are rare, and finally the failure
rate increases as the system reaches the end of its life. Parametric modelling
of such hazard rate functions can lead to unnecessarily restrictive
assumptions on the function shape, however the most common non-parametric
estimator (the Kaplan-Meier estimator) does not allow specification of the
requirement that it be bathtub shaped. In this paper we extend the Lo and Weng
(1989) approach and specify four non-parametric bathtub hazard rate functions
drawn from Gamma Process Priors. We implement and demonstrate simulation for
these four models.

Keywords: Reliability, Bayesian non-parametrics, Gamma Process, Bathtub hazard
rate function.

Table 1: Notation $\alpha$ |  Shape parameter of the Gamma Process Prior
---|---
$a$ |  Symmetry point of Lo-Weng Bathtub
$a_{1},a_{2}$ |  Parameters of hyperprior for $\alpha$
$\beta$ |  Scale parameter of the Gamma Process Prior
$b_{1},b_{2}$ |  Parameters of hyperprior for $\beta$
$\text{DPP}(\alpha H_{0})$ |  Dirichlet Process Prior
$\phi$ |  Parameter of exponential base probability measure $H_{0}$
$f_{1},f_{2}$ |  Parameters of hyperprior for $\phi$
$\text{Beta}(\cdot,\cdot)$ |  Beta distribution
$\text{Exp}(\cdot)$ |  Exponential distribution
$G(\cdot)$ |  Random measure on $\Theta$
$G_{K}(\cdot)$ |  Realisation of the truncated stick-breaking construction of $\text{GaPP}_{K}(\alpha H_{0},\beta)$
$\text{Ga}(\cdot,\cdot)$ |  Gamma distribution
$\text{GaPP}(\alpha H_{0},\beta)$ |  Gamma Process Prior
$\gamma$ |  Total mass of measure $G(\Theta)$
$H_{0}(\cdot)$ |  Probability measure on $\Theta$
$\text{IHPP}(\lambda(\cdot))$ |  Inhomogeneous Poisson Process with hazard rate function $\lambda(\cdot)$
$\kappa(t|u)$ |  Mixing kernel
$K$ |  Truncation limit for the stick-breaking construction
$\lambda(t|\cdot)$ |  Hazard rate function
$\Lambda(t|\cdot)$ |  Cumulative hazard rate function
$n$ |  Total sample size
$n_{0}$ |  Number of uncensored observations
$t$ |  Failure time
$\tau$ |  Censoring time
$\Theta$ |  Set on which random measure $G(\cdot)$ is defined
$\bm{\theta}$ |  Set of $K$ support points
$\theta_{k}$ |  Location of the $k^{\rm th}$ support point
v |  Set of $K-1$ stick-breaking values
$v_{k}$ |  $k^{\rm th}$ stick-breaking value
$w_{k}$ |  Weight of the $k^{\rm th}$ support point
$\tilde{w}_{k}$ |  Unscaled weight of the $k^{\rm th}$ support point

##  1 Introduction

Hazard rate functions of natural and manufactured systems often show a bathtub
shaped failure rate. A high early rate of failures is followed by an extended
period of useful working life where failures are rare, and finally the failure
rate increases as the system reaches the end of its life.

An example of such data is shown in Figure 1, where a histogram of failure
times of one Load-Haul-Dump machine from a Swedish mine (Kumar, Klefsjö, and
Granholm 1989). The associated empirical survival function is plotted
alongside. In this example the system is repairable, and subject to multiple
failures. These are more frequent earlier and later in the machine’s working
life, associated with steep decline in the survival function in the early and
late periods, and a flatter decline between.

Figure 1: (a) Failure times of Load Hauling Machine A from Kumar, Klefsjö, and
Granholm (1989). (b) Kaplan-Meier survival function of the machine failure
times.

Data of this nature can also arise from observations of a set of independent
items, where the length of the useful working life of each item is recorded.
In this setting there is a population of items, some proportion of which have
manufacturing defects, and fail early. The remainder have low failure
probability until late times, at which point they begin failing in large
numbers.

In either scenario the hazard rate starts high, decreases to some minimum, and
then climbs again. A ‘U-shaped’ hazard rate is one where the low hazard rate
period is non-existent or very short, and a ‘bathtub shaped’ hazard rate is
one where there is an extended period where the failure rate is low.

In both cases, the modelling of the associated hazard rate function may be
carried out in a number of ways. Inserting a polynomial dependence on time
into simpler baseline failure distributions is one approach, and another is to
use a two component mixture of parametric distributions (Glaser 1980). The
three phases of the bathtub function can also be built from a piecewise
combination of parametric functions (e.g. Peng, Liu, and Wang 2016). A further
approach is to exponentiate a baseline failure distribution (e.g. Abbasi et
al. 2019; Iqbal et al. 2021). Although such approaches allow a wide range of
behaviours, parametric modelling can lead to unnecessarily restrictive
assumptions on the hazard rate function shape. Conversely, the most common
non-parametric estimator of the hazard rate function (the Kaplan-Meier
estimator, Kaplan and Meier (1958)) is not sufficiently constrained, in that
it does not allow specification of the requirement that it be bathtub shaped.
In this paper we specify and investigate non-parametric bathtub hazard rate
functions based on the Gamma Process Priors, as described by Lo and Weng
(1989). Although the Gamma Process Prior does embed a parametric baseline
failure distribution, it allows the posterior distribution of the failure
distribution to depart significantly from that baseline.

The Gamma Process has been applied to reliability in many settings (e.g
Dykstra and Laud 1981), and in particular it has been used as a conjugate
prior for hazard rate functions. Müller et al. (2015) and Phadia (2016)
provide relevant introductions to Bayesian non-parametrics. The theory of
completely random measures by Kingman (1967) led to the development of the
Dirichlet Process Prior (Ferguson 1974), and in particular its implementations
in the Stick-Breaking and the Chinese Restaurant Processes (Müller et al.
2015; Sethuraman 1994; Paisley 2010).
These advances, and the Stick-Breaking representations of the Beta (Paisley et
al. 2010; Paisley, Blei, and Jordan 2012) and Gamma Processes (Roychowdhury
and Kulis 2014, 2015) have made Bayesian non-parametric approaches feasible
for inference problems.

The work by Dykstra and Laud (1981) in Bayesian non-parametric approaches to
reliability defined a prior over hazard rates using the extended gamma process
and were able to derive posterior distributions for monotonic hazard rate
function in analytic form. Lo and Weng (1989) used a variety of kernel
functions to weight a draw from a Gamma Process and thus generate hazard rate
functions of various forms. The weighted gamma process was also used by Ho and
Lo (2001) as a prior for hazard rate functions, and they implemented a Markov
Chain Monte Carlo method for deriving posterior properties of the hazard rate.

Hayakawa et al. (2001) used the Lo and Weng (1989) formulation to set up a
Bayesian hypothesis test for non-decreasing hazard rates. They also used a
Monte Carlo method, based on the weighted Chinese restaurant process (Lo,
Brunner, and Chan 1998).

Arnold, Chukova, and Hayakawa (2020) demonstrated the Lo and Weng (1989)
approach for hazard rate functions that were increasing (IFR) using the Gamma
Process Prior, and implemented a fully Bayesian non-parametric approach to
inference for such functions. In this paper we extend the results of that
paper to four specifications of bathtub failure rate hazard rate functions,
all based on Gamma Process Priors. We use a gamma-scaled Dirichlet Process
prior to implement the Gamma Process prior, and demonstrate the methodology
and the properties of the models using simulations.

Section 2 introduces the basic concepts of reliability modelling, the Gamma
Process Prior, and then lays out the definitions of six hazard rate functions.
Section 3 gives specific details of models conditional on draws from the Gamma
Process Prior, and protocols for simulation. We also specify the priors
required for each model. Section 4 demonstrates the properties of the models,
showing specific draws from the Gamma Process Prior and the properties of the
hazard rate functions and the failure time distributions that result. Brief
concluding remarks close the paper in Section 5, including a proposal for
future work, including inference, which we will carry out in a subsequent
paper. We list our notation in Table 1.

##  2 Model

###  2.1 Basic Model Specification

We assume that a common hazard rate function $\lambda(t)$ generates failures
in $n$ identical items. The failure time distribution has density
$f(t|\lambda(\cdot))$ and cumulative distribution function
$F(t|\lambda(\cdot))$. The cumulative hazard function is

$\Lambda(t)=\int_{0}^{t}\lambda(u)\,{\rm d}u\ .$ (1)

with survival function

$\bar{F}(t|\lambda(\cdot))=1-F(t|\lambda(\cdot))=e^{-\Lambda(t)}\ .$ (2)

If observations are censored at time $\tau$, then the likelihood of $n$
observations, $n_{0}$ observed and $n-n_{0}$ censored, is

$L(\lambda(\cdot)|\\{t_{i}\\}_{i=1}^{n})=\left[\prod_{i=1}^{n_{0}}\lambda(t_{i})\right]e^{-\sum_{i=1}^{n_{0}}\Lambda(t_{i})-(n-n_{0})\Lambda(\tau)}$
(3)

Our interest is in the estimation of $\lambda(t)$.

In the construction of Lo and Weng (1989) the hazard rate function is formed
by the following weighted integral

$\lambda(t\,|\,G)=\int_{\Theta}\kappa(t\,|\,u)G(du)$ (4)

where $G(\cdot)$ is a member of the space of locally finite measures on the
space $\Theta$ and is the hyperparameter specifying the Gamma Process Prior
for $G$ (see below for its specification). Three particular forms for the
kernel $\kappa(t\,|\,u)$ were proposed by Lo and Weng (1989) as follows:

$\displaystyle\kappa(t\,|\,u)$ $\displaystyle=$ $\displaystyle I(0<u<t)$
$\displaystyle\kappa(t\,|\,u)$ $\displaystyle=$ $\displaystyle I(0<t<u)$ (5)
$\displaystyle\kappa(t\,|\,u)$ $\displaystyle=$ $\displaystyle
I(0<u<|t-a|),\,\,\,a>0$

These define an Increasing Failure Rate (IFR), Decreasing Failure Rate (DFR)
and Bathtub failure rate respectively. The latter case we refer to as the ‘Lo-
Weng Bathtub’ (LWB) (see below). We generalise this kernel approach to define
three further bathtub hazard rate functions.

###  2.2 Specific Models

We now specify six models using draws from one or more Gamma Process Priors,
and supplemented by additional parameters as necessary.

Model 1. Increasing Failure Rate (IFR) (After Lo and Weng 1989):

$\lambda(t|\lambda_{0},G)=\lambda_{0}+\int_{0}^{t}G({\rm d}u)$ (6)

Since $G(\cdot)$ is a non-negative measure, the integral in (6) is a non-
decreasing function of $t$, giving an increasing hazard rate. A constant
offset background failure rate is given by the parameter $\lambda_{0}\geq 0$.

Model 2. Decreasing Failure Rate (DFR) (After Lo and Weng 1989):

$\lambda(t|\lambda_{0},G)=\lambda_{0}+\int_{t}^{\infty}G({\rm d}u)$ (7)

The integral in (7) is a non-increasing function of $t$, giving a decreasing
hazard rate, declining to the background rate $\lambda_{0}\geq 0$.

Model 3. Lo-Weng Bathtub (LWB) (After Lo and Weng 1989):

$\lambda(t|\lambda_{0},a,G)=\lambda_{0}+\int_{0}^{|t-a|}G({\rm d}u)$ (8)

This model combines the IFR and DFR behaviours of Models 1 and 2. Here the
hazard rate is DFR between $t=0$ and $t=a\geq 0$, takes its minimum value of
$\lambda_{0}\geq 0$ at $t=a$ and is IFR thereafter. The hazard is symmetric
around the minimum $t=a$.

Model 4. Superposition Bathtub (SBT): To avoid the artificial symmetry of the
Lo-Weng Bathtub Arnold, Chukova, and Hayakawa (2020) proposed the
superposition of two independent DFR and IFR functions.

$\lambda(t|\lambda_{0},G_{1},G_{2})=\lambda_{0}+\int_{t}^{\infty}G_{1}({\rm
d}u)+\int_{0}^{t}G_{2}({\rm d}u)$ (9)

Here $G_{1}$ and $G_{2}$ are random measures drawn independently from Gamma
Process Priors. The support for $G_{1}$ places mass at early times, and gives
DFR behaviour at those times. Late time IFR behaviour is provided by $G_{2}$
which places mass only at later times.

Model 5. Mixture Bathtub (MBT): Here we adopt the standard specification of a
finite mixture model, and specify the survival function:

$\bar{F}(t|\lambda(\cdot|\pi,G_{1},G_{2}))=\pi\bar{F}_{\rm
DFR}(t|G_{1})+(1-\pi)\bar{F}_{\rm IFR}(t|G_{2})$ (10)

Where $0<\pi<1$ and $\bar{F}_{\rm DFR}(t|G)$ and $\bar{F}_{\rm IFR}(t|G)$ are
the respective survival functions from the IFR (Model 1) and DFR (Model 2)
cases defined above. Although this specification does not technically lead to
a bathtub hazard rate function (see the discussion in Section 4 below), it
does allow for the high rates of early and late failures.

Model 6. Log-convex (LCV): If a bathtub hazard rate function is continuous
then its derivative is negative at $t=0$ and is non-decreasing. This suggests
that its derivative could be modelled using the Gamma Process to give non-
negative increments. In order to ensure that the hazard rate function is
itself non-negative we apply these non-negative incremements to the derivative
of the log hazard rate:

$\frac{d\log\lambda(t|G)}{dt}=w_{0}+\int_{0}^{t}G({\rm d}u)$ (11)

which implies that

$\lambda(t|G)=\lambda_{0}e^{w_{0}t+\int_{0}^{t}(t-u)G({\rm d}u)}$ (12)

Here the constant $w_{0}\geq 0$ gives an IFR distribution, but $w_{0}<0$ gives
a bathtub.

###  2.3 Gamma Process Prior

In the hazard rate specifications above we assume $G$ (or $G_{1}$, $G_{2}$ in
the SBT and MBT cases) is a draw from a Gamma Process prior
$G\sim\text{GaPP}(\alpha H_{0},\beta)$ with shape parameter $\alpha>0$, rate
parameter $\beta>0$ and base probability measure $H_{0}(\cdot)$ defined on the
space $\Theta=\mathbb{R}^{+}$. A draw from the Gamma Process Prior can be
formed by drawing first from a Dirichlet Process Prior $\text{DPP}(\alpha
H_{0})$ and then scaling the resulting weights by an independent Gamma random
variable drawn from $\gamma\sim\text{Gamma}(\alpha,\beta)$. To form the draw
from the Dirichlet Process Prior we use the stick-breaking construction of
Sethuraman (1994). We note that an alternative stick-breaking construction for
the Gamma Process prior was developed by Roychowdhury and Kulis (2014) and
Roychowdhury and Kulis (2015).

In the gamma-scaled Dirichlet Process Prior method, a draw from the prior
$\text{GaPP}(\alpha H_{0},\beta)$ is

$G(du)=\gamma\sum_{k=1}^{\infty}v_{k}\prod_{\ell<k}(1-v_{\ell})\delta_{\theta_{k}}(du)=\sum_{k=1}^{\infty}w_{k}(\gamma,\textbf{v})\delta_{\theta_{k}}(du)$
(13)

where $\gamma\sim\mbox{Ga}(\alpha,\beta)$,
$v_{k}\overset{\text{iid}}{\sim}\mbox{Beta}(1,\alpha)$, and
$\theta_{k}\overset{\text{iid}}{\sim}H_{0}$. A draw from this prior is
discrete: there is a countably infinite set of locations
$\\{\theta_{k}\\}_{k=1}^{\infty}$ and each is associated with a weight
$\\{w_{k}\\}_{k=1}^{\infty}$. These weights are the product of the Gamma
random draw $\gamma$ and the usual Dirichlet Process Prior (DPP) weights
$\\{\tilde{w}_{k}\\}_{k=1}^{\infty}$. These unscaled DPP weights add to 1, the
scaled weights $w_{k}$ sum up to the total mass $\gamma$, and so we have

$w_{k}(\gamma,\textbf{v})=\gamma\tilde{w}_{k}(\textbf{v})=\gamma
v_{k}\prod_{\ell<k}(1-v_{\ell})$ (14)

for $k=1,2,\ldots$. The weights $w_{k}$ are stochastically decreasing in $k$
which means that if we truncate the sum in (13) at some sufficiently large
finite $K$ then the finite sum of the weights $\sum_{k=1}^{K}w_{k}$ will be
very close to $\gamma$. In practice we draw $K$ locations $\theta_{k}$, and
construct the first $K-1$ weights using (14). The weight of the last location
is assigned to be $w_{K}=\gamma-\sum_{k=1}^{K-1}w_{k}$ so that the weights
exactly add to $\gamma$. If $K$ is sufficiently large then the effect of this
approximation is negligible.

Conditional on this truncation a draw $G$ from the prior
$G\sim\text{GaPP}_{K}(\alpha H_{0},\beta)$ can be generated as follows:

$\displaystyle\theta_{k}|H_{0}$ $\displaystyle\overset{\rm iid}{\sim}H_{0}$
for $k=1,2,\ldots,K$ (15) $\displaystyle v_{k}|\alpha$
$\displaystyle\overset{\rm iid}{\sim}\text{Be}(1,\alpha)$ for
$k=1,2,\ldots,K-1$ $\displaystyle\tilde{w}_{k}$
$\displaystyle=v_{k}\prod_{\ell<k}(1-v_{\ell})$ for $k=1,2,\ldots,K-1$
$\displaystyle\tilde{w}_{K}$
$\displaystyle=1-\sum_{k=1}^{K-1}\tilde{w}_{k}=\prod_{\ell=1}^{K-1}(1-v_{\ell})$
$\displaystyle\gamma$ $\displaystyle\sim\text{Ga}(\alpha,\beta)$
$\displaystyle G({\rm d}u)$
$\displaystyle=\gamma\sum_{k=1}^{K}\tilde{w}_{k}\delta_{\theta_{k}}({\rm d}u)$

Integrals over $G(\cdot)$ needed in the evaluation of the hazard rate and
cumulative hazard rate are then

$\displaystyle\int_{0}^{t}G({\rm d}u)$
$\displaystyle=\gamma\sum_{k=1}^{K}\tilde{w}_{k}I(\theta_{k}<t)$ (16)
$\displaystyle\int_{t}^{\infty}G({\rm d}u)$
$\displaystyle=\gamma\sum_{k=1}^{K}\tilde{w}_{k}I(\theta_{k}>t)$

and

$\displaystyle\int_{0}^{t}\int_{0}^{u}G({\rm d}v)$
$\displaystyle=\gamma\sum_{k=1}^{K}\tilde{w}_{k}\text{max}(t-\theta_{k},0)$
(17) $\displaystyle\int_{0}^{t}\int_{u}^{\infty}G({\rm d}v)$
$\displaystyle=\gamma\sum_{k=1}^{K}\tilde{w}_{k}\text{min}(t,\theta_{k})$

A fully hierarchical Bayesian model for $G$ can be completed by choosing a
form for $H_{0}$, with parameters $\phi$, and specifying suitable priors for
$\alpha$, $\beta$, $\gamma$ and $\phi$. For example if
$H_{0}=\text{Exp}(\phi)$ then we can set

$\displaystyle\alpha$ $\displaystyle\sim$
$\displaystyle\text{Ga}(a_{1},a_{2})$ (18) $\displaystyle\beta$
$\displaystyle\sim$ $\displaystyle\text{Ga}(b_{1},b_{2})$ $\displaystyle\phi$
$\displaystyle\sim$ $\displaystyle\text{Ga}(f_{1},f_{2})$

for non-negative constants $\\{a_{1},a_{2},b_{1},b_{2},f_{1},f_{2}\\}$.

##  3 Hazard Rate Functions and Simulation

If we combine the definitions of the hazard rate functions from §2.2 with the
implementation of the Gamma Process Priof in §2.3, we can obtain expressions
for the hazard and cumulative hazard rate functions. We do this for each model
in turn, including suggestions for priors to fully specify the model, and at
the same time provide a protocol for simulation of failure times.
We note that if $U$ is a draw from a Uniform$(0,1)$ distribution, then the
solution $T$ of the equation

$\Lambda(T)=-\log U$

is a draw from the failure distribution with cumulative hazard rate function
$\Lambda(t)$.

In what follows we write Write $w_{k}=\gamma\tilde{w}_{k}$, and also define
$k^{\ast}(t)=\text{max}\\{k\in(1,\ldots,K)\,:\,\theta^{\ast}_{k}\leq t\\}$
which implies that

$\theta^{\ast}_{k^{\ast}(t)}\leq t<\theta^{\ast}_{k^{\ast}(t)+1}$ (19)

We occasionally need to re-index $\\{(w_{k},\theta_{k})\\}_{k=1}^{K}$ as
$\\{(w^{\ast}_{k},\theta_{k}^{\ast})\\}_{k=1}^{K}$ such that

$\theta^{\ast}_{0}\equiv
0<\theta^{\ast}_{1}<\theta^{\ast}_{2}<\ldots<\theta^{\ast}_{K}$

In doing so we note that if $H_{0}$ is chosen to be absolutely continuous then
the $\theta_{k}$ values are almost surely distinct.

We also define the following partial sums of the ordered weights and
locations:

$\displaystyle C^{\ast}_{\ell}$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{\ell}w^{\ast}_{k}=\gamma-\sum_{k=\ell+1}^{K}w_{k}^{\ast}$
(20) $\displaystyle D^{\ast}_{\ell}$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{\ell}w^{\ast}_{k}\theta^{\ast}_{k}$ (21)

#### Model 1. Increasing Failure Rate (IFR)

The hazard and cumulative hazard rate functions are:

$\displaystyle\lambda(t|\lambda_{0},\gamma,\mathbf{\theta},\mathbf{v})$
$\displaystyle=$
$\displaystyle\lambda_{0}+\gamma\sum_{k=1}^{K}\tilde{w}_{k}I(\theta_{k}\leq
t)$ (22)
$\displaystyle\Lambda(t|\lambda_{0},\gamma,\mathbf{\theta},\mathbf{v})$
$\displaystyle=$
$\displaystyle\lambda_{0}t+\gamma\sum_{k=1}^{K}\tilde{w}_{k}\text{max}(0,t-\theta_{k})$
(23)

If we define $\Lambda^{\ast}_{k}=\Lambda(\theta^{\ast}_{k})$ with
$\theta^{\ast}_{0}=0$ and $\Lambda^{\ast}_{0}=0$, then $\Lambda(t)$ is
piecewise linear between the points
$\\{(\theta^{\ast}_{k},\Lambda^{\ast}_{k})\\}_{k=0}^{K}$. If
$U\sim\text{Uniform}(0,1)$ then $T$ is a draw from this distribution if

$\displaystyle T=\frac{-\log
U+D^{\ast}_{k^{\ast}}}{\lambda_{0}+C^{\ast}_{k^{\ast}}}$ (24)

where
$k^{\ast}=\text{max}\\{k\in\\{1,\ldots,K\\}\,:\,\Lambda^{\ast}_{k}\leq-\log
U\\}$. A suitable prior for $\lambda_{0}|G$ is
$\lambda_{0}|\gamma\sim\text{Exp}(\nu/\gamma)$ for some constant $\nu>0$.

#### Model 2. Decreasing Failure Rate (DFR)

The hazard and cumulative hazard rate functions are:

$\displaystyle\lambda(t|\lambda_{0},\gamma,\mathbf{\theta},\mathbf{v})$
$\displaystyle=$
$\displaystyle\lambda_{0}+\gamma\sum_{k=1}^{K}\tilde{w}_{k}I(\theta_{k}>t)$
(25) $\displaystyle\Lambda(t|\lambda_{0},\gamma,\mathbf{\theta},\mathbf{v})$
$\displaystyle=$
$\displaystyle\lambda_{0}t+\gamma\sum_{k=1}^{K}\tilde{w}_{k}\text{min}(t,\theta_{k})$
(26)

Following the same protocol as in the IFR case: if we define
$\Lambda^{\ast}_{k}=\Lambda(\theta^{\ast}_{k})$ with $\theta^{\ast}_{0}=0$ and
$\Lambda^{\ast}_{0}=0$, then $\Lambda(t)$ is piecewise linear between the
points $\\{(\theta^{\ast}_{k},\Lambda^{\ast}_{k})\\}_{k=0}^{K}$. If
$U\sim\text{Uniform}(0,1)$ then $T$ is a draw from this distribution if

$\displaystyle T=\frac{-\log
U-D^{\ast}_{k^{\ast}}}{\lambda_{0}+\gamma-C^{\ast}_{k^{\ast}}}$ (27)

where
$k^{\ast}=\text{max}\\{k\in\\{1,\ldots,K\\}\,:\,\Lambda^{\ast}_{k}\leq-\log
U\\}$. The same prior for $\lambda_{0}|G$ as in the IFR case can be used:
$\lambda_{0}|\gamma\sim\text{Exp}(\nu/\gamma)$ for some constant $\nu>0$.

#### Model 3. Lo-Weng Bathtub (LWB)

The hazard and cumulative hazard rate functions are:

$\displaystyle\lambda(t|\lambda_{0},a,\gamma,\mathbf{\theta},\mathbf{v})$
$\displaystyle=$
$\displaystyle\lambda_{0}+\gamma\sum_{k=1}^{K}\tilde{w}_{k}I(0<\theta_{k}<|t-a|)$
(28) $\displaystyle=$
$\displaystyle\lambda_{0}+\begin{cases}\gamma\sum_{k=1}^{K}\tilde{w}_{k}I(t<a-\theta_{k})\;\;&\text{if
$t<a$}\\\ \gamma\sum_{k=1}^{K}\tilde{w}_{k}I(t\geq a+\theta_{k})\;\;&\text{if
$t\geq a$}\end{cases}$
$\displaystyle\Lambda(t|\lambda_{0},a,\gamma,\mathbf{\theta},\mathbf{v})$
$\displaystyle=$
$\displaystyle\lambda_{0}t+\begin{cases}\gamma\sum_{k=1}^{K}\tilde{w}_{k}I(\theta_{k}<a)\text{min}(t,a-\theta_{k})\;\;&\text{if
$t<a$}\\\
\gamma\sum_{k=1}^{K}\tilde{w}_{k}\left[\text{max}(0,a-\theta_{k})+\text{max}(t-a-\theta_{k},0)\right]\;\;&\text{if
$t\geq a$}\end{cases}\qquad$ (29)

We combine the $\\{(\theta_{k},w_{k})\\}_{k=1}^{K}$ values with one further
location, weight pair: $(a,0)$.
We order the locations, and their associated weights, forming the set
$\\{(w^{\ast\ast}_{k},\theta^{\ast\ast}_{k})\\}_{k=1}^{K+1}$, and then compute
the cumulative hazard function values
$\Lambda^{\ast\ast}_{k}=\Lambda(\theta^{\ast\ast}_{k})$. If
$U\sim\text{Uniform}(0,1)$ then $T$ is a draw from this distribution if

$\displaystyle
T=\theta^{\ast\ast}_{k^{\ast}}+\frac{\theta^{\ast\ast}_{k^{\ast}+1}-\theta^{\ast\ast}_{k^{\ast}}}{\Lambda^{\ast\ast}_{k^{\ast}+1}-\Lambda^{\ast\ast}_{k^{\ast}}}\left(-\log
U-\Lambda^{\ast\ast}_{k^{\ast}}\right)$ (30)

where
$k^{\ast}=\text{max}\\{k\in\\{1,\ldots,K_{1}+K_{2}\\}\,:\,\Lambda^{\ast\ast}_{k}\leq-\log
U\\}$. The prior $\lambda_{0}|\gamma\sim\text{Exp}(\nu/\gamma)$ for some
constant $\nu>0$ is again suitable.

#### Model 4. Superposition Bathtub (SBT)

We draw $G_{1}$ and $G_{2}$ from separate Gamma Process priors:
$G_{1}\sim\text{GaPP}_{K}(\alpha_{1}H_{01},\beta_{1})$ and
$G_{2}\sim\text{GaPP}_{K}(\alpha_{2}H_{02},\beta_{2})$, and add a prior for
$\lambda_{0}|G_{2}\sim\text{Exp}(\nu/\gamma_{2})$. The hazard and cumulative
hazard rate functions are then:

$\displaystyle\lambda(t|\lambda_{0},\gamma_{1},\mathbf{\theta}_{1},\mathbf{v}_{1},\gamma_{2},\mathbf{\theta}_{2},\mathbf{v}_{2})$
$\displaystyle=$
$\displaystyle\lambda_{0}+\gamma_{1}\sum_{k=1}^{K_{1}}\tilde{w}_{1k}I(\theta_{1k}>t)+\gamma_{2}\sum_{k=1}^{K_{2}}\tilde{w}_{2k}I(\theta_{2k}\leq
t)$ (31)
$\displaystyle\Lambda(t|\lambda_{0},\gamma_{1},\mathbf{\theta}_{1},\mathbf{v}_{1},\gamma_{2},\mathbf{\theta}_{2},\mathbf{v}_{2})$
$\displaystyle=$
$\displaystyle\lambda_{0}t+\gamma_{1}\sum_{k=1}^{K_{1}}\tilde{w}_{1k}\text{min}(t,\theta_{1k})+\gamma_{2}\sum_{k=1}^{K_{2}}\tilde{w}_{2k}\text{max}(0,t-\theta_{2k})\qquad$
(32)

Extend the IFR/DFR approach: pool the $\theta_{1k}$ and $\theta_{2k}$
locations and order them, and their associated weights, forming the set
$\\{(w^{\ast\ast}_{k},\theta^{\ast\ast}_{k})\\}_{k=1}^{K_{1}+K_{2}}$. Compute
the cumulative hazard function values
$\Lambda^{\ast\ast}_{k}=\Lambda(\theta^{\ast\ast}_{k})$. Then if
$U\sim\text{Uniform}(0,1)$ then $T$ is a draw from this distribution if

$\displaystyle T$ $\displaystyle=$
$\displaystyle\theta^{\ast\ast}_{k^{\ast}}+\frac{\theta^{\ast\ast}_{k^{\ast}+1}-\theta^{\ast\ast}_{k^{\ast}}}{\Lambda^{\ast\ast}_{k^{\ast}+1}-\Lambda^{\ast\ast}_{k^{\ast}}}\left(-\log
U-\Lambda^{\ast\ast}_{k^{\ast}}\right)$ (33)

where
$k^{\ast}=\text{max}\\{k\in\\{1,\ldots,K_{1}+K_{2}\\}\,:\,\Lambda^{\ast\ast}_{k}\leq-\log
U\\}$.

#### Model 5. Mixture Bathtub (MBT)

We draw $G_{1}$ and $G_{2}$ as in the Superposition Bathtub case, and then the
hazard and cumulative hazard rate functions are:

$\displaystyle\mbox{}\lambda(t|\pi,\lambda_{01},\gamma_{1},\mathbf{\theta}_{1},\mathbf{v}_{1},\lambda_{02},\gamma_{2},\mathbf{\theta}_{2},\mathbf{v}_{2})$
$\displaystyle=$
$\displaystyle\frac{f(t|\pi,\lambda_{01},\gamma_{1},\mathbf{\theta}_{1},\mathbf{v}_{1},\lambda_{02},\gamma_{2},\mathbf{\theta}_{2},\mathbf{v}_{2})}{\bar{F}(t|\pi,\lambda_{01},\gamma_{1},\mathbf{\theta}_{1},\mathbf{v}_{1},\lambda_{02},\gamma_{2},\mathbf{\theta}_{2},\mathbf{v}_{2})}$
$\displaystyle=$ $\displaystyle\frac{\pi f_{\rm
DFR}(t|\lambda_{01},\gamma_{1},\mathbf{\theta}_{1},\mathbf{v}_{1})+(1-\pi)f_{\rm
IFR}(t|\lambda_{02},\gamma_{2},\mathbf{\theta}_{2},\mathbf{v}_{2})}{\pi\bar{F}_{\rm
DFR}(t|\lambda_{01},\gamma_{1},\mathbf{\theta}_{1},\mathbf{v}_{1})+(1-\pi)\bar{F}_{\rm
IFR}(t|\lambda_{02},\gamma_{2},\mathbf{\theta}_{2},\mathbf{v}_{2})}$
$\displaystyle\mbox{}\Lambda(t|\pi,\lambda_{01},\gamma_{1},\mathbf{\theta}_{1},\mathbf{v}_{1},\lambda_{02},\gamma_{2},\mathbf{\theta}_{2},\mathbf{v}_{2})$
$\displaystyle=$
$\displaystyle-\log(\bar{F}(t|\pi,\lambda_{01},\gamma_{1},\mathbf{\theta}_{1},\mathbf{v}_{1},\lambda_{02},\gamma_{2},\mathbf{\theta}_{2},\mathbf{v}_{2}))$
$\displaystyle=$ $\displaystyle-\log\left(\pi\bar{F}_{\rm
DFR}(t|\lambda_{01},\gamma_{1},\mathbf{\theta}_{1},\mathbf{v}_{1})+(1-\pi)\bar{F}_{\rm
IFR}(t|\lambda_{02},\gamma_{2},\mathbf{\theta}_{2},\mathbf{v}_{2})\right)$
$\displaystyle=$ $\displaystyle-\log\left(\pi e^{-\Lambda_{\rm
DFR}(t|\lambda_{01},\gamma_{1},\mathbf{\theta}_{1},\mathbf{v}_{1})}+(1-\pi)e^{-\Lambda_{\rm
IFR}(t|\lambda_{02},\gamma_{2},\mathbf{\theta}_{2},\mathbf{v}_{2})}\right)$

Simulate as a mixture: draw $c\sim\text{Cat}(\\{1,2\\};(\pi,1-\pi))$ then if
$c=1$ draw from the DFR component (parameters
$(\lambda_{01},\gamma_{1},\mathbf{\theta}_{1},\mathbf{v}_{1})$), and if $c=2$
draw from the IFR component (parameters
$(\lambda_{02},\gamma_{2},\mathbf{\theta}_{2},\mathbf{v}_{2})$).

The priors for the two baseline rates $\lambda_{01}$ and $\lambda_{02}$ can be
taken to be $\text{Exp}(\nu/\gamma_{1})$ and $\text{Exp}(\nu/\gamma_{2})$
respectively.

#### Model 6. Log convex (LCV)

The hazard and cumulative hazard rate functions implied by (12) are:

$\displaystyle\lambda(t|\lambda_{0},w_{0},\gamma,\mathbf{\theta},\mathbf{v})$
$\displaystyle=$
$\displaystyle\lambda_{0}\exp\left(w_{0}t+\sum_{k=1}^{K}w_{k}\text{max}(0,t-\theta_{k})\right)$
(36)
$\displaystyle\Lambda(t|\lambda_{0},w_{0},\gamma,\mathbf{\theta},\mathbf{v})$
$\displaystyle=$
$\displaystyle\sum_{\ell=0}^{k^{\ast}(t)-1}\lambda_{0}e^{-D^{\ast}_{\ell}}\left\\{\frac{e^{C^{\ast}_{0\ell}\theta^{\ast}_{\ell+1}}-e^{C^{\ast}_{0\ell}\theta^{\ast}_{\ell}}}{C^{\ast}_{0\ell}}\right\\}$
(37)
$\displaystyle+\lambda_{0}e^{-D^{\ast}_{k^{\ast}(t)}}\left\\{\frac{e^{C^{\ast}_{0k^{\ast}(t)}t}-e^{C^{\ast}_{0k^{\ast}(t)}\theta^{\ast}_{k^{\ast}(t)}}}{C^{\ast}_{0k^{\ast}(t)}}\right\\}$

where we are using the re-indexed set of weights and locations
$\\{(w^{\ast}_{k},\theta^{\ast}_{k})\\}_{k=1}^{K}$, and where

$\displaystyle
C^{\ast}_{0\ell}=w_{0}+C^{\ast}_{\ell}=\sum_{k=0}^{\ell}w^{\ast}_{k}$ (38)

Note that $C^{\ast}_{0k}$ is non-decreasing in $k$, and can in some
circumstances take the value zero. Where $C^{\ast}_{0k}$ is zero the two
quantities in the braces in the expression for
$\Lambda(t|\lambda_{0},w_{0},\gamma,\theta,{\bf v})$ should be replaced by

$\theta^{\ast}_{k+1}-\theta^{\ast}_{k}$

and

$t-\theta^{\ast}_{k^{\ast}(t)}$

respectively.

Again following the same protocol as in the IFR case we define
$\Lambda^{\ast}_{k}=\Lambda(\theta^{\ast}_{k})$ with $\theta^{\ast}_{0}=0$ and
$\Lambda^{\ast}_{0}=0$. Then in the interval $\theta^{\ast}_{k^{\ast}(t)}\leq
t<\theta^{\ast}_{k^{\ast}(t)}$ the function $\Lambda$ has the exponential form

$\displaystyle\Lambda(t|\lambda_{0},w_{0},\gamma,\mathbf{\theta},\mathbf{v})$
$\displaystyle=$
$\displaystyle\Lambda^{\ast}_{k^{\ast}(t)}+\frac{\lambda_{0}e^{-D^{\ast}_{k^{\ast}(t)}}}{C^{\ast}_{0l^{\ast}}}\left\\{e^{C^{\ast}_{0k^{\ast}(t)}t}-e^{C^{\ast}_{0k^{\ast}(t)}\theta^{\ast}_{k^{\ast}(t)}}\right\\}$
(39)

It follows that If $U\sim\text{Uniform}(0,1)$ then $T$ is a draw from this
distribution if

$\displaystyle
T=\frac{1}{C^{\ast}_{0k^{\ast}}}\log\left[e^{C^{\ast}_{0k^{\ast}}\theta^{\ast}_{k^{\ast}}}+C^{\ast}_{0k^{\ast}}\frac{-\log
U-\Lambda^{\ast}_{k^{\ast}}}{\lambda_{0}e^{-D^{\ast}_{k^{\ast}}}}\right]$ (40)

where
$k^{\ast}=\text{max}\\{k\in\\{1,\ldots,K\\}\,:\,\Lambda^{\ast}_{k}\leq-\log
U\\}$. In the case where $C^{\ast}_{0k^{\ast}}=0$ instead we have:

$T=\theta^{\ast}_{k^{\ast}}+\frac{-\log
U-\Lambda^{\ast}_{k^{\ast}}}{\lambda_{0}e^{-D^{\ast}_{k^{\ast}}}}$

Priors are needed for $\lambda_{0}$ and $w_{0}$. Suitable priors are

$\displaystyle\log\lambda_{0}$ $\displaystyle\sim$
$\displaystyle\text{Normal}(0,(\gamma/\nu)^{2})$ (41) $\displaystyle w_{0}$
$\displaystyle\sim$ $\displaystyle\text{Normal}(0,(\gamma/\nu)^{2})$ (42)

##  4 Demonstration

We now demonstrate the properties of the six models specified in Section 2.
Figures 3-8 show example draws from each of the models, and in each case we
show the location/weight pairs $\\{(\theta_{k},w_{k})\\}$, the hazard rate
function $\lambda(t)$, the cumulative hazard rate $\Lambda(t)$, the density
function $f(t)=\lambda(t)\exp(-\Lambda(t))$ and the survival function
$\bar{F}(t)=\exp(-\Lambda(t))$. We also show a histogram of $n=1000$ random
draws from the failure time distribution.

We note that any draw from the Gamma Process prior is concentrated on a set of
discrete points, and that the weights $w_{k}$ decrease stochastically in $k$.
This is what enables the truncation of the otherwise infinite sums in §2.3.
Consider a draw $G$ from $G\sim\text{GaPP}_{K}(\alpha H_{0},\beta)$ as given
in (15). We set $\alpha=3$, $\beta=1$ and choose the baseline distribution
$H_{0}$ as $\text{Exp}(\phi)$ with $\phi=1$. Figure 2 shows the first $K=100$
location/weight pairs $\\{(\theta_{k},w_{k})\\}_{k=1}^{K}$. The first 4
weights (labelled individually in Figure 2) account for 87.2% of the total
weight, and the first 40 account for all but 0.0001% of the total weight.

Figure 2: A draw from the Gamma Process Prior. The locations and weights are
plotted as filled circles, except the first four points $k=1,2,3,4$ are
labelled and shown as white circles.

The discrete nature of the draw $G$ is visible in Figures 3-8 as multiple step
changes in the hazard rate function $\lambda(t)$, and in sharp peaks and
troughs in the failure time distribution $f(t|\lambda(\cdot))$.

To aid comparison we have used the same single draw $G$ from Figure 2 for the
IFR, DFR, LWB and LCV models, and have used this same draw (though rescaled)
as $G_{1}$ in the SBT and MBT models. The draw $G_{2}$ in the SBT and MBT
models has the same $\alpha=3$ and $\beta=1$ parameter values, but the
baseline failure rate distribution $H_{0}$ is a Normal$(2,1^{2})$.

The Increasing Failure Rate case (Model 1, IFR, Figure 3) is a simple case
where the hazard rate $\lambda(t)$ increases by the corresponding weight value
at each value of $\theta$. The hazard rate thus jumps upwards at these
discrete points and is constant between them. The point with the largest
weight in Figure 3(a) is responsible for the single largest jump in
$\lambda(t)$, and corresponds to the sharp peak in the failure time density
$f(t)$ in Figure 3(d), and also to the sudden sharp decline in the survival
function $\bar{F}(t)$ in Figure 3(e).

Figure 3: Model 1. Increasing Failure Rate (IFR): A single draw from the
prior, and a sample of simulated data. (a) Locations and weights
$\\{(\theta_{k},w_{k})\\}$; (b) Hazard rate function $\lambda(t)$; (c)
Cumulative hazard rate function $\Lambda(t)$; (d) Probability density $f(t)$;
(e) Survival Function $\bar{F}(t)$; (f) Sample of $n=1000$ failure time
observations.

The Decreasing Failure Rate case (Model 2, DFR, Figure 4) shows a hazard rate
function $\lambda(t)$ the opposite behaviour to the IFR case: namely at each
weight location the hazard rate jumps downwards. With the hazard rate so large
at early times the failure time distribution is more strongly concentrated
towards zero than in the IFR case.

Figure 4: Model 2. Decreasing Failure Rate (DFR): A single draw from the
prior, and a sample of simulated data. (a) Locations and weights
$\\{(\theta_{k},w_{k})\\}$; (b) Hazard rate function $\lambda(t)$; (c)
Cumulative hazard rate function $\Lambda(t)$; (d) Probability density $f(t)$;
(e) Survival Function $\bar{F}(t)$; (f) Sample of $n=1000$ failure time
observations.

For the Lo-Weng Bathtub (Model 3, LWB, Figure 5) the hazard rate function has
reflectional symmetry about $t=a$, which we have set to be at $a=0.6$. There
are downward (DFR) jumps of $w_{k}$ at times $a-\theta_{k}$ and equal upward
(IFR) jumps at times $a+\theta_{k}$. This bathtub hazard rate is responsible
for the bimodal failure time distribution $f(t)$ in Figure 5(d).

Figure 5: Model 3. Lo-Weng Bathtub (LWB): A single draw from the prior, and a
sample of simulated data. (a) Locations and weights
$\\{(\theta_{k},w_{k})\\}$; (b) Hazard rate function $\lambda(t)$; (c)
Cumulative hazard rate function $\Lambda(t)$; (d) Probability density $f(t)$;
(e) Survival Function $\bar{F}(t)$; (f) Sample of $n=1000$ failure time
observations. The location of the inflection point $a$ is shown by vertical
dashed lines.

The Superposition Bathtub (Model 4, SBT, Figure 6) has a hazard rate function
that is the linear combination of a DFR and IFR hazard rate function. In
Figure 6(b) the sum of these two monotonic hazard rate functions leads to the
bathtub shape of $\lambda(t)$, and again to the bimodal failure time
distribution $f(t)$. Figure 6(a) shows the two separate sets of weights and
locations for the two draws
$G_{1}\sim\text{GaPP}_{K}(\alpha_{1}H_{01},\beta_{1})$ and
$G_{2}\sim\text{GaPP}_{K}(\alpha_{2}H_{02},\beta_{2})$ where $H_{01}$ is an
exponential distribution and $H_{02}$ is a Normal distribution.

Figure 6: Model 4. Superposition Bathtub (SBT): A single draw from the prior,
and a sample of simulated data. (a) Locations and weights
$\\{(\theta_{k},w_{k})\\}$ with the DFR component weights shown as open
circles, and the IFR component as ’+’ symbols; (b) Hazard rate function
$\lambda(t)$; (c) Cumulative hazard rate function $\Lambda(t)$; (d)
Probability density $f(t)$; (e) Survival Function $\bar{F}(t)$; (f) Sample of
$n=1000$ failure time observations.

As we noted earlier, the Mixture Bathtub (Model 5, MBT, Figure 7) case does
not have an overall bathtub hazard function, but does show the characteristic
U-shape at early times (a decrease followed by an increase). This leads to the
bimodal failure time distribution $f(t)$. In this example the two components
have been given equal weight in the mixture: $\pi=0.5$.

Figure 7: Model 5. Mixture Bathtub (MBT): A single draw from the prior, and a
sample of simulated data. (a) Locations and weights $\\{(\theta_{k},w_{k})\\}$
with the DFR component weights shown as open circles, and the IFR component as
’+’ symbols; (b) Hazard rate function $\lambda(t)$; (c) Cumulative hazard rate
function $\Lambda(t)$; (d) Probability density $f(t)$; (e) Survival Function
$\bar{F}(t)$; (f) Sample of $n=1000$ failure time observations.

The Log Convex case (Model 6, LCV, Figure 8) differs somewhat from the others
due to the random measure $G$ contributing linearly to the logarithm of the
hazard rate function. The corresponding hazard rate function $\lambda(t)$ is
piecewise curved rather than being piecewise constant. The example shown in
Figure 8 shows a failure time distribution with a monotonically decreasing
rather than bimodal density. With appropriate choice of parameters this
monotonic behaviour can be seen in all of Models 3-6, and is a consequence of
a large integrated hazard at early times making it less likely that failures
occur at late times. The increasing hazard in such situations serves to
shorten the tail of $f(t)$, rather than creating a mode at late times.

Figure 8: Model 6. Log Convex (LCV): A single draw from the prior, and a
sample of simulated data. (a) Locations and weights
$\\{(\theta_{k},w_{k})\\}$; (b) Hazard rate function $\lambda(t)$; (c)
Cumulative hazard rate function $\Lambda(t)$; (d) Probability density $f(t)$;
(e) Survival Function $\bar{F}(t)$; (f) Sample of $n=1000$ failure time
observations.

##  5 Conclusions

This paper lays out the details of simulation for four different bathtub
hazard rate functions (Models 3-6), based on draws from the Gamma Process
Prior. These non-parametric specifications are highly flexible, and include
only minimal assumptions about the form of the hazard rate function.

The examples in Section 4 demonstrate the range of behaviours that the various
models can exhibit, and we expect to be able to find suitable models to match
many real data sets.

The next step is thus inference from data, which we intend to address in a
future paper. Arnold, Chukova, and Hayakawa (2020) have demonstrated how
inference can be carried out for the IFR model using MCMC samplers in a fully
Bayesian framework. Some of the parameter updates are Gibbs however updates of
the weights $\\{w_{k}\\}$ and locations $\\{\theta_{k}\\}$ require Metropolis-
Hastings proposals. Arnold, Chukova, and Hayakawa (2020) noted that their
samplers may need further development to ensure that they converge
efficiently, and made various suggestions for improvements which may be
especially relevant in the case of censored data.

We anticipate that the model selection question, choosing which among these
models is best, may be carried out using a criterion such as WAIC (Watanabe
2010), with model checking using posterior predictive distributions.

### Acknowledgements

This work was supported by: Waseda University, Grant for Special Research
Projects (2018K-383); JSPS KAKENHI Grant-in-Aid for Scientific Research (C)
Grant Number 18K04621; Waseda Institute for Advanced Study Visiting Scholars
2018; FY2018 and FY2022 Grant Program for Promotion of International Joint
Research, Waseda University. Fulbright New Zealand: Fulbright Scholar Award
2018.

## References

reAbbasi, Jamal N. Al, Mundher A. Khaleel, Moudher Kh. Abdal-hammed, Yue Fang
Loh, and Gamze Ozel. 2019. “A new uniform distribution with bathtub-shaped
failure rate with simulation and application.” _Mathematical Science_ 13:
105–14.

preArnold, Richard, Stefanka Chukova, and Yu Hayakawa. 2020. “Nonparametric
Bayesian Analysis of Hazard Rate Functions using the Gamma Process Prior.” In
_2020 Asia-Pacific International Symposium on Advanced Reliability and
Maintenance Modeling (APARM)_.

preDykstra, R. L., and Purushottam Laud. 1981. “A Bayesian Non-parametric
Approach to Reliability.” _Annals of Statistics_ 9: 356–67.

preFerguson, Thomas S. 1974. “Prior Distributions on Spaces of Probability
Measures.” _The Annals of Statistics_ 2: 615–29.

preGlaser, Ronald E. 1980. “Bathtub and Related Failure Rate
Characterizations.” _Ournal of the American Statistical Association_ 75:
667–72.

preHayakawa, Yu, Jonathan Zukerman, Sue Paul, and Tony Vignaux. 2001.
“Bayesian non-parametric testing of constant versus non-decreasing hazard
rates.” In _System and Bayesian Reliability_ , edited by Yu Haykawa, Telba Z.
Irony, and Min Xie, 5:391–406. Series on Quality, Reliability and Engineering
Statistics. Singapore: World Scientific.

preHo, Man-wai, and A. Y. Lo. 2001. “Bayesian non-parametric estimation of a
monotone hazard rate.” In _System and Bayesian Reliability_ , edited by Yu
Haykawa, Telba Z. Irony, and Min Xie, 5:301–14. Series on Quality, Reliability
and Engineering Statistics. Singapore: World Scientific.

preIqbal, Alya Al Al MutairiMuhammad Zafar, Muhammad Zafar Iqbal, Zeeshan
Arshad, Badr Alnssyan, Hazem Al-Mofleh, and Ahmed Z. Afify. 2021. “A New
Extended Model with Bathtub Shaped Failure Rate: Properties, Inference,
Simulation, and Applications.” _Mathematics_ 9: 2024:1–32.
https://doi.org/10.3390/math9172024.

preKaplan, E. L., and Paul Meier. 1958. “Nonparametric Estimation from
Incomplete Observations.” _Journal of the American Statistical Association_
53: 457–81.

preKingman, J. F. C. 1967. “Completely random measures.” _Pacific Journal of
Mathematics_ 21: 59–78.

preKumar, U., B. Klefsjö, and S. Granholm. 1989. “Reliability investigation
for a fleet of load haul dump machines in a Swedish mine.” _Reliability
Engineering and System Safety_ 26: 341–61.

preLo, A. Y., L. J. Brunner, and A. T. Chan. 1998. “Weighted Chinese
restaurant processes and Bayesian mixture models (Revision 1.1).” Department
of Information; Systems Management, Hong Kong University of Science;
Technology, Hong Kong: Hong Kong University of Science; Technology.

preLo, A. Y., and C-S. Weng. 1989. “On a class of Bayesian non-parametric
estimates II. Hazard rate estimates.” _Annals of the Institute of Statistical
Mathematics_ 41: 221–45.

preMüller, Peter, Fernando Andrés Quintana, Alejandro Jara, and Tim Hanson.
2015. _Bayesian Nonparametric Data Analysis_. New York, NY: springer.

prePaisley, J. 2010. “A Simple Proof of the Stick-Breaking Construction of the
Dirichlet Process.” Department of Computer Science, Princeton University,
Princeton, NJ: Princeton University.
http://www.columbia.edu/~jwp2128/Teaching/E6892/papers/SimpleProof.pdf.

prePaisley, J., D. M. Blei, and M. I. Jordan. 2012. “Stick-Breaking Beta
Processes and the Poisson Process.” In _Proceedings of the 15th International
Conference on Artificial Intelligence and Statistics (AISTATS) 2012, La Palma,
Canary Islands_ , 850–58.

prePaisley, J., A. Zaas, C. W. Woods, G. S. Ginsburg, and L. Carin. 2010. “A
Stick-Breaking Construction of the Beta Process.” In _Proceedings of the 26th
International Conference on Machine Learning, Haifa, Israel, 2010_ , 1–8.

prePeng, Chong, Guangpeng Liu, and Lun Wang. 2016. “Piecewise modelling and
parameter estimation of repairable system failure rate.” _SpringerPlus_ 5:
1477:1–14.

prePhadia, Eswar G. 2016. _Prior Processes and Their Applications:
Nonparametric Bayesian Estimation_. 2nd ed. New York, NY: Springer.

preRoychowdhury, Anirban, and Brian Kulis. 2014. “Gamma Processes, Stick-
Breaking, and Variational Inference.” https://arxiv.org/abs/1410.1068.

pre———. 2015. “Gamma Processes, Stick-Breaking, and Variational Inference.” In
_Proceedings of the 18th International Conference on Artificial Intelligence
and Statistics (AISTATS)_ , 800–808.

preSethuraman, Jayaram. 1994. “A constructive definition of Dirichlet priors.”
_Statistica Sinica_ 4: 639–50.

preWatanabe, Sumio. 2010. “Asymptotic Equivalence of Bayes Cross Validation
and Widely Applicable Information Criterion in Singular Learning Theory.”
_Journal of Machine Learning Research_ 11: 3571–94.

p