# State estimation for aoristic models

M.N.M. van Lieshout  Centrum Wiskunde & Informatica (CWI), P.O. Box 94079,
1090 GB Amsterdam, The Netherlands  Department of Applied Mathematics,
University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands  R.L.
Markwitz  Centrum Wiskunde & Informatica (CWI), P.O. Box 94079, 1090 GB
Amsterdam, The Netherlands  Department of Applied Mathematics, University of
Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands

ABSTRACT: Aoristic data can be described by a marked point process in time in
which the points cannot be observed directly but are known to lie in
observable intervals, the marks. We consider Bayesian state estimation for the
latent points when the marks are modelled in terms of an alternating renewal
process in equilibrium and the prior is a Markov point point process. We
derive the posterior distribution, estimate its parameters and present some
examples that illustrate the influence of the prior distribution.

KEYWORDS: alternating renewal process, aoristic data, marked temporal point
process, Markov chain Monte Carlo methods, Markov point process state
estimation.

##  1 Introduction

Inference for point processes on the real line has been dominated by a dynamic
approach based on the stochastic intensity or hazard function (Brémaud, 1972;
Karr, 1991; Last & Brandt, 1995). Such an approach is quite natural, is
amenable to likelihood-based inference and allows the utilisation of powerful
tools from martingale theory. However, it breaks down completely when
censoring breaks the orderly progression of time. In such cases, state
estimation techniques are needed that are able to fill in the gaps (Brix &
Diggle, 2001; Van Lieshout, 2016).

In this paper, we concentrate on aoristic data (Ratcliffe & McCullagh, 1998)
in which the points may not be observed directly but upper and lower bounds
exist. Such data are commonplace in criminology. Suppose, for example, that a
working person leaves his place of residence early in the morning and returns
late in the afternoon to discover that the residence has been burgled. Then
the exact time of the break-in cannot be determined, but it must have happened
during the absence of the resident. In rare cases, a burglar may also be
caught in the act, in which case the time of break-in coincides with the
bounds. The actual break-in times tend to be estimated by ad hoc, naive
approaches, e.g. the mid-point of the reported interval (Helms, 2008) or the
length-weighted empirical probability mass function of the interval lengths.
An obvious disadvantage of such methods is that dependencies between offence
times, such as the near-repeat effect (Bernasco, 2009), are ignored.

The focus of this research is to develop a Bayesian inference framework for
aoristic data that is able to infer missing information and takes into account
expert knowledge and interaction. Specifically, in Section 2 we formalise
aoristic censoring as a marked point process in time in which the points
cannot be observed directly but are known to lie in observable intervals, the
marks. Upon employing a Markov point process prior (Van Lieshout, 2000), the
posterior distribution of the point locations is derived in Section 3. In
Section 4 we turn to Monte Carlo based inference. The paper is concluded by
simulated examples that demonstrate the influence of the prior (Section 5) and
some reflections on future research.

##  2 Marked point process formulation

###  2.1 Alternating renewal processes for censoring

In this paper, we consider a censoring mechanism based on an alternating
renewal process. Let $C_{1},C_{2},\,...$ be a sequence of random $2$-vectors
such that $C_{i}=(Y_{i},Z_{i})$, $i\in{\mathbb{N}}$, are independent and
identically distributed (Asmussen, 2003; Ross, 1996). Furthermore, assume that
$C_{i}$ has joint probability density function $f$ on
$({\mathbb{R}}^{+})^{2}$. Introduce $T_{i}$ = $Y_{i}+Z_{i}$, set $S_{0}=0$ and
let, for $n\in{\mathbb{N}}$, $S_{n}=\sum_{i=1}^{n}T_{i}$ be the time of the
$n$th renewal. Note that no renewal occurs at the end of a $Y$-phase.
Furthermore, assume that $0<{\mathbb{E}}T_{1}<\infty$. Then, by the strong law
of large numbers,

$N(t)=\sup\left\\{n\in{\mathbb{N}}_{0}:S_{n}\leq t\right\\},\quad t\geq 0,$
(2.1)

is well-defined and the supremum is attained with probability one.
Furthermore, the renewal function

$M(t)={\mathbb{E}}N(t)=\sum_{n=1}^{\infty}{\mathbb{P}}(S_{n}\leq t),\quad
t\geq 0,$

is finite and absolutely continuous with respect to Lebesgue measure (Ross,
1996, Chapter 3).

$S_{i}$$S_{i+1}$$Y_{i+1}$$Z_{i+1}$$(1)$$t$$S_{i}$$S_{i+1}$$t$$Y_{i+1}$$Z_{i+1}$$(2)$
Figure 1: View of two alternating renewal processes. In $(1)$, a point $t$ of
the point process $X$ falls in a $Y$-phase. In criminological terms, this may
correspond to a person being away from home and being burgled; the exact time
of burglary is unknown. Hence the whole interval is recorded. In $(2)$, a
point $t$ of the point process $X$ falls in a $Z$-phase. This corresponds to a
person being home and being burgled. In this case, the exact time of burglary
is known, thus the exact point is recorded. Renewal times are denoted by
$S_{i},\,S_{i+1}$ with $S_{0}=0$.

An alternating renewal process can be used for censoring in the following way.
Let $X$ be a temporal point process on ${\mathbb{R}}^{+}$ (Brémaud, 1972;
Daley & Vere-Jones, 2003, 2008; Last & Brandt, 1995) and let each point $t\in
X$ be associated with an alternating renewal process, independently of other
points of $X$. Now, if $t$ happens to fall in some $Z$-phase, $t$ is observed
perfectly, whereas $t$ is observed aoristically if it falls in a $Y$-phase.
The censoring mechanism is illustrated in Figure 1.

###  2.2 Age and excess distribution

Aoristic data generated by the censoring mechanism described in Section 2.1
can be expressed in terms of the age and excess (also referred to as residual
lifetime) with respect to the $Y$-process. Recall that for an alternating
renewal process, the age with respect to the $Y$-process is defined as

$\displaystyle A(t)$
$\displaystyle=(t-S_{N(t)})\,{\mathbf{1}}\\{S_{N(t)}+Y_{N(t)+1}>t\\},\quad
t\geq 0,$

the excess with respect to the $Y$-process as

$\displaystyle B(t)$
$\displaystyle=(S_{N(t)+1}-Z_{N(t)+1}-t)\,{\mathbf{1}}\\{S_{N(t)}+Y_{N(t)+1}>t\\},\quad
t\geq 0.$

Let us parametrise an interval on the real line by its left-most point and
length. In other words, a pair $(a,l)\in{\mathbb{R}}\times{\mathbb{R}}^{+}$
corresponds to the closed interval $[a,a+l]$. Then, the recorded interval of a
latent point $t\geq 0$ can be written as $t+[-A(t),B(t)]$. Note that
$[-A(t),B(t)]$ is parametrised by
$I(t)=(-A(t),A(t)+B(t))\in{\mathbb{R}}\times{\mathbb{R}}^{+}$.

First, let us consider the joint distribution of age and excess with respect
to the $Y$-process.

######  Proposition 2.1.

Let $N$ be an alternating renewal process as in (2.1). Assume that $T_{1}$ is
absolutely continuous with respect to Lebesgue measure and that
$0<{\mathbb{E}}T_{1}<\infty$. Then, for $t\geq 0$, the joint distribution of
$(A(t),B(t))$ has an atom at $(0,0)$ of size

$c(t)=F_{Y}(t)-\int_{0}^{t}\left[1-F_{Y}(t-s)\right]dM(s)$

and, for $0\leq u\leq t,v\geq 0$,

${\mathbb{P}}(A(t)\leq u;B(t)\leq
v)=c(t)+\left[F_{Y}(t+v)-F_{Y}(t)\right]{\mathbf{1}}\\{u=t\\}+\int_{t-u}^{t}\left[F_{Y}(t+v-s)-F_{Y}(t-s)\right]dM(s).$

Here $F_{Y}$ denotes the cumulative distribution function of $Y_{1}$ and $M$
is the renewal function.

###### Proof..

Write $F_{n}$ for the cumulative distribution function of $S_{n}$,
$n\in{\mathbb{N}}$. By partitioning over the number of renewals up to time $t$
and upon noting that $N(t)=n$ if and only if $S_{n}\leq t$ and
$S_{n}+Y_{n+1}+Z_{n+1}>t$, one obtains that

$\displaystyle{\mathbb{P}}(A(t)\leq u;B(t)\leq v)$ $\displaystyle=$
$\displaystyle c(t)+{\mathbb{P}}(t-S_{N(t)}\leq u;\,t<S_{N(t)}+Y_{N(t)+1}\leq
t+v)$ $\displaystyle=$ $\displaystyle
c(t)+\sum_{n=0}^{\infty}{\mathbb{P}}(t-u\leq
S_{N(t)};\,t<S_{N(t)}+Y_{N(t)+1}\leq t+v;\,N(t)=n)$ $\displaystyle=$
$\displaystyle c(t)+{\mathbb{P}}(t-u\leq S_{N(t)};\,t<S_{N(t)}+Y_{N(t)+1}\leq
t+v;\,N(t)=0)$ $\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\int_{t-u}^{t}{\mathbb{P}}(t-s<Y_{n+1}\leq
t+v-s)\,dF_{n}(s).$

The claim follows by an application of Fubini’s theorem for the last term, the
observation that

${\mathbb{P}}(t-u\leq S_{N(t)};\,t<S_{N(t)}+Y_{N(t)+1}\leq
t+v;\,N(t)=0)={\mathbb{P}}(t<Y_{1}\leq t+v)$

if $u=t$ and zero otherwise, and because

$\displaystyle c(t)$ $\displaystyle=$ $\displaystyle
1-{\mathbb{P}}(S_{N(t)}+Y_{N(t)+1}>t)$ $\displaystyle=$ $\displaystyle
1-{\mathbb{P}}(Y_{1}>t)-\sum_{n=1}^{\infty}\int_{0}^{t}{\mathbb{P}}(Y_{n+1}>t-s)\,dF_{n}(s)=F_{Y}(t)-\int_{0}^{t}\left[1-F_{Y}(t-s)\right]dM(s).$

∎

The long-run behaviour as time goes to infinity can be obtained by appealing
to the key renewal theorem (for example found in Ross (1996, Theorem 3.4.2)).
Specialising to the parameter vector $I(t)$, the following theorem holds.

######  Theorem 2.2.

Let $N$ be an alternating renewal process as in (2.1). Assume that $Y_{1}$ and
$T_{1}$ are absolutely continuous with respect to Lebesgue measure with
probability density functions $f_{Y}$ and $f$ and that
$0<{\mathbb{E}}T_{1}<\infty$. Then $(-A(t),A(t)+B(t))$ tends in distribution
to $\nu$, the mixture of an atom at $(0,0)$ and an absolutely continuous
component that has probability density function $f_{Y}(l)/{\mathbb{E}}Y_{1}$
on $\\{(a,l)\in{\mathbb{R}}\times{\mathbb{R}}^{+}:a\leq 0\leq a+l\\}$. The
mixture weights are, respectively, ${\mathbb{E}}Z_{1}/{\mathbb{E}}T_{1}$ and
${\mathbb{E}}Y_{1}/{\mathbb{E}}T_{1}$.

###### Proof..

First, let us consider the limit behaviour of the joint cumulative
distribution function of $A(t)$ and $B(t)$ as $t\to\infty$. With the notation
of Proposition 2.1, by Theorem 3.4.4 in (Ross, 1996), $c(t)$ converges to
${\mathbb{E}}Z_{1}/{\mathbb{E}}T_{1}$. Also, for $t>u$, the second term in the
joint cumulative distribution function of $A(t)$ and $B(t)$ is zero. For the
last term, note that for $v\geq 0$ the function
$h_{v}:{\mathbb{R}}^{+}\to{\mathbb{R}}$ defined by $h_{v}(s)=1-F_{Y}(v+s)$ is
non-negative, monotonically non-increasing and integrable. Hence the key
renewal theorem implies that

$\lim_{t\to\infty}\int_{0}^{t}\left[1-F_{Y}(t+v-s)\right]dM(s)=\frac{1}{{\mathbb{E}}T_{1}}\int_{0}^{\infty}\left[1-F_{Y}(v+s)\right]ds.$

Analogously, for fixed $u\geq 0$, $t-u\to\infty$ if and only if $t\to\infty$.
Writing $s=t-u$,

$\lim_{s\to\infty}\int_{0}^{s}\left[1-F_{Y}(s+u+v-r)\right]dM(r)=\frac{1}{{\mathbb{E}}T_{1}}\int_{0}^{\infty}\left[1-F_{Y}(u+v+r)\right]dr$
$=\lim_{t\to\infty}\int_{0}^{t-u}\left[1-F_{Y}(t+v-r)\right]dM(r).$

We conclude that $Q(u,v)=\lim_{t\to\infty}{\mathbb{P}}(A(t)\leq u;B(t)\leq v)$
exists and equals

$Q(u,v)=\frac{{\mathbb{E}}Z_{1}}{{\mathbb{E}}T_{1}}+\frac{1}{{\mathbb{E}}T_{1}}\int_{0}^{u}\left[F_{Y}(v+s)-F_{Y}(s)\right]ds.$

Note that $Q$ is the cumulative distribution function of the mixture of an
atom at $(0,0)$ and an absolutely continuous component with probability
density function $f_{Y}(u+v)/{\mathbb{E}}Y_{1}$ on $({\mathbb{R}}^{+})^{2}$.
By Helly’s continuity theorem, $(A(t),B(t))$ converges in distribution.

Turning to the parametrisation $I(t)=(-A(t),A(t)+B(t))$, its limit
distribution inherits an atom at $(0,0)$ from $Q$ of size
${\mathbb{E}}Z_{1}/{\mathbb{E}}T_{1}$. By the change of variables bijection
$h:({\mathbb{R}}^{+})^{2}\to{\mathbb{R}}^{-}\times{\mathbb{R}}^{+}$ defined by
$h(u,v)=(-u,u+v)$, since $h$ is differentiable, the absolutely continuous part
has probability density function
$f_{Y}(h^{-1}(a,l))\,|\text{det}\,J_{h^{-1}}(a,l)|/{\mathbb{E}}Y_{1}=f_{Y}(-a+a+l)/{\mathbb{E}}Y_{1}$,
where $J_{h^{-1}}$ is the Jacobian of $h^{-1}$. ∎

###  2.3 Complete model formulation

We are now ready to formulate a model. Let ${\cal{X}}$ be an open subset of
the positive half-line ${\mathbb{R}}^{+}$. The state space of $X$, denoted by
${\cal{N}}_{\cal{X}}$, consists of finite sets
$\\{t_{1},\dots,t_{n}\\}\subset{\cal{X}}$, $n\in{\mathbb{N}}_{0}$, which we
equip with the Borel $\sigma$-algebra of the weak topology (Daley & Vere-
Jones, 2003, 2008). We will assume that the distribution of $X$ is specified
in terms of a probability density function $p_{X}$ with respect to the
distribution of a unit rate Poisson process on ${\cal{X}}$ (Van Lieshout,
2000).

Upon labelling the points of $X$ independently with a mark according to the
mixture distribution of Theorem 2.2, denoted by $\nu$, the complete model $W$
is obtained. Its realisations are sets
$\\{(t_{1},I_{1}),\dots,(t_{n},I_{n})\\}\subset{\cal{X}}\times({\mathbb{R}}\times{\mathbb{R}}^{+})$.
For parametrisation $I_{j}=(a_{j},l_{j})$, the pair $(t_{j},I_{j})$ defines an
interval $[t_{j}+a_{j},t_{j}+a_{j}+l_{j}]$. The ensemble of all realisations
is denoted by
${\cal{N}}_{{\cal{X}}\times({\mathbb{R}}\times{\mathbb{R}}^{+})}$ and equipped
with the Borel $\sigma$-algebra of the weak topology. Note that $W$ has
probability density function $p_{X}$ with respect to the distribution of a
Poisson process on ${\cal{X}}\times({\mathbb{R}}\times{\mathbb{R}}^{+})$ with
intensity measure $\ell\times\nu$ where $\ell$ is Lebesgue measure.

Due to the censoring, one does not observe the complete model $W$ but rather
the set

$U=\bigcup_{(t,I)\in W}((t,0)+I).$ (2.2)

Our aim is to reconstruct $X$ or $W$ from $U$. In order to do so, the
posterior distribution of $X$ or $W$ given $U$ is needed. This will be the
topic of the next section.

##  3 The Bayesian framework

In a Bayesian framework, the posterior distribution updates prior forms in the
light of data gathered (Gamerman & Lopes, 2006). Heuristically,

$\displaystyle p_{X|U}(\mathbf{x}\,|\,\mathbf{u})\propto
p_{U|X}(\mathbf{u}\,|\,\mathbf{x})\,\,p_{X}(\mathbf{x})$ (3.1)

through the use of Bayes’ theorem. The term
$p_{U|X}({\mathbf{u}}\,|\,{\mathbf{x}})$ describes the likelihood that the
points of ${\mathbf{x}}$ generate the intervals in ${\mathbf{u}}$. In the
literature this term is referred to as a forward term, forward density or
forward model (Van Lieshout, 1995; Van Lieshout & Baddeley, 2002). The term
$p_{X}({\mathbf{x}})$ captures prior beliefs about the geometry of
${\mathbf{x}}$. In our context, since the forward model is a mixture of
discrete and absolutely continuous components, some care is required in
handling (3.1).

######  Theorem 3.1.

Let $W$ be a point process on the open set ${\cal{X}}\subset{\mathbb{R}}$ with
probability density function $p_{X}$ with respect to the distribution of a
unit rate Poisson process on ${\cal{X}}$ marked independently with mark
distribution $\nu$ defined in Theorem 2.2 . Write $X$ for the ground process
of locations in ${\cal{X}}$ and consider the forward model (2.2). Let
${\mathbf{u}}$ be a realisation of $U$ that consists of an atomic part
$\\{(a_{1},0),\dots,(a_{m},0)\\}$, $m\in{\mathbb{N}}_{0}$, and a non-atomic
part $\\{(a_{m+1},l_{m+1}),\,\dots,(a_{n},l_{n})\\}$, $n\geq m$. Then the
posterior distribution of $X$ given $U={\mathbf{u}}$ satisfies, for $A$ in the
Borel $\sigma$-algebra of the weak topology on ${\cal{N}}_{\cal{X}}$,

${\mathbb{P}}(X\in A\mid
U={\mathbf{u}})=c({\mathbf{u}})\int_{{\cal{X}}^{n-m}}p_{X}(\\{a_{1},\dots,a_{m},x_{1},\dots,x_{n-m}\\}){\mathbf{1}}_{A}(\\{a_{1},\dots,a_{m},x_{1},\dots,x_{n-m}\\})\times$
$\times\left(\sum_{\begin{subarray}{c}D_{1},\dots,D_{n-m}\\\
\cup_{i}\\{D_{i}\\}=\\{1,\dots,n-m\\}\end{subarray}}\prod_{i=1}^{n-m}{\mathbf{1}}\\{x_{D_{i}}\in[a_{m+i},a_{m+i}+l_{m+i}]\\}\right)\prod_{i=1}^{n-m}dx_{i}$

provided that $c({\mathbf{u}})^{-1}$ defined by

$\int_{{\cal{X}}^{n-m}}p_{X}(\\{a_{1},\dots,a_{m},x_{1},\dots,x_{n-m}\\})\left(\sum_{\begin{subarray}{c}D_{1},\dots,D_{n-m}\\\
\cup_{i}\\{D_{i}\\}=\\{1,\dots,n-m\\}\end{subarray}}\prod_{i=1}^{n-m}{\mathbf{1}}\\{x_{D_{i}}\in[a_{m+i},a_{m+i}+l_{m+i}]\\}\right)\prod_{i=1}^{n-m}dx_{i}$

exists in $(0,\infty)$.

###### Proof..

We must show that for each $A$ in the Borel $\sigma$-algebra of
${\cal{N}}_{\cal{X}}$ with respect to the weak topology and each $F$ in the
Borel $\sigma$-algebra of the weak topology on
${\cal{N}}_{{\mathbb{R}}\times{\mathbb{R}}^{+}}$ the following identity holds:

${\mathbb{E}}\left[{\mathbf{1}}_{F}(U)\,{\mathbb{P}}(X\in A\mid
U)\right]={\mathbb{E}}\left[{\mathbf{1}}_{F}(U)\,{\mathbf{1}}_{A}(X)\right].$
(3.2)

Let

$\displaystyle
q_{x}(a,l)=\frac{f_{Y}(l)}{{\mathbb{E}}Y_{1}}\,{\mathbf{1}}\\{a\leq x\leq
a+l\\}$

describe a probability density function for parametrisations of intervals
generated by $x\in{\cal{X}}$, noting that it is jointly measurable as a
function on ${\cal{X}}\times({\mathbb{R}}\times{\mathbb{R}}^{+})$. Then,
denoting the cardinality of a set by $|\cdot|$ and Lebesgue measure by $\ell$,
${\mathbb{E}}[{\mathbf{1}}_{F}(U)\,{\mathbf{1}}_{A}(X)]$ can be expanded as

$\sum_{n=0}^{\infty}\frac{e^{-\ell({\cal{X}})}}{n!}\sum_{C_{0}\subset\\{1,...,n\\}}\left(\frac{{\mathbb{E}}Z_{1}}{{\mathbb{E}}T_{1}}\right)^{|C_{0}|}\left(\frac{{\mathbb{E}}Y_{1}}{{\mathbb{E}}T_{1}}\right)^{n-|C_{0}|}\int_{{\cal{X}}^{n}}{\mathbf{1}}_{A}(\\{x_{1},\dots,x_{n}\\})\,p_{X}(\\{x_{1},...,x_{n}\\})\sum_{\begin{subarray}{c}C_{1},...,C_{n-|C_{0}|}\\\
\cup_{j}\\{C_{j}\\}=\\{1,...,n\\}\setminus C_{0}\end{subarray}}\,$
$\frac{1}{(n-|C_{0}|)!}\left(\int_{({\mathbb{R}}\times{\mathbb{R}}^{+})^{n-|C_{0}|}}\prod_{j=1}^{n-|C_{0}|}q_{x_{C_{j}}}(u_{j})\,{\mathbf{1}}_{F}(\\{u_{1},\dots,u_{n-|C_{0}|}\\}\cup\\{(x_{k},0):k\in
C_{0}\\})\prod_{j=1}^{n-|C_{0}|}\,du_{j}\right)\prod_{i=1}^{n}\,dx_{i}.$ (3.3)

For the left-hand side of (3.2), expanding as before, we get

$\sum_{n=0}^{\infty}\frac{e^{-\ell({\cal{X}})}}{n!}\sum_{C_{0}\subset\\{1,...,n\\}}\left(\frac{{\mathbb{E}}Z_{1}}{{\mathbb{E}}T_{1}}\right)^{|C_{0}|}\left(\frac{{\mathbb{E}}Y_{1}}{{\mathbb{E}}T_{1}}\right)^{n-|C_{0}|}\int_{{\cal{X}}^{n}}p_{X}(\\{x_{1},...,x_{n}\\})\sum_{\begin{subarray}{c}C_{1},...,C_{n-|C_{0}|}\\\
\cup_{j}\\{C_{j}\\}=\\{1,...,n\\}\setminus
C_{0}\end{subarray}}\,\frac{1}{(n-|C_{0}|)!}$
$\left(\int_{({\mathbb{R}}\times{\mathbb{R}}^{+})^{n-|C_{0}|}}\prod_{j=1}^{n-|C_{0}|}q_{x_{C_{j}}}(u_{j})\,{\mathbb{P}}\left(X\in
A\mid U=\\{u_{1},\dots,u_{n-|C_{0}|}\\}\cup\\{(x_{k},0):k\in
C_{0}\\}\right)\times\right.$
$\times\left.{\mathbf{1}}_{F}(\\{u_{1},\dots,u_{n-|C_{0}|}\\}\cup\\{(x_{k},0):k\in
C_{0}\\})\prod_{j=1}^{n-|C_{0}|}\,du_{j}\right)\prod_{i=1}^{n}\,dx_{i}.$

Plugging in the claimed expression for the conditional expectation, one
obtains

$\sum_{n=0}^{\infty}\frac{e^{-\ell({\cal{X}})}}{n!}\sum_{C_{0}\subset\\{1,...,n\\}}\left(\frac{{\mathbb{E}}Z_{1}}{{\mathbb{E}}T_{1}}\right)^{|C_{0}|}\left(\frac{{\mathbb{E}}Y_{1}}{{\mathbb{E}}T_{1}}\right)^{n-|C_{0}|}\int_{{\cal{X}}^{n}}p_{X}(\\{x_{1},...,x_{n}\\})\sum_{\begin{subarray}{c}C_{1},...,C_{n-|C_{0}|}\\\
\cup_{j}\\{C_{j}\\}=\\{1,...,n\\}\setminus
C_{0}\end{subarray}}\,\frac{1}{(n-|C_{0}|)!}$
$\int_{({\mathbb{R}}\times{\mathbb{R}}^{+})^{n-|C_{0}|}}\prod_{j=1}^{n-|C_{0}|}q_{x_{C_{j}}}(u_{j})\,c(\\{u_{1},\dots,u_{n-|C_{0}|}\\}\cup\\{(x_{k},0):k\in
C_{0}\\})$
$\left(\int_{{\cal{X}}^{n-|C_{0}|}}p_{X}(\\{y_{1},\dots,y_{n-|C_{0}|}\\}\cup\\{x_{k}:k\in
C_{0}\\})\,{\mathbf{1}}_{A}(\\{y_{1},\dots,y_{n-|C_{0}|}\\}\cup\\{x_{k}:k\in
C_{0}\\})\sum_{\begin{subarray}{c}D_{1},...,D_{n-|C_{0}|}\\\
\cup_{j}\\{D_{j}\\}=\\{1,...,n\\}\setminus C_{0}\end{subarray}}\right.$
$\left.\prod_{k=1}^{n-|C_{0}|}{\mathbf{1}}\\{y_{D_{k}}\in[u_{k,1},u_{k,1}+u_{k,2}]\\}\,dy_{k}\right){\mathbf{1}}_{F}(\\{u_{1},\dots,u_{n-|C_{0}|}\\}\cup\\{(x_{k},0):k\in
C_{0}\\})\,\prod_{j=1}^{n-|C_{0}|}\,du_{j}\prod_{i=1}^{n}dx_{i}.$

Note that in order to cancel terms, the order of integration must be changed.
As evidently the first term in $q_{x}(a,l)$, that is
${f_{Y}(l)}/{{\mathbb{E}}Y_{1}}$, does not depend on $x$, by Fubini’s theorem,

$\displaystyle{\mathbb{E}}[1_{F}(U)\,{\mathbb{P}}(X\in A\mid U)]$
$\displaystyle=\sum_{n=0}^{\infty}\frac{e^{-\ell({\cal{X}})}}{n!}\sum_{C_{0}\subset\\{1,...,n\\}}\left(\frac{{\mathbb{E}}Z_{1}}{{\mathbb{E}}T_{1}}\right)^{|C_{0}|}\left(\frac{{\mathbb{E}}Y_{1}}{{\mathbb{E}}T_{1}}\right)^{n-|C_{0}|}$
$\int_{{\cal{X}}^{n}}p_{X}(\\{y_{1},\dots,y_{n-|C_{0}|}\\}\cup\\{x_{k}:k\in
C_{0}\\})\,{\mathbf{1}}_{A}(\\{y_{1},\dots,y_{n-|C_{0}|}\\}\cup\\{x_{k}:k\in
C_{0}\\})$ $\sum_{\begin{subarray}{c}D_{1},...,D_{n-|C_{0}|}\\\
\cup_{j}\\{D_{j}\\}=\\{1,...,n\\}\setminus
C_{0}\end{subarray}}\frac{1}{(n-|C_{0}|)!}\int_{({\mathbb{R}}\times{\mathbb{R}}^{+})^{n-|C_{0}|}}{\mathbf{1}}_{F}(\\{u_{1},\dots,u_{n-|C_{0}|}\\}\cup\\{(x_{k},0):k\in
C_{0}\\})\prod_{j=1}^{n-|C_{0}|}q_{y_{D_{j}}}(u_{j})\,$
$\left(\int_{{\cal{X}}^{n-|C_{0}|}}p_{X}(\\{x_{1},\dots,x_{n}\\})\sum_{\begin{subarray}{c}C_{1},...,C_{n-|C_{0}|}\\\
\cup_{j}\\{C_{j}\\}=\\{1,...,n\\}\setminus
C_{0}\end{subarray}}\prod_{i=1}^{n-|C_{0}|}{\mathbf{1}}\\{x_{C_{i}}\in[u_{i,1},u_{i,1}+u_{i,2}]\\}\,dx_{C_{i}}\right)\times$
$\times c(\\{u_{1},\dots,u_{n-|C_{0}|}\\}\cup\\{(x_{k},0):k\in
C_{0}\\})\prod_{j=1}^{n-|C_{0}|}\,du_{j}\prod_{k=1}^{n-|C_{0}|}\,dy_{k}\prod_{i\in
C_{0}}dx_{i}={\mathbb{E}}[{\mathbf{1}}_{F}(U)\,{\mathbf{1}}_{A}(X)]$

after cancelling and rearranging terms and noting that the term in between
brackets cancels out against the normalisation constant
$c(\\{u_{1},\dots,u_{n-|C_{0}|}\\}\cup\\{(x_{k},0):k\in C_{0}\\})$. ∎

Theorem 3.1 states that the posterior distribution of $X$ given
$U={\mathbf{u}}$ is the union of $m$ atoms combined with $n-m$ points that are
distributed on ${\cal{X}}^{n-m}$ according to the symmetric probability
density function

$c({\mathbf{u}})\,p_{X}(\\{a_{1},\dots,a_{m},x_{1},\dots,x_{n-m}\\})\sum_{\begin{subarray}{c}D_{1},\dots,D_{n-m}\\\
\cup_{i}\\{D_{i}\\}=\\{1,\dots,n-m\\}\end{subarray}}\prod_{i=1}^{n-m}{\mathbf{1}}\\{x_{D_{i}}\in[a_{m+i},a_{m+i}+l_{m+i}]\\}$
(3.4)

with respect to Lebesgue measure.

######  Corollary 3.1.1.

In the framework of Theorem 3.1, the conditional distribution of the mark
assignments $D_{1},\dots,D_{n-m}$ for non-atomic marks is as follows. For
$d_{1},\dots d_{n-m}\in\\{1,\dots,n-m\\}$ such that
$\\{d_{1},\dots,d_{n-m}\\}$ $=\\{1,\dots,n-m\\}$,

${\mathbb{P}}(D_{1}=d_{1},\dots,D_{n-m}=d_{n-m}|X=\\{a_{1},\dots,a_{m},x_{1},\dots,x_{n-m}\\},U={\mathbf{u}})$
$=\frac{\prod_{i=1}^{n-m}{\mathbf{1}}\\{x_{d_{i}}\in[a_{m+i},a_{m+i}+l_{m+i}]\\}}{\sum_{\begin{subarray}{c}C_{1},\dots,C_{n-m}\\\
\cup_{i}\\{C_{i}\\}=\\{1,\dots,n-m\\}\end{subarray}}\prod_{i=1}^{n-m}{\mathbf{1}}\\{x_{C_{i}}\in[a_{m+i},a_{m+i}+l_{m+i}]\\}}$

provided that $x_{i}\in[a_{m+i},a_{m+i}+l_{m+i}]$ for $i=1,\dots,n-m$ and zero
otherwise.

As a special case, let us consider an inhomogeneous Poisson process with
integrable intensity function $\lambda:{\cal{X}}\to{\mathbb{R}}^{+}$. Then,
under the posterior distribution, $X$ consists of $n$ independent points, one
in each interval of ${\mathbf{u}}$, with probability density function

$\frac{\lambda(x)}{\int_{[a_{i},a_{i}+l_{i}]\cap{\cal{X}}}\lambda(s)ds}$

on $[a_{i},a_{i}+l_{i}]\cap{\cal{X}}$ for intervals with $l_{i}>0$. To see
this, recall that for a Poisson process (Van Lieshout, 2000)

$p_{X}(\\{a_{1},\dots,a_{m},x_{1},\dots,x_{n-m}\\})=\exp\left[\int_{{\cal{X}}}(1-\lambda(s))ds\right]\prod_{j=1}^{m}\lambda(a_{j})\prod_{i=1}^{n-m}\lambda(x_{i})$

which factorises over terms associated with each interval. Hence (3.4) is
proportional to

$\sum_{\begin{subarray}{c}D_{1},\dots,D_{n-m}\\\
\cup_{i}\\{D_{i}\\}=\\{1,\dots,n-m\\}\end{subarray}}\prod_{i=1}^{n-m}\lambda(x_{D_{i}})\,{\mathbf{1}}\\{x_{D_{i}}\in[a_{m+i},a_{m+i}+l_{m+i}]\\}.$

##  4 Statistical inference

In this section we will consider statistical inference for aoristically
censored data. Our main aim is to reconstruct the latent point process $X$
from observed parametrised intervals $U$, that may or may not be censored. In
tandem, the censoring probability as well as the parameters $\eta$ of the
distribution of the non-degenerate intervals must be estimated. Parameters of
the prior distribution may either be treated as fixed or subject to
estimation.

###  4.1 Forward model parameters

Suppose that we observe a realisation
${\mathbf{u}}=\\{(a_{1},0),\dots,(a_{m},0),(a_{m+1},l_{m+1}),\dots,(a_{n},l_{n})\\}$
of $U$, where $a_{i}\in{\mathbb{R}}$, $l_{i}>0$ and $n\neq 0$. Our first aim
is to estimate the parameters $\eta$ of the mark distribution $\nu$ (cf.
Theorem 2.2). The parameter vector $\eta$ comprises the parameters $\zeta$ of
the probability density function $f_{Y}$ as well as any other parameters
$\theta$ involved in the joint distribution of the random vector
$C_{1}=(Y_{1},Z_{1})$ that defines the alternating renewal process (cf.
Section 2.1).

The likelihood function can be obtained from the proof of Theorem 3.1 by
taking $A$ equal to ${\cal{N}}_{\cal{X}}$ in equation (3.3). On a logarithmic
scale,

$L(\eta;{\mathbf{u}})=m\log\left(\frac{{\mathbb{E}}[Z_{1};\zeta,\theta]}{{\mathbb{E}}[T_{1};\zeta,\theta]}\right)+(n-m)\log\left(\frac{{\mathbb{E}}[Y_{1};\zeta]}{{\mathbb{E}}T_{1}[\zeta,\theta]}\right)+\sum_{i=1}^{n-m}\log\left(\frac{f_{Y}(l_{i};\zeta)}{{\mathbb{E}}[Y_{1};\zeta]}\right)$
(4.1)

upon ignoring terms that do not depend on $\eta$.

Equation (4.1) simplifies greatly if we assume that the mixture weight
$p={\mathbb{E}}[Z_{1};\zeta,\theta]/{\mathbb{E}}[T_{1};\zeta,\theta]$ does not
depend on $\zeta$. Then $\eta=(p,\zeta)$ and

$L(p,\zeta;{\mathbf{u}})=m\log
p+(n-m)\log(1-p)+\sum_{i=1}^{n-m}\log\left(\frac{f_{Y}(l_{i};\zeta)}{{\mathbb{E}}[Y_{1};\zeta]}\right).$

The atom probability $p$ may be estimated by $m/n$, the fraction of atoms in
the sample ${\mathbf{u}}$. For $\zeta$, we need the following result.

######  Proposition 4.1.

Let $\nu$ be as in Theorem 2.2. Then the distribution of the lengths of non-
degenerate intervals is given by the length-weighted marginal distribution
$f(l)={lf_{Y}(l)}/{{\mathbb{E}}Y_{1}}$ and the left-most points are,
conditionally on $L=l$, uniformly distributed on $[-l,0]$.

###### Proof..

Let $f(l)$ be the marginal distribution of the length $l$. Evaluating,

$\displaystyle f(l)$
$\displaystyle=\int\frac{f_{Y}(l)}{{\mathbb{E}}Y_{1}}\,{\mathbf{1}}\\{a\leq
0\leq
a+l\\}\,da=\int_{-l}^{0}\frac{f_{Y}(l)}{{\mathbb{E}}Y_{1}}\,da=\frac{lf_{Y}(l)}{{\mathbb{E}}Y_{1}}.$

Let $f_{A|L=l}(a)$ be the conditional probability density function for the
left-most point $A$ of an interval given its length $L$. Using the definition
of conditional density,

$\displaystyle
f_{A\,|L=l}(a)=\frac{f(a,l)}{f(l)}=\frac{\frac{f_{Y}(l)}{{\mathbb{E}}Y_{1}}\,{\mathbf{1}}\\{a\leq
0\leq
a+l\\}}{\frac{lf_{Y}(l)}{{\mathbb{E}}Y_{1}}}=\frac{1}{l}\,{\mathbf{1}}\\{a\in[-l,0]\\}.$

Thus $A\sim\text{Unif}[-l,0]$. ∎

When the censoring probability does not depend on $\zeta$, the latter may be
estimated by treating the non-degenerate intervals as an independent sample
from $f(l)$ and applying the maximum likelihood method. For example, if
$f_{Y}$ is the probability density function of a $\text{Gamma}(k,\lambda)$
distribution with shape parameter $k>0$ and rate parameter $\lambda>0$, $f(l)$
is the probability density of a Gamma distribution with parameters $k+1$ and
$\lambda$.

###  4.2 State estimation

Since the posterior distribution of $X$ or $W$ given $U$ (cf. Theorem 3.1) is
intractable because of the normalisation constant $c({\mathbf{u}})$, we will
use Markov chain Monte Carlo methods (Brooks, Gelman, Jones & Meng, 2011;
Møller & Waagepetersen, 2003) for simulation. These methods construct a Markov
chain in such a way that the stationary distribution of the chain is exactly
the posterior distribution. Of these methods, a Metropolis-Hastings algorithm
with a fixed number of points will be used. Since the transition probabilities
depend on likelihood ratios, the benefit is that one can sample from
unnormalised densities.

Let us return to the framework of Theorem 3.1. Note that sampling from the
posterior distribution of $X$ given $U$ is cumbersome due to the presence of
the permutation sum term in (3.4). Therefore our approach is to sample from
the posterior distribution of the complete model $W$ and project on its ground
process of locations. Doing so avoids attributing points to intervals and
therefore avoids the intractable sum. Moreover, as we saw in Section 2.3, $W$
has probability density function $p_{X}$ with respect to a unit rate Poisson
process on ${\cal{X}}\times({\mathbb{R}}\times{\mathbb{R}}^{+})$ with
intensity measure $\ell\times\nu$. Upon observing $U={\mathbf{u}}$ for
${\mathbf{u}}=\\{(a_{1},0),\dots,(a_{m},0),(a_{m+1},l_{m+1}),\dots,(a_{n},l_{n})\\}$,
by (3.4) this means that we must sample from a probability density function
$\pi$ on ${\cal{X}}^{n-m}$ that is proportional to
$p_{X}(\\{a_{1},\dots,a_{m},x_{1},\dots,x_{n-m}\\})$. The ordering of the
points inherent in working on ${\cal{X}}^{n-m}$ represents the unique
correspondence between points in $X$ and intervals in $U$ in the complete
model. We will use the notation $\bar{\mathbf{x}}$ to indicate that we look at
vectors rather than sets ${\mathbf{x}}$. In the special case that $n=m$, all
points are observed perfectly and there is no need for any simulation.

We will use the Metropolis–Hastings algorithm (Brooks, Gelman, Jones & Meng,
2011) when $n>m$, i.e. when there are density-admitting points. The state
space is given by

$\displaystyle\overline{E}({\mathbf{u}})=\\{(x_{1},\dots,x_{n-m})\in{\cal{X}}^{n-m}:x_{i}\in{\cal{X}}\cap[a_{m+i},a_{m+i}+l_{m+i}],\,p_{X}(\\{a_{1},\dots,a_{m},x_{1},\dots,x_{n-m}\\})>0\\}.$

From now on, we shall assume that the state space is non-degenerate in the
sense that

$\int_{\overline{E}({\mathbf{u}})}p_{X}(\\{a_{1},\dots,a_{m},x_{1},\dots,x_{n-m}\\})\,dx_{1}\dots
dx_{n-m}>0.$ (4.2)

Now, the Metropolis–Hastings algorithm is defined as follows. Let
$q:\overline{E}({\mathbf{u}})\times\overline{E}({\mathbf{u}})\to{\mathbb{R}}^{+}$
be a Markov kernel. Iteratively, if the current state is
$\bar{\mathbf{x}}\in\overline{E}({\mathbf{u}})$, propose a new state
$\bar{\mathbf{y}}\in\overline{E}({\mathbf{u}})$ according to the probability
density function $q(\bar{\mathbf{x}},\cdot)$ and accept the proposal to move
to $\bar{\mathbf{y}}$ with probability

$\alpha(\bar{\mathbf{x}},\bar{\mathbf{y}})=\left\\{\begin{array}[]{ll}1&\mbox{
if
}p_{X}(\\{a_{1},\dots,a_{m},y_{1},\dots,y_{n-m}\\})\,q(\bar{\mathbf{y}},\bar{\mathbf{x}})\\\
&\quad\quad\geq
p_{X}(\\{a_{1},\dots,a_{m},x_{1},\dots,x_{n-m}\\})\,q(\bar{\mathbf{x}},\bar{\mathbf{y}});\\\
\frac{p_{X}(\\{a_{1},\dots,a_{m},y_{1},\dots,y_{n-m}\\})\,q(\bar{\mathbf{y}},\bar{\mathbf{x}})}{p_{X}(\\{a_{1},\dots,a_{m},x_{1},\dots,x_{n-m}\\})\,q(\bar{\mathbf{x}},\bar{\mathbf{y}})\\}}&\mbox{
otherwise.}\end{array}\right.$ (4.3)

When the proposal is rejected, stay in the current state $\bar{\mathbf{x}}$.
The choice of $q$ depends on $p_{X}$. In our simulations in Section 5, we will
use the following algorithm which is valid when the prior density function
$p_{X}$ is strictly positive.

######  Algorithm 4.2.

Supppose that $p_{X}>0$ and $n>m$. Iteratively, if the current state is
$\bar{\mathbf{x}}\in\overline{E}({\mathbf{u}})$,

  * •

pick an interval $[a_{m+i},a_{m+i}+l_{m+i}]$, $i=1,\dots,n-m$, uniformly at
random from the non-degenerate ones;

  * •

generate a uniformly randomly distributed point $y_{i}$ on
${\cal{X}}\cap[a_{m+i},a_{m+i}+l_{m+i}]$ and propose to update $x_{i}$ to
$y_{i}$;

  * •

accept the proposal with probability

$\alpha_{i}((x_{1},\dots,x_{n-m}),y_{i})=\min\left(1,\frac{p_{X}((\\{a_{1},\dots,a_{m},x_{1},\dots,x_{n-m}\\}\setminus\\{x_{i}\\})\cup\\{y_{i}\\})}{p_{X}(\\{a_{1},\dots,a_{m},x_{1},\dots,x_{n-m}\\})}\right)$
(4.4)

and otherwise stay in the current state.

A few remarks are in order. First, note that since ${\cal{X}}$ is open, the
intersection with closed intervals that contain a point in ${\cal{X}}$ is also
non-degenerate when $l_{i}>0$. Secondly, when $p_{X}$ may take the value zero,
the proposal mechanism in Algorithm 4.2 might result in a new state that does
not belong to $\overline{E}({\mathbf{u}})$, even when $\bar{\mathbf{x}}$ does.
Moreover, only changing one component at a time might lead to non-irreducible
Markov chains. For example, if ${\mathbf{u}}$ contains the parametrisations of
the intervals $[0,1]$ and $[0.1,1]$ and $p_{X}({\mathbf{x}})=0$ for
realisations ${\mathbf{x}}$ that contain components separated by a distance
less than $0.55$, then states such as $\bar{\mathbf{x}}=(0.3,0.9)$ and
$\bar{\mathbf{y}}=(0.9,0.3)$ cannot be reached from one another.

Let the target distribution be $\pi$ given by (3.4) interpreted as a
probability density on $\overline{E}({\mathbf{u}})$. In the next propositions,
basic properties of the algorithm are considered. The proofs are modifications
to our context of classic Metropolis-Hastings proofs found in, for example,
(Mengersen & Tweedie, 1996), (Roberts & Smith, 1994) or (Møller &
Waagepetersen, 2003, Chapter 7).

We will write $Y_{i}$ for subsequent states and denote by
$P(\bar{\mathbf{x}},F)=P(Y_{i+1}\in F\,|\,Y_{i}=\bar{\mathbf{x}})$ the
transition probability from state
$\bar{\mathbf{x}}\in\overline{E}({\mathbf{u}})$ into
$F\subset\overline{E}({\mathbf{u}})$.

######  Proposition 4.3.

Consider the set-up of Theorem 3.1 with $n>m$ and assume that condition (4.2)
is met. Then, the Metropolis-Hastings algorithm defined by Markov kernel $q$
on $\overline{E}({\mathbf{u}})$ and acceptance probabilities (4.3) is
reversible with respect to $\pi$.

###### Proof..

Take $\bar{\mathbf{x}},\bar{\mathbf{y}}$ in $\overline{E}({\mathbf{u}})$ and
assume that $\pi(\bar{\mathbf{y}})\,q(\bar{\mathbf{y}},\bar{\mathbf{x}})$
$>\pi({\mathbf{x}})\,q(\bar{\mathbf{x}},\bar{\mathbf{y}})\geq 0$. Then

$\displaystyle\pi(\bar{\mathbf{x}})\,q(\bar{\mathbf{x}},\bar{\mathbf{y}})\,\alpha(\bar{\mathbf{x}},\bar{\mathbf{y}})=c({\mathbf{u}})\,p_{X}(\\{a_{1},\dots,a_{m},x_{1},\dots,x_{n-m}\\})\,q(\bar{\mathbf{x}},\bar{\mathbf{y}})=$
$\displaystyle
c({\mathbf{u}})\,p_{X}(\\{a_{1},\dots,a_{m},y_{1},\dots,y_{n-m}\\})\,q(\bar{\mathbf{y}},\bar{\mathbf{x}})\,\frac{p_{X}(\\{a_{1},\dots,a_{m},x_{1},\dots,x_{n-m}\\})\,q(\bar{\mathbf{x}},\bar{\mathbf{y}})}{p_{X}(\\{a_{1},\dots,a_{m},y_{1},\dots,y_{n-m}\\})\,q(\bar{\mathbf{y}},\bar{\mathbf{x}})}=\pi({\mathbf{y}})\,q(\bar{\mathbf{y}},{\mathbf{x}})\,\alpha(\bar{\mathbf{y}},{\mathbf{x}})$

writing $c({\mathbf{u}})$ for the normalisation constant. We conclude that the
chain is in detailed balance and therefore reversible with respect to $\pi$. ∎

Recall that the Markov chain is called $\pi$-irreducible (Meyn & Tweedie,
2009) if for every $\bar{\mathbf{x}}\in\overline{E}({\mathbf{u}})$ and every
$F\subset\overline{E}({\mathbf{u}})$ with $\pi(F)>0$ there exists some natural
number $\tau$ such that $P^{\tau}(\bar{\mathbf{x}},F)>0$.

######  Proposition 4.4.

Consider the set-up of Theorem 3.1 with $n>m$ and assume that condition (4.2)
is met. Let $Q$ be the one-step transition kernel of the Markov chain on
$\overline{E}({\mathbf{u}})$ generated by Markov kernel
$q:\overline{E}({\mathbf{u}})\times\overline{E}({\mathbf{u}})\to{\mathbb{R}}^{+}$
in which every proposal is accepted. If the chain defined by $Q$ is
$\pi$-irreducible and $q(\bar{\mathbf{x}},\bar{\mathbf{y}})=0$ if and only if
$q(\bar{\mathbf{y}},\bar{\mathbf{x}})=0$, then the Metropolis-Hastings
algorithm defined by $q$ and (4.3) is $\pi$-irreducible. In particular, the
chain of Algorithm 4.2 is $\pi$-irreducible when $p_{X}>0$.

###### Proof..

The first part follows from (Roberts & Smith, 1994, Theorem 3.ii).

For Algorithm 4.2, $q(\bar{\mathbf{x}},\bar{\mathbf{y}})>0$ only if
$\bar{\mathbf{x}}$ and $\bar{\mathbf{y}}$ in $\overline{E}({\mathbf{u}})$
differ in at most a single component. Thus, assume that $x_{j}=y_{j}$ for all
$j\neq i\in\\{1,\dots,n-m\\}$ and $x_{i}\neq y_{i}$. Then
$q(\bar{\mathbf{x}},\bar{\mathbf{y}})=q(\bar{\mathbf{y}},\bar{\mathbf{x}})$ so
they are strictly positive or zero together. Write $Q^{\tau}$ for the
$\tau$-step transition kernel of the aways-accept chain. Then, for
$\bar{\mathbf{x}}$, $\bar{\mathbf{y}}\in\overline{E}({\mathbf{u}})$,

$q^{n-m}(\bar{\mathbf{x}},\bar{\mathbf{y}})\geq\left(\frac{1}{n-m}\right)^{n-m}\,\prod_{i=1}^{n-m}\frac{1}{\ell({\cal{X}}\cap[a_{m+i},a_{m+i}+l_{m+i}])}>0$

by changing each component in turn. We conclude that the Markov chain of
Algorithm 4.2 is $\pi$-irreducible. ∎

Recall that a $\pi$-irreducible Markov chain is called aperiodic (Meyn &
Tweedie, 2009) if the state space ($\overline{E}({\mathbf{u}})$ in our case)
cannot be partitioned into measurable sets $B_{0}$, $B_{1},\dots B_{d-1}$ such
that $\pi(\overline{E}({\mathbf{u}})\setminus\cup_{j=0}^{d-1}B_{j})=0$ and
$P(\bar{\mathbf{x}},B_{j+1\,\text{mod}\,d})$ $=1$ for all $\bar{\mathbf{x}}\in
B_{j}$ (for some $d>1$, the period). By (Møller & Waagepetersen, 2003,
Proposition 7.6), a $\pi$-irreducible Markov chain with invariant probability
distribution $\pi$ is aperiodic if and only if for some small set $D$ with
$\pi(D)>0$ and some $\tau\in{\mathbb{N}}$, the following holds:
$P^{i}(\bar{\mathbf{x}},D)>0$ for all $\bar{\mathbf{x}}\in D$ and $i\geq\tau$.

######  Proposition 4.5.

Consider the set-up of Theorem 3.1 with $n>m$ and assume that condition (4.2)
is met. If $0<p_{X}(\\{a_{1},\dots,$ $a_{m},$ $x_{1},\dots,x_{n-m})\leq\delta$
for some $\delta>0$ and all $\bar{\mathbf{x}}\in\overline{E}({\mathbf{u}})$,
then the Metropolis-Hastings Markov chain of Algorithm 4.2 is aperiodic.

###### Proof..

Let $\xi$ be the point on ${\cal{X}}\cap[a_{m+1},a_{m+1}+l_{m+1}]$ that
replaces $x_{1}$. By (4.2), there exist $x_{2},\dots,x_{n-m}$ such that

$\int_{{\cal{X}}\cap[a_{m+1},a_{m+1}+l_{m+1}]}c({\mathbf{u}})\,p_{X}(\\{a_{1},\dots,a_{m},\xi,x_{2},\dots,x_{n-m})\,d\xi=\int_{{\cal{X}}\cap[a_{m+1},a_{m+1}+l_{m+1}]}\pi(\xi,x_{2},\dots,x_{n-m})\,d\xi$

is strictly positive. Define a measure $\mu$ on the Borel $\sigma$-algebra on
$\overline{E}({\mathbf{u}})$ by

$\mu(F)=\int_{{\cal{X}}\cap[a_{m+1},a_{m+1}+l_{m+1}]}{\mathbf{1}}\\{(\xi,x_{2},\dots,x_{n-m})\in
F\\}\,\pi(\xi,x_{2},\dots,x_{n-m})\,d\xi$

and note that $\mu(\overline{E}({\mathbf{u}}))>0$. Set
$C=\\{(\xi,x_{2},\dots,x_{n-m}):\xi\in{\cal{X}}\cap[a_{m+1},a_{m+1}+l_{m+1}]\\}$.
We claim that $C$ is small with respect to $\mu$. To see this, take
$\bar{\mathbf{y}}=(y,x_{2},\dots,x_{n-m})\in C$ and note that for
$F\subset\overline{E}({\mathbf{u}})$, the transition probability
$P(\bar{\mathbf{y}},F)$ is at least

$\frac{1}{n-m}\frac{1}{\ell({\cal{X}}\cap[a_{m+1},a_{m+1}+l_{m+1}])}\int_{{\cal{X}}\cap[a_{m+1},a_{m+1}+l_{m+1}]}{\mathbf{1}}_{F}(\xi,x_{2},\dots,x_{n-m})\,\alpha_{1}((y,x_{2},\dots,x_{n-m}),\xi)\,d\xi.$

If $\pi(\xi,x_{2},\dots,x_{n-m})>\pi(y,x_{2},\dots,x_{n-m})$, then
$\alpha_{1}((y,x_{2},\dots,x_{n-m}),\xi)=1>\pi(\xi,x_{2},\dots,x_{n-m})/(c({\mathbf{u}})\,\delta)$,
again writing $c({\mathbf{u}})$ for the normalisation constant. Otherwise,
$\alpha_{1}((x_{1},\dots,x_{n-m}),\xi)=\pi(\xi,x_{2},\dots,x_{n-m})/$
$\pi(x_{1},x_{2},\dots,x_{n-m})>\pi(\xi,x_{2},\dots,x_{n-m})/(c({\mathbf{u}})\,\delta)$.
In summary,

$P(\bar{\mathbf{x}},F)>\frac{1}{n-m}\frac{1}{\ell({\cal{X}}\cap[a_{m+1},a_{m+1}+l_{m+1}])}\,\frac{\mu(F)}{c({\mathbf{u}})\,\delta}$

so $C$ is small with respect to $\mu$. Moreover,
$\pi(C)=\mu(\overline{E}({\mathbf{u}}))>0$. Iterating the above argument one
notices that $P^{\tau}(\bar{\mathbf{x}},F)$ is at least as large as the
$\tau$-th power of the bound above, an observation that completes the proof. ∎

In conclusion, from almost all initial states, Algorithm 4.2 converges in
total variation to the invariant probability distribution. Conditions for
general proposal kernels $q$ can be found in (Møller & Waagepetersen, 2003,
Chapter 7) or (Roberts & Smith, 1994, Theorem 3).

###  4.3 Prior model parameters

In Section 4.2, we discussed Monte Carlo methods to sample from the posterior
distribution of $W$ or $X$ given $U={\mathbf{u}}$. This distribution is
defined in terms of the prior probability density function $p_{X}$. Typically,
$p_{X}$ is given in unnormalised form and depends on a parameter vector
$\theta$, that is,
$p_{X}({\mathbf{x}};\theta)=c(\theta)h_{X}({\mathbf{x}};\theta)$ for an
explicit function $h_{X}:{\cal{N}}_{\cal{X}}\to{\mathbb{R}}^{+}$. When
$\theta$ is treated as an unknown, since the likelihood function for $\theta$
contains the latent marked point process $W$, we call on techniques from
missing data analysis.

The likelihood function $l(\theta)$ is obtained from the proof of Theorem 3.1
by taking $A$ equal to ${\cal{N}}_{\cal{X}}$ in equation (3.3). Disregarding
terms that do not depend on $\theta$, one obtains
$l(\theta;{\mathbf{u}})=c(\theta)c(\theta|{\mathbf{u}})^{-1}$ where

$c(\theta)^{-1}=\sum_{k=0}^{\infty}\frac{e^{-\ell({\cal{X}})}}{k!}\int_{{\cal{X}}^{k}}h_{X}(\\{x_{1},\dots,x_{k}\\};\theta)\,dx_{1}\dots
dx_{k}$

and $c(\theta|{\mathbf{u}})^{-1}$ is given by

$\int_{{\cal{X}}^{n-m}}h_{X}(\\{a_{1},\dots,a_{m},x_{1},\dots,x_{n-m}\\};\theta)\left(\sum_{\begin{subarray}{c}D_{1},\dots,D_{n-m}\\\
\cup_{i}\\{D_{i}\\}=\\{1,\dots,n-m\\}\end{subarray}}{\mathbf{1}}\\{x_{D_{i}}\in[a_{m+i},a_{m+i}+l_{m+i}]\\}\right)\prod_{i=1}^{n-m}dx_{i}.$

One observes that $c(\theta|{\mathbf{u}})$ is equal to the normalisation
constant on $\overline{E}({\mathbf{u}})$ of the non-atomic part of the
posterior distribution of $X$ given $U={\mathbf{u}}$.

To handle the two normalisation constants, it is necessary to look at the log
relative likelihood $L(\theta)$ of $U$ with respect to some fixed and user-
selected reference parameter $\theta_{0}$, Then, as in (Gelfand & Carlin,
1993; Geyer, 1999),

$L(\theta)=\log\left[\frac{c(\theta)\,c(\theta_{0}|{\mathbf{u}})}{c(\theta_{0})\,c(\theta|{\mathbf{u}})}\right]=\log{\mathbb{E}}_{\theta_{0}}\left[\frac{h_{X}(X;\theta)}{h_{X}(X;\theta_{0})}\mid
U={\mathbf{u}}\right]-\log{\mathbb{E}}_{\theta_{0}}\left[\frac{h_{X}(X;\theta)}{h_{X}(X;\theta_{0})}\right].$

Being expressible in terms of expectations under the reference parameter, the
log likelihood ratio can be approximated by Markov chain Monte Carlo methods.
Note that two samples are required: one from the posterior distribution of $X$
and one from the prior. For the latter, provided $p_{X}$ is locally stable,
classic Metropolis–Hastings methods based on births and deaths apply (Geyer,
1999; Møller & Waagepetersen, 2003). If the conditional intensity is monotone,
exact simulation can be carried out (Kendall & Møller, 2000; Van Lieshout &
Baddeley, 2002).

##  5 Examples

In this section, we present a few examples to illustrate how the choice of
prior affects state estimation. Calculations were carried out using the C++
marked point process library MPPLIB, developed by Steenbeek et al. For $p_{X}$
we choose the area-interaction point process (Widom & Rowlinson, 1970;
Baddeley & Van Lieshout, 1995), a model that favours clustered, regular and
random realisations depending on parameter values. Specifically, this model
has probability density function

$\displaystyle
p(\mathbf{x})=\alpha\beta^{n(\mathbf{x})}\exp\left[-\log\gamma\,\ell({\cal{X}}\cap
U_{r}(\mathbf{x}))\right]$ (5.1)

with respect to a unit rate Poisson process on $\mathcal{X}$. Here
$U_{r}(\mathbf{x})=\bigcup_{i=1}^{n}B(x_{i},r)$ where $B(x_{i},r)$ is the
closed interval $[x_{i}-r,x_{i}+r]$. When $\gamma<1$, realisations tend to be
regular, for $\gamma>1$ clustered. When $\gamma=1$, one has a Poisson process
with intensity $\beta$. The scalar $\alpha$ is a normalisation constant.
Realisations can be obtained by Kendall’s dominated coupling from the past
(CFTP) algorithm (Kendall, 1998) developed initially from the perfect
simulation methods of Propp and Wilson for coupled Markov chains (Propp &
Wilson, 1996).

###  5.1 Toy example

Consider data ${\mathbf{u}}=\\{(0.45,0.4),(0.51,0),(0.58,0)\\}$ that consist
of two atoms and a single non-degenerate interval. By the discussion at the
end of Section 3, for a Poisson prior ($\gamma=1$), the posterior distribution
of the location $X_{3}$ in ${\cal{X}}=(0,1)$ that generated the non-degenerate
interval is uniformly distributed. To see the effect of informative priors,
Figure 2 plots the posterior distribution of $X_{3}$ when the prior is an
area-interaction model with $\eta=2r\log\gamma=1.2$ and $r=0.1$. Note that
mass is shifted to the left side of the interval due to the presence of atoms.
For $\eta=-1.2$ and $r=0.1$, the atoms repel $X_{3}$, resulting in mass being
shifted to the right side of the interval (cf. Figure 3). To carry out the
state estimation, we ran Algorithm 4.2 with a burn-in of 10,000 steps and
calculated the histograms based on the subsequent 100,000 steps.

Figure 2: Locations of two atoms and a spanning interval together with a
histogram of point locations. The interval start and end points as well as the
atom locations are marked on the histogram $x$-axis.

Figure 3: Locations of two atoms and a spanning interval together with a
histogram of point locations. The interval start and end points as well as the
atom locations are marked on the histogram $x$-axis.

###  5.2 Area-interaction gamma model

The left-most panel in Figure 4 shows a simulation in ${\cal{X}}=(0,1)$ from
$U$ in a model where $X$ is an area-interaction process with parameters
$\beta=12$, $\eta=2r\log\gamma=0$ and $r=0.05$ marked by a mixture
distribution $\nu$ in which the atom probability is $p=0.2$ and $f_{Y}$ is the
probability density function of a Gamma distribution with shape parameter
$k=2.5$ and rate parameter $\lambda=0.07$. The points shown as black dots are
the points of $X$ in the simulated pattern, the red points constitute a
realisation from the posterior distribution of $X$ given $U$ obtained by
running Algorithm 4.2 for 10,000 steps. The points seem to settle in a random
manner within the intervals.

Figure 4: Plots of the simulated and new locations of a random area-
interaction point process. Parameter values:
$(\beta,\eta,r,\lambda,k,p)=(12,0,0.05,0.07,2.5,0.2)$.

A simulation using a prior favouring clustering can be found in black in
Figure 5. The parameter settings were as before except that $\eta=1.2$. In
red, a realisation from the posterior distribution is shown, obtained after
10,000 steps from the Metropolis–Hastings algorithm. Figure 5 shows the effect
of the complex underlying geometry when choosing proposal points within
parametrised intervals. The algorithm tends to move proposed times to areas
where multiple intervals intersect, leading to clustering within these
regions.

Figure 5: Parameter values:
$(\beta,\eta,r,\lambda,k,p)=(12,1.2,0.05,0.07,2.5,0.2)$.

Figure 6 shows, in black, a simulation from $U$ for a regular area-interaction
prior with $\eta=2r\log\gamma=-1.2$. A realisation from the posterior
distribution of $X$ is shown in red. The structure of the prior point process
is maintained in the posterior, with points being spread out from each other.

Figure 6: Parameter values:
$(\beta,\eta,r,\lambda,k,p)=(12,-1.2,0.05,0.07,2.5,0.2)$.

##  6 Conclusion

In this work, a Bayesian inference framework for aoristic data was introduced
in which an alternating renewal process is used to interval censor temporal
data, converting it into a marked point process model. A prospective point,
which cannot be observed directly, was paired with an interval within which
the point surely lies. State estimation was then applied to best estimate the
location of this point. Theory was developed regarding the distribution of
these marks based on this renewal framework and the posterior distribution
deduced. The fact that the forward model allows for a mixture of discrete and
absolutely continuous components makes this process nontrivial. A state
estimation procedure was outlined in the form of a Metropolis-Hastings
algorithm for a fixed number of points, after which ergodicity properties were
verified. Using an area-interaction prior, this procedure was applied to
sample from the posterior distribution. Effects of the prior are clearly
present when sampling from the complete model.

Throughout, we assumed that all intervals corresponding to a point in
${\cal{X}}$ were observed. Returning to a criminology context, sampling bias
may arise since the data may contain only intervals whose right-most point is
in a given interval. Additionally, a random labelling regime was assumed. It
might be more realistic to have location-dependent independent marking, for
example based on a semi-Markov process rather than an alternating renewal
process. Furthermore, spatial aspects were completely ignored. These
generalisations will form the topic for our future research.

## Acknowledgements

This research was funded by NWO, the Dutch Research Council (grant
OCENW.KLEIN.068).

## References

  * Asmussen (2003) Asmussen, S. (2003). Applied probability and queues. Springer.
  * Baddeley & Van Lieshout (1995) Baddeley, A. & Lieshout, M.N.M. van. (1995). Area-interaction point processes. Annals of the Institute of Statistical Mathematics, 47, 601–619.
  * Bernasco (2009) Bernasco, W. (2009). Burglary. In Tonry, M. (Ed.), Oxford handbook of crime and public policy. Oxford University Press (pp. 165–190).
  * Brémaud (1972) Brémaud, P. (1972). A martingale approach to point processes (PhD thesis). University of Berkeley.
  * Brix & Diggle (2001) Brix, A. & Diggle, P.J. (2001). Spatio-temporal prediction for log-Gaussian Cox processes. Journal of the Royal Statistical Society: Series B, 63, 823–841.
  * Brooks, Gelman, Jones & Meng (2011) Brooks, S., Gelman, A., Jones, G. & Meng, X.L. (2011). Handbook of Markov chain Monte Carlo. CRC Press.
  * Daley & Vere-Jones (2003, 2008) Daley, D.J. & Vere–Jones, D. (2003, 2008). An introduction to the theory of point processes (2nd ed., in two volumes). Springer.
  * Dempster, Laird & Rubin (1997) Dempster, A.P., Laird, N.M. & Rubin, D.B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B, 39, 1–38.
  * Gamerman & Lopes (2006) Gamerman, D. & Lopes, H.F. (2006). Markov chain Monte Carlo: stochastic simulation for Bayesian inference. CRC Press.
  * Gelfand & Carlin (1993) Gelfand, A.E. & Carlin, B.P. (1993). Maximum-likelihood estimation for constrained- or missing-data models. Canadian Journal of Statistics, 21, 303–311.
  * Geyer (1999) Geyer, C.J. (1999). Likelihood inference for spatial point processes. In O. Barndorff–Nielsen, W.S. Kendall and M.N.M. van Lieshout (Eds.), Stochastic geometry, likelihood, and computation. CRC Press (pp. 141–172).
  * Helms (2008) Helms, D. (2008). Temporal analysis. In Gwinn, S.L., Bruce, C., Cooper, J.P. & Hick, S. (Eds.), Exploring crime analysis: readings on essential skills (2nd ed.). International Association of Crime Analysts (pp. 214–257).
  * Kendall (1998) Kendall, W. (1998). Perfect simulation for the area-interaction point process. In Accardi, L. & Heyde, C. (Eds.), Probability towards the year 2000. Springer (pp. 218–234).
  * Kendall & Møller (2000) W.S. Kendall & J. Møller (2000). Perfect simulation using dominating processes on ordered spaces, with application to locally stable point processes. Advances in Applied Probability, 32, 844–865.
  * Karr (1991) Karr, A.F. (1991). Point processes and their statistical inference (2nd ed.). Marcel Dekker.
  * Last & Brandt (1995) Last, G. & Brandt, A. (1995). Marked point processes on the real line. The dynamic approach. Springer.
  * Van Lieshout (1995) Lieshout, M.N.M. van (1995). Stochastic geometry models in image analysis and spatial statistics. CWI Tracts.
  * Van Lieshout (2000) Lieshout, M.N.M. van (2000). Markov point processes and their applications. Imperial College Press.
  * Van Lieshout (2016) Lieshout, M.N.M. van (2016). Likelihood based inference for partially observed renewal processes. Statistics and Probability Letters, 118, 190–196.
  * Van Lieshout & Baddeley (2002) Lieshout, M.N.M. van & Baddeley, A.J. (2002). Extrapolating and interpolating spatial patterns. In Lawson, A.B. and Denison, D.G.T. (Eds.), Spatial cluster modelling. CRC Press (pp. 61–86).
  * Mengersen & Tweedie (1996) Mengersen, K.L. & Tweedie, R.L. (1996). Rates of convergence of the Hastings and Metropolis algorithms. Annals of Statistics, 24, 101–121.
  * Meyn & Tweedie (2009) Meyn, S. & Tweedie, R.L. (2009). Markov chains and stochastic stability (2nd ed.). Cambridge University Press.
  * Mulder, Ruiter & Klugkist (2018) Mulder, K., Ruiter, S. & Klugkist, I. (2018). Analysis of circular interval-censored data motivated by aoristic data in criminology. Proceedings CMStatistics.
  * Møller & Waagepetersen (2003) Møller, J. & Waagepetersen, R.P. (2003). Statistical inference and simulation for spatial point processes. CRC Press.
  * Propp & Wilson (1996) Propp, J.G. & Wilson, D.B. (1996). Exact sampling with coupled Markov chains and applications to statistical mechanics. Random Structures and Algorithms, 9, 223–252.
  * Ratcliffe & McCullagh (1998) Ratcliffe, J.H. & McCullagh, M.J. (1998). Aoristic crime analysis. International Journal of Geographical Information Science, 12, 751–764.
  * Ripley & Kelly (1977) Ripley, B.D. & Kelly, F.P. (1977). Markov point processes. Journal of the London Mathematical Society, 15, 188–192.
  * Roberts & Smith (1994) Roberts, G.O. & Smith, A.F.M. (1994). Simple conditions for the convergence of the Gibbs sampler and Metropolis–Hastings algorithms. Stochastic Processes and their Applications, 49, 207–216.
  * Ross (1996) Ross, S.M. (1996). Stochastic processes (2nd ed.). Wiley.
  * Widom & Rowlinson (1970) Widom, B. & Rowlinson, J.S. (1970). A new model for the study of liquid-vapor phase transitions. Journal of Chemical Physics, 52, 1670–1684.


# A Soft Continuum Robot with Self-Controllable Variable Curvature

Xinran Wang1, Qiujie Lu1,2, Dongmyoung Lee1, Zhongxue Gan2, and Nicolas Rojas1
Manuscript received: September, 25, 2023; Revised December, 7, 2023; Accepted
January, 1, 2024.This paper was recommended for publication by Editor Yong-Lae
Park upon evaluation of the Associate Editor and Reviewers’ comments. This
work was supported in part by the Natural Science Foundation of China under
Grant 52305013, and by the Shanghai Pujiang Program (22PJD006).1REDS Lab,
Dyson School of Design Engineering, Imperial College London, 25 Exhibition
Road, London, SW7 2DB, UK (xinran.wang20, q.lu17, d.lee20, and
<EMAIL_ADDRESS>for Engineering and Technology, Fudan
University, 200433, Shanghai, China. (qj_lu<EMAIL_ADDRESS>
###### Abstract

This paper introduces a new type of soft continuum robot, called SCoReS, which
is capable of self-controlling continuously its curvature at the segment
level; in contrast to previous designs which either require external forces or
machine elements, or whose variable curvature capabilities are
discrete—depending on the number of locking mechanisms and segments. The
ability to have a variable curvature, whose control is continuous and
independent from external factors, makes a soft continuum robot more adaptive
in constrained environments, similar to what is observed in nature in the
elephant’s trunk or ostrich’s neck for instance which exhibit multiple
curvatures. To this end, our soft continuum robot enables reconfigurable
variable curvatures utilizing a variable stiffness growing spine based on
micro-particle granular jamming for the first time. We detail the design of
the proposed robot, presenting its modeling through beam theory and FEA
simulation—which is validated through experiments. The robot’s versatile
bending profiles are then explored in experiments and an application to grasp
fruits at different configurations is demonstrated. A narrated video detailing
the work can be seen at https://youtu.be/H6SCK0NjGpE.

###### Index Terms:

Continuous stiffness regulation, variable curvature, soft robot applications,
soft robot materials and design.

##  I Introduction

Soft continuum robots are in general flexible and adaptive. This feature
allows them to be more easily actuated with less actuators than degrees of
freedom and achieve more complex motions in constrained environments. However,
the common approach for soft continuum robots to have more bending profiles to
increase dexterity is adding more sections or segments [1][2][3][4]. This
strategy is effective but usually implies adding more actuators to control the
continuum robot, while increasing design complexity.

Exploring stiffness change inside soft robots has become a popular approach to
achieve a more dexterous workspace and higher external loads. The stiff-flop
manipulator, for instance, explored granular jamming using coffee powder to
allow stiffness change on silicone rubber-based soft continuum robots [5].
Tendon-driven continuum robots can also utilize granular jamming to alter its
curvature [6]. Other alternative options for achieving uniform stiffness
changes for each section of a soft continuum robot (i.e., the same stiffness
across a section) include fiber jamming [7], scale jamming [8], shape memory
alloy (SMA) [9], layer jamming [10], and pressurized fluidic compartments
[11].

Figure 1: (A) Illustration of controllable curvatures when the robot is
pressurized at 250kPa. (B) Inner structure: The soft continuum robot has a
hollow channel that allows the growing spine to travel inside. The jammed
spine then creates a high stiffness inside while the remaining length is
flexible. With the same control inputs, it creates different bending profiles.
(C) Realizing concept: The growing spine filled with granules changes length
using positive pressure and length control. Using negative pressure, it
increases stiffness by granular jamming.

Recently, there has been a greater focus on a non-uniform stiffness change
regardless of the number of actively actuated sections on the continuum robot,
allowing an even more dexterous motion thanks to the resulting variable
curvatures. Methods reported include using position locking with SMA, magnets,
screws, or tendons along the robots [12][13][14][15][16][17]. Other approaches
use pre-programmed springs or disks and add them to the continuum robot to
alter its local stiffness [18] [19]. These methods for achieving variable
curvatures have limitations due to the number of locking mechanisms or the
need for manual configuration in advance. Soft robots with continuously
changing stiffness are still an area of limited research, and current
approaches often rely on rigid materials such as rods or steel tubes with
patterns to alter stiffness [20][21].

Figure 2: (A) Fabrication of the soft continuum robot: silicone mixing,
degassing, winding strings on the rod, molding silicone for each tube,
assembling all the tubes for a single section continuum robot. (B) Dimension
of the soft continuum robot. It has nine chambers, and each three of them are
connected.

Soft continuum robots that take inspiration from vine robots [22], known for
their steering ability and growth capabilities, have recently been developed.
These continuum robots feature an additional controllable parameter—length
[23]. The benefit of such soft continuum robot comes from the material
everting from inside, allowing it to travel under different environments and
locations while being robust to friction. Much progress has been made to
adjust the stiffness change on such robots as well [24][25]. Utilizing layer
jamming on the robot for instance, stiffness can vary along the body [26].
Despite these advantages, these lightweight growing continuum robots are
designed for 2D ground navigation with limited 3D workspace. Our 3D soft
continuum robot with variable stiffness and curvature capability stands out
from other similar designs, as showcased in the advantages presented in Table
I.

In this paper, we then propose a Soft Continuum Robot with Self-controllable
variable curvature (SCoReS) that leverages the growing characteristics of vine
robots and micro-particle granular jamming effect, which is based on
lightweight glass bubbles, to create a new type of soft continuum robot.
Figure 1 presents the central concept to achieve variable bending curvature
profiles; The main contributions of this research work as follows:

  * •

We proposed a novel self-controllable variable curvature soft continuum robot
design, which allows continuous stiffness change at the segment level.

  * •

Develop a growing spine that can grow and adjust its stiffness based on micro-
particle granular jamming, along with a mechanism that allows the growing
spine to control its length automatically.

  * •

The bending profiles of SCoReS have been simulated, and experimentally
validated.

TABLE I: Comparison between existing soft continuum robots in 3D with variable stiffness or variable curvature capabilities | Variable Stiffness Capability
---|---
Research | Uniform | Localized | Continuous | Deployment
[5][6][7] | Y | N | N | Various
[11][13]
[18][19] | Y | Y | N | Pre-configure
Our work | Y | Y | Y | Self-controlled

##  II Robot Design

This section will discuss the comprehensive design of SCoReS robot. It has two
main components: one variable stiffness growing spine and one soft continuum
robot. The variable stiffness growing spine incorporates the ability to grow
inside the hollow continuum robot body and varying degrees of stiffness
automatically through granular jamming. By doing so, the single-section soft
continuum robot can achieve the effect of multiple bending curvatures,
allowing for greater flexibility and adaptability.

Figure 3: System detailed design (A) 1. Experimental Setup, with a motion
tracking system to track its end-effector positions and rotations. 2. Cross-
section area to show the inner structure of the robot. Inside the steel box,
airtight fabric divides the volume into two. The bottom volume contains glass
bubbles and a length control mechanism. The filter paper on top of the
pegboard prevents glass bubbles from going into the vacuum chamber when
jamming. 3. Detailed length control mechanism: using a stepper motor for
precise length control. (B) The state transitions to allow the SCoReS to
reconfigure.

###  II-A Continuum Robot Design

The continuum robot comprises nine air chambers, with three of them
interconnected for seamless three-dimensional movement. An electric
proportional pressure valve regulates the pressure of each chamber.
Furthermore, the robot features an inner channel that enables a growing robot
to pass through.

We chose DragonSkin 10 (Smooth-on Inc.) as the material for our soft continuum
robot, as it offers excellent stiffness levels comparable to EcoFlex 30
(Smooth-on Inc.) , and similar elongation properties.

To prevent radical expansion of the air chamber, we have implemented separate
constraints for each chamber. This way, the soft continuum robot can maintain
a consistent inner diameter, regardless of its configuration, whether it’s
bent or straight. This feature ensures that the soft continuum robot has
minimal inner radius change.

The process of creating the soft continuum robot is illustrated in Figure
2(A). To begin, we mix two components of DragonSkin 10 silicone and put them
in a vacuum container to degas. Then, we wind kevlar strings tightly around an
aluminum rod in Step III. Next, we pour silicone onto the outside of the
aluminum rod in Step IV. Once the silicone is cured, we remove the aluminum
rod, leaving the kevlar strings inside the molded silicone. In Steps V and VI,
we place a smaller aluminum rod in the center of the shell and pour silicone
into the mold. After curing, we remove the small center aluminum rod and are
left with the silicone air tube and perfectly embedded Kevlar strings.

By integrating nine circular air channels into a new mold and merging them
together with silicone, we can produce a soft continuum robot with an inner
hollow channel. This process was previously only used for small surgical
applications [27], but we have now successfully implemented it in a larger
continuum robot measuring 40cm in length. Our design features an inner channel
with an inner diameter of 5.8cm and an outer diameter of 10cm as the dimension
of the continuum robot is in Figure 2(B).

###  II-B Design of Variable Stiffness Growing Spine

Creating multiple curvatures and providing variable stiffness in a single-
section continuum robot has always been challenging.

We are proposing a design of a self-controllable growing spine with variable
stiffness inside an enclosed steel box. Our project first upgrades growing
robots with a new length control unit for accurate position control.

Figure 3(3) clearly illustrates the various components of the length control
mechanism. These include a pair of shafts with gears, a pre-tensioned
torsional spring holder for storing the growing robot material, and a NEMA23
stepper motor that comes with a planetary gearbox and an L-shaped connector.

By applying positive pressure to the chamber, we create a constant exerting
force on the air-tight ripstop TPU-coated fabric of the growing robot. As the
stepper motor rotates, the fabric is released through the fixed-diameter
geared shafts, enabling the growing spine to grow from the inside out without
any high friction. By regulating the rotation of the stepper motor, we can
effortlessly control the length of the growing robot (as seen in Video from
0:0:14s to 0:0:53s).

The material is wrapped around a pre-tensioned torsional spring, which allows
for retraction of the fabric while the stepper motor pulls it back. The step
motor has been applied with air-tight glue and wrapped in TPU-coated fabric to
prevent particle obstruction during operation. This mechanism ensures that the
spine grows smoothly and consistently.

Glass bubble material has been utilized to enhance fluidity and prevent volume
shrinking in the jamming process when applied with growing robots. Unlike
conventional granular jamming techniques that use large and heavy particles,
glass bubble consists of small and lightweight (0.2g/cc) round micro-size
particles with a hollow spherical structure that promotes easy flow and
movement while Capable of withstanding intense pressure during the jamming
process [28].

In Figure 3(A-2), the steel box is divided into two separate volumes by a
flexible TPU-coated ripstop fabric. As shown in the pink section of Figure
3(A-2), the pink volume is either pressurized or connected to the atmosphere.
The blue-colored volume contains a length control mechanism, which is
displayed in Figure 3(A-3), along with glass bubble materials. To secure the
length control mechanism in place, a pegboard is installed in the blue volume.
Further, layers of industrial-grade filter paper and steel nets are placed on
top of the pegboard to prevent the glass bubble from falling out.

Finally, the green volume serves as a vacuum chamber, it connects with blue
volume and growing spine, ensuring the safe and particle-free operation of the
vacuum pump. When performing granular jamming, the growing spine stiffness
will be significantly increased.

Figure 4: (A) Control diagram of the overall system. (B) Detailed physical
components: Electric Proportional Regulators, Step Motor Driver, Digital
Analog Converter (DAC) for Proportional Regulators, Arduino and Solenoid
Valves .

###  II-C Overall System and Self-Controllable workflow

Figure 3(A-1) displays the overall experimental setup. The soft continuum
robot is positioned beneath the top-growing robot, and a laser photoelectric
switch is placed between them for precise position calibration during the
growing and retracting of the top-growing robot. The Opti-track motion
tracking System is employed to track the end-effector constantly, providing
real-time position and rotation information.

SCoReS undergoes three states in sequence: idle, growing, and jammed. In the
idle state, the entire volume inside the cubic steel box is connected to the
atmosphere. The robot neither grows nor retracts during this state.

The steel box in its cubic shape, is fully pressurized to 10kPa to facilitate
the growth of the growing spine. As the spine grows, gravity helps the glass
bubble to occupy all the empty spaces inside the growing spine. The internal
pressure of the spine prevents any buckling during retraction. Furthermore,
the fabric is pressurized to distribute glass bubbles evenly across the
growing spine, ensuring they are not concentrated at the bottom. This allows
more space for the glass bubble to move freely inside the growing spine.

When in the jammed state, the cubic steel box’s top volume is connected to the
atmosphere, while the rest of the volume is vacuumed through a vacuum pump.
This effectively jams the growing robot in place and significantly increases
its stiffness. When we control the soft continuum robot, it will bend in
various curvatures due to alterations in its stiffness (as seen in the
supplementary video from 0:0:54s to 0:01:25s).

The SCoReS has the capability to adjust stiffness at various lengths and
produce multiple bending curvatures by following the three states in order and
utilizing different inner lengths during growth.

###  II-D Robot Control

Figure 4(A) illustrates that the control command is sent to the
microcontroller via serial communication from a computer. To control of the
bending of the soft continuum robot in three dimensions, we have employed
three electric-proportional regulators.

Four solenoid valves are in the system to control the two separate volume
pressures inside the steel box. Valve 1 is linked to the atmosphere, while
valve 2 is connected to an air pump via a proportional pressure regulator.
Valve 3 is linked to the pink volume in the steel box, and valve 4 is linked
to a vacuum pump.

When in the idle state, the steel box is connected to the atmosphere through
the open of solenoid valve 1 while valve 2 remains closed. Valve 3 is open,
valve 4 is closed.

During the growing state, the solenoid valve 1 and 4 are closed while valve 2
and 3 are open. The proportional regulator is adjusted to maintain a pressure
of 10kPa to pressurize the steel box.

In the jammed state, the solenoid valve 1 is open to connect to the atmosphere
while valve 2 and valve 3 are closed. Valve 4 is open to vacuum the steel box.

The computer receives pressure feedback from a pressure sensor measure range
from -100kPa to 100kPa (SMC PSE543A-R06). The electric proportional pressure
regulators we used are SMC ITV 2030.

##  III Modelling and Simulation

###  III-A Stiffness Modelling of Jammed growing Spine

In order to accurately assess the stiffness of the jammed growing spine, we
have employed the Euler-Bernoulli beam theory for mathematical modeling
purposes. The beam model has been employed previously to check the variable
stiffness capability in granular jamming [29]. Presently, we are evaluating
this model on the jammed growing spine and integrating it with Finite Element
Analysis for a comprehensive system evaluation.

The jammed growing spine has an airtight fabric outside and inside filled with
material glass bubble. The jammed growing spine’s stiffness can be represented
as a cantilever beam problem for characterization purposes. In this scenario,
one end of the beam is fixed while a point force acts upon the free end. The
fundamental equation in beam theory serves as the basis for this analysis:

$M=-E\times I\times\frac{d^{2}y}{dx^{2}}$ (1)

Where: $M$ is the bending moment. $E$ is the Young’s Modulus. $I$ is the
moment of inertia of the beam’s cross-sectional area. $y$ represents the
deflection in the direction perpendicular to the beam length. $x$ represents
the position length from the fixed end along the beam.

The bending moment, in this case, is caused by a point load on the free end,
for a cantilever length of $L$, is:

$M=-F\times(L-x)$ (2)

Now combining equation 1 and 2, we get relationship between point load $F$ and
the deflection $y$:

$\frac{d^{2}y}{dx^{2}}=\frac{F\times(L-x)}{E\times I}$ (3)

Having two integrations over position $x$, and equations are constricted by
two boundary conditions, on the fixed end, the $\frac{dy}{dx}=0$ and when
$x=0,y=0$, the equation can be simplified as:

$y=\frac{F\times(3L-x)x^{2}}{6\times E\times I}$ (4)

The jammed growing spine has a circular shape, $I=\frac{\pi r^{4}}{4}$,
rearranging equation 4 and substituting $I$, we can get $E$ when position is
equal to the length of the beam $x=L$ :

$E=\frac{4\times F\times L^{3}}{3\times\pi r^{4}\times y}$ (5)

Figure 5: (A) Experimental setup for bending experiments with load cell
mounted on UR5. (B) Experimental setup for measuring length control with
transparent tube and ruler. Figure 6: (A) Bending Experiment during loading
and unloading phrase for jammed growing robot at -70kPa and various lengths.
(B) Calculated Young’s modulus in kPa for the jammed growing spine at various
lengths.

The overall experimental setup to measure the stiffness is provided in the
5(A). For a jammed growing spine measured in length $L$, we have used a
Universal Robots UR5 to push a distance of $Y$ in the direction perpendicular
to the length direction of the growing robot. This model provides a
fundamental basis for the following simulation analysis.

TABLE II: Length Growing Experiments Reference (cm) | 5 | 10 | 15 | 20 | 25 | 30
---|---|---|---|---|---|---
Experiment (cm) | 5.25 | 10.53 | 15.77 | 21.06 | 26.42 | 31.43
Std Dev (cm) | 0.30 | 0.32 | 0.31 | 0.34 | 0.39 | 0.36

The growing spine is able to grow at different lengths. To characterize the
stiffness property of the jammed top-growing robot, we discretize the
continuous length control problem into several quantitative segments in length
for better comparison and understanding. We decided to evaluate the properties
of SCoReS at these growing robot length segments: 5cm, 10cm, 15cm, 20cm, 25cm,
and 30cm.

Figure 7: FEA simulation and experimental results comparison for the SCoReS at
different configurations. (A) The SCoReS overall bending profile comparison
for different configurations at the pressure of 250kPa (B) Box plot for
bending angle values comparison for the SCoReS pressurized at 50 kPa, 100 kPa,
150 kPa, 200 kPa, and 250 kPa with jammed growing spine segments ranging from
0 cm to 30 cm.

For a jammed growing spine growing at different length $L$, UR5 push the
jammed robot for a certain distance, and the load cell (ROBOTiq ft300-s)
records the forces during the loading and unloading phase. The jammed spine is
at -70kPa pressure. Figure 6(A) shows the bending experiments results for the
growing spine at different lengths. Based on the equation 5, we are able to
obtain the stiffness parameter for the jammed growing spine at different
length segments.

The averaged stiffness values at different positions are $0.318e^{3}$ kPa,
$1.323e^{3}$ kPa, $2.032e^{3}$ kPa, $3.069e^{3}$ kPa, $3.763e^{3}$ kPa,
$4.389e^{3}$ kPa for length segments of 5cm, 10cm, 15cm, 20cm, 25cm and 30cm.
We can find as the length gets longer, the stiffness property also increases
from Figure 6(B). The main reason is that gravity also has an effect on the
overall stiffness of the jammed-growing robot. The SCoReS is at a large scale.
The length difference here is significant. While the glass bubble material is
light (0.2g/cc), the pressure at the bottom of the growing robot increases as
the length gets longer. This, in return, increases the overall stiffness
value.

We have also evaluated the length control accuracy of the growing spine
through experiments. The setup in Figure 5(B) was utilized for testing the
length control accuracy at corresponding positions, and the reported results
are listed in Table II. The average error was found to be $5.18\%$, while the
maximum standard deviation error was 0.39 cm.

Figure 8: Experiment results for bending length difference for SCoReS robot of
different growing spine configurations (0cm, 5cm, 10cm, 15cm, 20cm, 25cm and
30cm growing spine length)

in A) X-axis direction, B) Y-axis direction.

###  III-B Finite Element Analysis of the Robot

Figure 9: Demonstrations of the SCoReS in various configurations in real-world
applications.: A) 20cm jammed growing spine pressurized from 250kPa to 0kPa.
B) 0cm jammed growing spine pressurized from 0kPa to 250kPa. C) Comparing a
20cm configuration pressurized at 250kPa to a 0cm configuration pressurized at
200kPa shows a difference in bending angle. D) Comparing a 20cm configuration
pressurized at 250kPa to a 0cm configuration pressurized at 250kPa shows a
difference in bending length.

We incorporated the stiffness property of the jammed growing robot into our
Finite Element Analysis (FEA) to model the interaction between the soft
continuum robot’s bending and the jammed growing robot.

In the previous section, we modeled the jammed top-growing robot as a
cantilever beam. Using this modeling information, we create a cylindrical
shape with the same diameter as the jammed-growing robot and place it inside
the soft continuum robot. The chamber of the soft continuum robot is
constrained similarly to the real-world case with kevlar strings inside the
FEA simulation. The coefficient we used for DragonSkin 10 material is Neo-
Hookean Model in Abaqus with $C_{1}$ = 42.5 kPa [30][31].

As shown in Figure 7(A), the cross-section area of the SCoReS has a cylinder
with different lengths inside the soft continuum robot in FEA. The image’s
color coding reveals the position displacement of the elements. At the 30cm
configuration with a jammed growing spine inside, we observed a significant
reduction in displacement and a shift in bending angle. With a portion of the
soft continuum robot remaining stiff and the rest flexible, the inner jammed
growing spine’s length can be configured to alter the SCoReS’s overall shape.
This design makes the robot highly versatile, offering multiple different
bending curvature profiles compared to the traditional single-section
continuum robot.

##  IV Experiment Evaluation

###  IV-A Bending Angle

In order to validate our mathematical modeling and corresponding FEA
simulation, we have conducted a comparison with our experimental results. As
depicted in Figure 7, our analysis focuses on the SCoReS with the inner jammed
growing spine configuration at various lengths ranging from 0cm to 30cm. By
comparing the bending angle from our experimental results with the FEA
simulation, we have observed a close match between the two. Specifically, we
found that the bending shape in FEA closely resembles the results of real
experiments, particularly when the SCoReS is pressurized from 50kPa to 200kPa.
It provides us with confidence in the accuracy of our mathematical modeling
using beam theory and corresponding FEA simulation.

Upon closer inspection, we noticed that the maximum error occurred at 250kPa
when the SCoReS’s growing spine at 0cm. FEA predicted a bending angle of 75.48
degrees, while the experiment measured 65.64 degrees, with a maximum error of
only 14.99%. The error occurs without the interference of the jammed growing
spines. Mentioning the material coefficients for DragonSKin 10 from the
literature has a discrepancy with the real-world case.

In real experiments, the SCoReS robot bent 41.50 degrees at 30cm and 250kPa
pressure, compared to 65.64 degrees at 0cm configuration, a difference of
24.13 degrees. These findings provide us insight into the varying bending
angles that the SCoReS robot can offer under different configurations.

###  IV-B End Effector Position

We have further examined the SCoReS robot bending profile with respect to the
end-effector positions. Figure 8 illustrates the variations in the SCoReS
robot extending length in both the X-axis and the Y-axis for different
configurations. In particular, we observed that in the X-axis, there is a
difference of 10.36cm in the distance between the robot with the inner growing
spine at 0cm segment (20.87cm) and the robot at 30cm segment (10.51cm) when at
250kPa.

In the Y-axis direction, there is a variation of 3.57 cm at 250kPa between the
robot at 0cm configuration and the robot at 30cm configuration. Although this
is less than the changes observed in the X-axis, the robot at the
configuration of 30cm exhibits a much more stable motion in the Y-axis
direction being pressurized from 0kPa to 250kPa.

###  IV-C Real World Application Demo

Furthermore, we showcase a demo of a real-world application by mounting a
gripper on the SCoReS robot. As illustrated in Figure 9, a clear comparison is
made between the robot with the inner growing spine at 0cm segment (Figure
9(A)) and the robot at the configuration of 20cm (Figure 9(B)). It is evident
from Figure 9(A) that half of the continuum robot remains rigid in contrast to
Figure 9(B) even when both are pressurized to 250kPa. Moreover, the end-
effector path highlights the distinct bending profiles.

In Figure 9(C), we compare the SCoReS robot with 0 cm configuration
pressurized to 200kPa and the robot with 20cm pressurized to 250kPa. The
robot’s end effectors remain in a similar position, but the robot displays a
more versatile bending profile with a 12-degree divergence in bending angle.
This proves the effectiveness of the robot design. Last but not least, at a
pressure of 250kPa, Figure 9(D) shows varying extension lengths for the two
mentioned configurations, making it valuable in restricted spaces.

##  V Conclusion

In this study, we proposed the design of SCoReS that can be self-controllable
to achieve a versatile bending profile. The jammed growing spine design has
both characteristics of growing and jamming, making it unique in the
application of variable curvature soft continuum robots. We used beam theory
to model the jammed growing spine and conducted FEA analysis on the SCoReS.
Our approach proved to be appropriate, as demonstrated by the close match
between FEA and real experiments. Our experiments also showed significant
changes in the bending profile, including bending angles and end-effector
positions. In addition, The SCoReS’s real-world demonstrations illustrate its
ability to execute complex motions, such as varying the bending angle at the
same end-effector location and shortening its bending length within a
constrained environment. In the future, we plan to test the capability of the
robot when jammed in different bending configurations instead of straight, and
improve the design by adding fabric sleeves to anchor transportation units on
the SCoReS. We’ll also explore control strategies based on the calibrated
silicone elastic models to fully utilize its bending profiles in constrained
environments.

## References

  * [1] J. Santoso and C. D. Onal, “An origami continuum robot capable of precise motion through torsionally stiff body and smooth inverse kinematics,” Soft Robotics, vol. 8, no. 4, pp. 371–386, 2021.
  * [2] W. McMahan, B. A. Jones, and I. D. Walker, “Design and implementation of a multi-section continuum robot: Air-octor,” in 2005 IEEE/RSJ international conference on intelligent robots and systems, pp. 2578–2585, IEEE, 2005.
  * [3] X. Wang and N. Rojas, “A data-efficient model-based learning framework for the closed-loop control of continuum robots,” in 2022 IEEE 5th International Conference on Soft Robotics (RoboSoft), pp. 247–254, IEEE, 2022.
  * [4] M. B. Wooten and I. D. Walker, “Environmental interaction with continuum robots exploiting impact,” IEEE Robotics and Automation Letters, vol. 7, no. 4, pp. 10136–10143, 2022.
  * [5] M. Cianchetti, T. Ranzani, G. Gerboni, T. Nanayakkara, K. Althoefer, P. Dasgupta, and A. Menciassi, “Soft robotics technologies to address shortcomings in today’s minimally invasive surgery: the stiff-flop approach,” Soft robotics, vol. 1, no. 2, pp. 122–131, 2014.
  * [6] W. R. Wockenfuß, V. Brandt, L. Weisheit, and W.-G. Drossel, “Design, modeling and validation of a tendon-driven soft continuum robot for planar motion based on variable stiffness structures,” IEEE Robotics and Automation Letters, vol. 7, no. 2, pp. 3985–3991, 2022.
  * [7] L. Arleo, L. Lorenzon, and M. Cianchetti, “Variable stiffness linear actuator based on differential drive fiber jamming,” IEEE Transactions on Robotics, 2023.
  * [8] S. H. Sadati, Y. Noh, S. Elnaz Naghibi, K. Althoefer, and T. Nanayakkara, “Stiffness control of soft robotic manipulator for minimally invasive surgery (mis) using scale jamming,” in Intelligent Robotics and Applications: 9th International Conference, ICIRA 2015, Portsmouth, UK, August 24–27, 2015, Proceedings, Part III, pp. 141–151, Springer, 2015.
  * [9] S. Jiang, B. Chen, F. Qi, Y. Cao, F. Ju, D. Bai, and Y. Wang, “A variable-stiffness continuum manipulators by an sma-based sheath in minimally invasive surgery,” The International Journal of Medical Robotics and Computer Assisted Surgery, vol. 16, no. 2, p. e2081, 2020.
  * [10] A. B. Clark and N. Rojas, “Malleable robots: Reconfigurable robotic arms with continuum links of variable stiffness,” IEEE Transactions on Robotics, vol. 38, no. 6, pp. 3832–3849, 2022.
  * [11] J. Fras and K. Althoefer, “A fluidic actuator with an internal stiffening structure inspired by mammalian erectile tissue,” in 2023 IEEE International Conference on Robotics and Automation (ICRA), pp. 662–668, IEEE, 2023.
  * [12] C. Pogue, P. Rao, Q. Peyron, J. Kim, J. Burgner-Kahrs, and E. Diller, “Multiple curvatures in a tendon-driven continuum robot using a novel magnetic locking mechanism,” in 2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 472–479, IEEE, 2022.
  * [13] P. Wang, S. Guo, X. Wang, and Y. Wu, “Design and analysis of a novel variable stiffness continuum robot with built-in winding-styled ropes,” IEEE Robotics and Automation Letters, vol. 7, no. 3, pp. 6375–6382, 2022.
  * [14] K. Nakano, M. Gunji, M. Ikeda, K. Or, M. Ando, K. Inoue, H. Mochiyama, K. Nakajima, R. Niiyama, and Y. Kuniyoshi, ““robostrich” manipulator: A novel mechanical design and control based on the anatomy and behavior of an ostrich neck,” IEEE Robotics and Automation Letters, 2023.
  * [15] C. Yang, S. Geng, I. Walker, D. T. Branson, J. Liu, J. S. Dai, and R. Kang, “Geometric constraint-based modeling and analysis of a novel continuum robot with shape memory alloy initiated variable stiffness,” The International Journal of Robotics Research, vol. 39, no. 14, pp. 1620–1634, 2020.
  * [16] P. Rao, C. Pogue, Q. Peyron, E. Diller, and J. Burgner-Kahrs, “Modeling and analysis of tendon-driven continuum robots for rod-based locking,” IEEE Robotics and Automation Letters, 2023.
  * [17] E. Amanov, T.-D. Nguyen, and J. Burgner-Kahrs, “Tendon-driven continuum robots with extensible sections—a model-based evaluation of path-following motions,” The International Journal of Robotics Research, vol. 40, no. 1, pp. 7–23, 2021.
  * [18] J. Zhang, Y. Li, Z. Kan, Q. Yuan, H. Rajabi, Z. Wu, H. Peng, and J. Wu, “A preprogrammable continuum robot inspired by elephant trunk for dexterous manipulation,” Soft Robotics, vol. 10, no. 3, pp. 636–646, 2023.
  * [19] K. Ma, X. Chen, J. Zhang, Z. Xie, J. Wu, and J. Zhang, “Inspired by physical intelligence of an elephant trunk: Biomimetic soft robot with pre-programmable localized stiffness,” IEEE Robotics and Automation Letters, vol. 8, no. 5, pp. 2898–2905, 2023.
  * [20] J. Kim, W.-Y. Choi, S. Kang, C. Kim, and K.-J. Cho, “Continuously variable stiffness mechanism using nonuniform patterns on coaxial tubes for continuum microsurgical robot,” IEEE Transactions on Robotics, vol. 35, no. 6, pp. 1475–1487, 2019.
  * [21] B. Zhao, L. Zeng, Z. Wu, and K. Xu, “A continuum manipulator for continuously variable stiffness and its stiffness control formulation,” Mechanism and Machine Theory, vol. 149, p. 103746, 2020.
  * [22] E. W. Hawkes, L. H. Blumenschein, J. D. Greer, and A. M. Okamura, “A soft robot that navigates its environment through growth,” Science Robotics, vol. 2, no. 8, p. eaan3028, 2017.
  * [23] J. D. Greer, T. K. Morimoto, A. M. Okamura, and E. W. Hawkes, “A soft, steerable continuum robot that grows via tip extension,” Soft robotics, vol. 6, no. 1, pp. 95–108, 2019.
  * [24] R. Jitosho, S. Simón-Trench, A. M. Okamura, and B. H. Do, “Passive shape locking for multi-bend growing inflated beam robots,” in 2023 IEEE International Conference on Soft Robotics (RoboSoft), pp. 1–6, IEEE, 2023.
  * [25] S. Wang, R. Zhang, D. A. Haggerty, N. D. Naclerio, and E. W. Hawkes, “A dexterous tip-extending robot with variable-length shape-locking,” in 2020 IEEE International Conference on Robotics and Automation (ICRA), pp. 9035–9041, IEEE, 2020.
  * [26] B. H. Do, V. Banashek, and A. M. Okamura, “Dynamically reconfigurable discrete distributed stiffness for inflated beam robots,” in 2020 IEEE International Conference on Robotics and Automation (ICRA), pp. 9050–9056, IEEE, 2020.
  * [27] J. Fraś, J. Czarnowski, M. Maciaś, J. Główka, M. Cianchetti, and A. Menciassi, “New stiff-flop module construction idea for improved actuation and sensing,” in 2015 IEEE international conference on robotics and automation (ICRA), pp. 2901–2906, IEEE, 2015.
  * [28] S. E. Bakarich, R. Miller, R. A. Mrozek, M. R. O’Neill, G. A. Slipher, and R. F. Shepherd, “Pump up the jam: Granular media as a quasi-hydraulic fluid for independent control over isometric and isotonic actuation,” Advanced Science, vol. 9, no. 15, p. 2104402, 2022.
  * [29] A. Jiang, T. Ranzani, G. Gerboni, L. Lekstutyte, K. Althoefer, P. Dasgupta, and T. Nanayakkara, “Robotic granular jamming: Does the membrane matter?,” Soft Robotics, vol. 1, no. 3, pp. 192–201, 2014.
  * [30] P. Polygerinos, K. Galloway, Z. Wang, F. Connolly, J. Overvelde, and H. Young, “Fiber reinforced actuators: Finite element modeling,” 2023.  Accessed: 2023-09.
  * [31] M. S. Xavier, A. J. Fleming, and Y. K. Yong, “Finite element modeling of soft fluidic actuators: Overview and recent developments,” Advanced Intelligent Systems, vol. 3, no. 2, p. 2000187, 2021.