Planar segment processes with reference mark distributions, modeling and estimation

TL;DR

This paper introduces statistical methods for modeling planar segment processes with reference distributions, combining parametric and non-parametric estimation techniques, and demonstrates their effectiveness through simulation studies.

Contribution

It develops a general theoretical framework and two specific models for segment processes with reference distributions, including estimation procedures and simulation validation.

Findings

Estimation methods effectively recover reference distributions.

Simulation studies show estimators' variability and accuracy.

Models accommodate inhomogeneous and Gibbs-type segment processes.

Abstract

The paper deals with planar segment processes given by a density with respect to the Poisson process. Parametric models involve reference distributions of directions and/or lengths of segments. These distributions generally do not coincide with the corresponding observed distributions. Statistical methods are presented which first estimate scalar parameters by known approaches and then the reference distribution is estimated non-parametrically. Besides a general theory we offer two models, first a Gibbs type segment process with reference directional distribution and secondly an inhomogeneous process with reference length distribution. The estimation is demonstrated in simulation studies where the variability of estimators is presented graphically.

Tables

Table 1. Table 1: Empirical mean, standard deviation (sd) and coefficient of variation (CV) of Takacs-Fiksel estimates of scalar parameters in the model having density ( 12 ) with reference directional distribution. It is based on 100 simulations, the two cases correspond to a = − 0.5 , a = − 3 . formulae-sequence 𝑎 0.5 𝑎 3 a=-0.5,\;a=-3.

	true	mean	sd	CV
$a$	-0.5	-0.496	0.071	0.14
$\tau$	1000	1011	154.7	0.15
	true	mean	sd	CV
$a$	-3	-3.03	0.356	0.12
$\tau$	1000	976	141.0	0.14

Table 2. Table 2: Empirical mean, standard deviation (sd) and coefficient of variation (CV) of the maximum likelihood estimates of scalar parameters in the model ( 21 ) with reference length density. Based on a sample of 60 simulated realizations of the segment process with parameters b = 3 , τ = 900 , α = 2 , β = 4 , e a = 1 . formulae-sequence 𝑏 3 formulae-sequence 𝜏 900 formulae-sequence 𝛼 2 formulae-sequence 𝛽 4 subscript 𝑒 𝑎 1 b=3,\tau=900,\alpha=2,\beta=4,e_{a}=1.

	true	mean	sd	CV
$b$	3	3.051	0.481	0.16
$\tau$	900	948.9	186.2	0.19

Equations

$$\mathrm{d}{P}_X(\mathbf{x})=p(\mathbf{x})\mathrm{d}{P}_{\eta }(\mathbf{x}),\mathbf{x}\in \mathbf{N},$$

$${\int }_{\mathbf{N}}p(\mathbf{x})\mathrm{d}{P}_{\eta }(\mathbf{x})=1.$$

$$\mathbb{P}(X\in A)=$$

$$=e^{-\lambda (B)}\sum _{n=0}^{\infty }\frac{1}{n!}{\int }_B\cdots {\int }_{B}1[(x_1,\dots ,x_n)\in A]p(x_1,\dots ,x_n)\lambda (\text{d}{x}_1)\dots \lambda (\text{d}{x}_n),$$

$${\lambda }^{\ast }(\mathbf{x},u)=\frac{p(\mathbf{x}\cup u)}{p(\mathbf{x})},\mathbf{x}\in \mathbb{N},u\notin \mathbf{x}.$$

$$\mathbb{E}\left[\sum _{u\in X}q(u,X\setminus u)\right]={\int }_B\mathbb{E}[{\lambda }^{\ast }(X,u)q(u,X)]\lambda (\text{d}u),$$

$$Y=B\times S,S=(0,b]\times [0,\pi ),$$

$$\lambda (\text{d}(y,r,\phi ))=\frac{1}{b\pi }\mathrm{d}y\text{d}r\text{d}\phi ,$$

$$p(\mathbf{x})=c{\tau }^{n(\mathbf{x})}\prod _{u\in \mathbf{x}}g(r,\phi )\exp (bD(\mathbf{x})+aE(\mathbf{x})),\mathbf{x}\in \mathbf{N},$$

$$D(\mathbf{x})=\sum _{u\in \mathbf{x}}d(u),E(\mathbf{x})=\sum _{\{u,v\}\subset \mathbf{x}}^{\ne }\Phi (u,v)$$

$${\lambda }^{\ast }(\mathbf{x},u)=\tau g(r,\phi )\exp \left(bd(u)+a\sum _{v\in \mathbf{x}\setminus \{u\}}\Phi (u,v)\right).$$

$$\mathbb{E}X(A)={\int }_A\rho (u)\text{d}u.$$

$$\rho (u)=\mathbb{E}{\lambda }^{\ast }(X,u).$$

$$\kappa (A)=\mathbb{E}\sum _{y\in \nu }{1}_A(y).$$

$$C(A\times D)=\mathbb{E}\sum _{(y,\xi )\in X}{1}_A(y)1_D(\xi ).$$

$${\omega }_D(A)={\int }_A{P}^y(D)\kappa (\text{d}y),$$

$$f_X^{(y)}(\xi )=\frac{\rho (y,\xi )}{{\int }_S\rho (y,\xi )\text{d}\xi }.$$

$$C(A\times D)={\int }_Y{1}_A(y)1_D(\xi )\rho (y,\xi )\text{d}\xi \text{d}y.$$

$$\kappa (A)={\int }_A{\int }_S\rho (y,\xi )\text{d}\xi \text{d}y.$$

$${\int }_A{P}^y(D){\int }_S\rho (y,\xi )\text{d}\xi \text{d}y={\int }_A{\int }_D\rho (y,\xi )\text{d}\xi \text{d}y$$

$$P^y(D)=\frac{{\int }_D\rho (y,\xi )\text{d}\xi }{{\int }_S\rho (y,\xi )\text{d}\xi }$$

$$Y=B\times [0,\pi ).$$

$$\lambda (\text{d}y)=\text{d}y\frac{1}{\pi }\text{d}\phi$$

$$p(\mathbf{x})=c\exp (aN(\mathbf{x})){\tau }^{n(\mathbf{x})}\prod _{u_i\in \mathbf{x}}g({\phi }_i),\mathbf{x}\in \mathbf{N},$$

$$\Phi (u,v)=1_{[u\uparrow v]},$$

$${\lambda }^{\ast }(\mathbf{x},u)=\tau g(\phi )\exp (a{N}_{\mathbf{x}}(u)),$$

$$g(\varphi )=h(\kappa )\exp (\kappa \cos (2(\varphi -\mu ))),\varphi \in [0,\pi ],h(\kappa )=\frac{1}{\pi I_0(\kappa )},$$

$$p(\mathbf{x})=c(\theta )\exp (⟨\theta ,G(\mathbf{x})⟩),\mathbf{x}\in \mathbf{N},$$

$$\theta =(a,\log (\tau h(\kappa )),\kappa ),$$

$$G(\mathbf{x})=(N(\mathbf{x}),n(\mathbf{x}),\sum _{x_i\in \mathbf{x}}\cos (2({\phi }_i-\mu ))),$$

Full text

Keywords: Conditional intensity, segment process, semiparametric estimation

AMS subject classification: 60D05, 60G55

1 Introduction

The present research addresses an important problem in the statistics of spatial marked point processes given by a density with respect to the Poisson process. Observing a realization of spatial data which shall be fitted to such a model we first estimate the parameters by a method of point estimation. However, among the quantities to be estimated there may appear also the reference mark distribution which need not coincide with the observed mark distribution of the process. Both the scalar parameters and the reference mark distribution are needed e.g. when we try to simulate the model. This distribution can be also parametrized by a subjective choice of model, cf. [4]. In the present paper the main aim is to estimate the reference mark distribution non-parametrically, i.e. in total to use a semiparametric approach instead of a fully parametric one.

An early paper [1] mentions parameter estimation of a marked point process by means of the maximum pseudolikelihood method but the authors do not identify our problem. Much more attention is paid to it in [12] where the marks form radii of circles centered at the points of a point process given by a density with respect to the Poisson process. The random set corresponding to the union of circles in a compact window is investigated. Since an exact method is not available the authors use an approximation what means that estimation of the distribution of radii is done by methods for a Boolean model. Then an MCMC maximum likelihood method ([13]) is used for the estimation of parameters of the point process. A recent paper by [9] deals with the same model as [12], their goal is to use the Takacs-Fiksel estimator instead of the computationally demanding maximum likelihood method.

In our work we present a solution of the problem for another random set. We deal with the planar segment process ([6], [15]) having a density of exponential form with respect to a Poisson process. We consider first a model with reference directional distribution and in the end a model with reference length distribution of segments. The difference in comparison to marked point processes presented in the literature is that when the directional distribution is present we are not completely on Euclidean spaces. Our main tool is the derived relation between the reference and observed mark distribution for the segment process. The basic asymptotic properties of the Takacs-Fiksel method of estimation are known, see [7]. In the present paper we pay attention to the computing of semiparametric estimators using simulated data and quantifying the variability for small sample size.

First some background from spatial point processes having a density is presented. Then in Section 3 a general formula for the mark distributions is derived. In Section 4 we present a Gibbs type segment process with reference directional distribution and fixed lengths of segments. Interactions enter the model by means of intersections, cf. [19], [18]. An approximation from [3] is used to avoid the problem with unknown moments. The Takacs-Fiksel estimator for this model is developed in Section 5 and tested on data simulated by MCMC algorithm from [10]. The estimator of the reference directional distribution is evaluated numerically. In Section 6 we address an old problem of the existence of a stationary process with given conditional intensity from the previous Section. Section 7 presents a model with reference length distribution of segments while the reference directional distribution is uniform. This model is an inhomogeneous Poisson process with a condition on segments to lie entirely in a circular window. The maximum likelihood method of estimation is available, see Section 8, we use the isotropy of the process to simplify the computation.

2 Spatial point process given by a density

Consider a bounded Borel set $B\subset{\mathbb{R}}^{d}$ with Lebesgue measure $|B|>0$ and a measurable space $({\mathbf{N}},{\mathcal{N}})$ of integer-valued finite measures on $B.$ ${\mathcal{N}}$ is the smallest $\sigma$-algebra which makes the mappings $\mathbf{x}\mapsto\mathbf{x}(A)$ measurable for all Borel sets $A\subset B.$ A random element having a.s. values in $({\mathbf{N}},{\mathcal{N}})$ is called a finite point process. Let a Poisson point process $\eta$ on $B$ have finite intensity measure $\lambda$ with no atoms and distribution $P_{\eta}$ on $\mathcal{N}.$ We consider a finite point process $X$ on $B$ given by a density $p$ w.r.t. $P_{\eta},$ i. e. with distribution $P_{X}$

[TABLE]

where $p:{\mathbf{N}}\rightarrow{\mathbb{R}}_{+}$ is measurable satisfying

[TABLE]

As described in [2], p. 61, integer-valued finite measures can be represented in this context by $n$-tuples of points corresponding to their support ($n$ is variable). We have

[TABLE]

[TABLE]

$A\in\mathcal{N}.$ For an hereditary density $p$ the distribution of the process is alternatively determined by the conditional intensity

[TABLE]

An important tool is the Georgii-Nguyen-Zessin (GNZ) formula

[TABLE]

valid for any measurable test function $q$ on $B\times{\mathbf{N}}.$

3 Segment process with reference mark distribution

In the paper we study random segment processes in the plane $\mathbb{R}^{2}.$ A segment is a closed set which will be parametrized by its centre $y=(y_{1},y_{2})\in\mathbb{R}^{2},$ length $r>0$ and direction $\varphi\in[0,\pi).$ There is a bijection $\iota$ between the parametric space $\{u=(y,r,\varphi)\}$ and the system of segments as a subsystem of the space of closed sets in $\mathbb{R}^{2}.$ Throughout the paper we use exclusively the parametric representation of segments, omitting the bijection $\iota$ in some expressions, e.g. $z\in u$ means a point $z$ of the segment $u,$ etc.

A segment process can be considered as a marked point process with two marks corresponding to the length and direction of a segment. Let $B\subset\mathbb{R}^{2}$ be bounded measurable,

[TABLE]

where $b>0$ is an upper bound for the segment length, $[0,\pi)$ is the manifold of axial directions. Further $({\mathbf{N}},{\mathcal{N}})$ is a measurable space of integer-valued finite measures on $Y.$ Let the Poisson process $\eta$ on $Y$ have intensity measure $\lambda,$

[TABLE]

where on the right hand side we have a multiple of Lebesgue measure on $Y.$ Let the segment process $X$ have an hereditary density $p$ with respect to $\eta$

[TABLE]

where $c>0$ is a normalizing constant, $b\in\mathbb{R},$ $a\leq 0,$ and $\tau>0$ are parameters, $g$ is the density of reference length-direction distribution, $n(\mathbf{x})$ is number of segments in the configuration $\mathbf{x}.$ Further

[TABLE]

(sum over pairs of different segments from $\mathbf{x}$) for some nonnegative functions $d,\;\Phi$ on $Y,\;Y^{2},$ respectively, $\Phi$ is called the pair potencial, $E$ energy function in the theory of Gibbs processes. The conditional intensity

[TABLE]

Let $\rho$ be the intensity function of the process $X,$ i.e. for $A\subset Y$

[TABLE]

Using $q(u,\mathbf{x})={\bf 1}_{[u\in A]}$ in the GNZ formula we obtain

[TABLE]

For $A\subset B$ measurable and the point process $\nu$ of segment centres, denote

[TABLE]

Let $D\subset S,$

[TABLE]

The measure $\omega_{D},\;\omega_{D}(A)=C(A\times D),$ is absolutely continuous with respect to $\kappa$ and the Radon-Nikodym density $P^{y}(D),$ such that

[TABLE]

is the Palm mark distribution of $X$ at $y.$ Let $f_{X}^{(y)}$ be the density of the Palm mark distribution of length and direction of a typical segment of the process $X$ at the location $y\in B.$

Proposition 1

For each $y\in B$ we have

[TABLE]

Proof: For Borel sets $A\subset B,\;D\subset S$ it holds using the Campbell theorem

[TABLE]

Specially for $D=S$ we have

[TABLE]

Then from (9) we have

[TABLE]

for all Borel sets $A\subset B,$ finally

[TABLE]

and (10) follows. $\Box$

4 Segment process with reference directional distribution

In this Section we assume that the segment length $r>0$ is fixed and for a bounded Borel set $B\subset{\mathbb{R}}^{2}$ we deal with the

[TABLE]

A segment $u=(y,\varphi)\in Y$ has center $y$ and axial orientation $\varphi.$ Consider a measurable space $({\mathbf{N}},{\mathcal{N}})$ of integer-valued finite measures on $Y,$ alternatively of the supports of these measures.

We deal with the unit Poisson segment process $\eta$ with the intensity measure

[TABLE]

on $Y.$ Let the segment process $X$ have a density $p$ with respect to $\eta,$ we consider

[TABLE]

where $N(\mathbf{x})$ is the total number of intersections between segments, $g$ a reference probability density on $[0,\pi),$ $\varphi_{i}$ direction of $i-$th segment $u_{i},$ $a\leq 0,\;\tau>0$ are parameters, $c$ a normalizing constant. That means in the general model (6) we have $b=0,$ $g$ depends only on $\varphi$ and the pair potential is

[TABLE]

i.e. indicator of the event that segment $u$ hits segment $v.$ This model belongs to the more general class of facet processes ([20]). The conditional intensity is

[TABLE]

where $N_{\mathbf{x}}(u)$ is the number of segments of $\mathbf{x}\setminus\{u\}$ hit by the segment $u.$

Example: Let the direction density $g$ be that of von Mises distribution on $[0,\pi)$ with parameters $\kappa\geq 0,\;\mu\in\mathbb{R},$

[TABLE]

$I_{0}(\kappa)$ is the modified Bessel function of the first kind and order 0. Then

[TABLE]

where

[TABLE]

[TABLE]

segments $x_{i}=(u_{i},\varphi_{i})$ have centres $u_{i}\in B$ and directions $\varphi_{i}\in[0,\pi).$ The normalizing constant $c(\theta)$ is defined as

[TABLE]

Also

[TABLE]

is the largest set of $\theta$ such that the density (15) is well defined.

In the present paper we want to relax the unimodality assumption of the reference directional distribution and deal with a general density $g.$ Let $f_{X}^{(y)}$ be the density of the Palm mark distribution $P^{y}$ of direction of a typical segment of the process $X$ at the location $y\in B.$ For $u=(y,\varphi)$ we get from (14) that

[TABLE]

is the intensity function of $X,$ cf. (8). The expectation $\mathbb{E}_{X}$ with respect to the distribution of $X$ is not analytically tractable, therefore [3] suggest an approximation

[TABLE]

where $\eta(\rho)$ is a Poisson process with intensity function $\rho.$

Proposition 2

For $u=(y,\varphi)\in Y$ we have

[TABLE]

where $M_{u}=\{x\in Y;\;x\uparrow u\}$.

Proof: Using (2) we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

$\Box$

From Proposition 1 we have

[TABLE]

for normalizing constants $C_{y},\;y\in B.$ Assuming that there exists a stationary segment process in $\mathbb{R}^{2}$ with given conditional intensity (this time the conditional intensity cannot be defined by means of densities, but from the energy function) corresponding to (14), we have that $C=C_{y},\;f_{X}^{(y)}=f_{X}$ do not depend on $y$ in the extension of $B$ onto the whole $\mathbb{R}^{2}.$ Under this assumption, discussed in Section 6, we estimate ($\beta$ is the direction of segment $x$)

[TABLE]

Using (16) we can then express the desired density $g$ approximately as

[TABLE]

where

[TABLE]

5 Semiparametric estimation using Takacs-Fiksel approach

In this section we suggest a method of estimation of parameters $C,a,\tau$ and the density $g$ from the previous section using the Takacs-Fiksel method. From formula (3) we obtain innovation equations

[TABLE]

and solve them for various test functions $q.$ We take $\lambda^{*}$ from (14) where we insert approximation (18) for unknown $g.$ First the density $f_{X}$ is estimated using a kernel estimator for directional data [11]. Then put

[TABLE]

and we estimate $C,a$ from the system of Takacs-Fiksel equations:

[TABLE]

[TABLE]

Here in the innovations equations we take score functions $q(u,\mathbf{x})=N_{\mathbf{x}}(u)$ and $q=1,$ respectively, the integrals in the second term of (19) are evaluated by Monte Carlo method using $J$ independent simulations of segments uniformly distributed in $Y.$ Then we plug the estimators of $C,a$ in a formula obtained by integrating (18):

[TABLE]

and finally estimate $g$ from (18).

A numerical study is based on twice 100 simulated realizations of segment process with parameters $\kappa=1,\mu=0,\tau=1000,r=0.12$ on $[0,1]^{2}\times[-\frac{\pi}{2},\frac{\pi}{2}]$. The two cases I, II investigated are $a=-0.5,\;a=-3,$ respectively. The results are in Table 1, Fig. 1 and Fig. 2.

In Table 1 we observe a small difference between the true and mean values for both estimates of $a$ and $\tau$. The coefficient of variation

[TABLE]

is also comparable despites the fact that the model II involves more interactions (inhibition of intersections) than the model I. In Fig. 1 we can observe how the kernel estimator of the observed (Palm mark) directional distribution differs from the true reference directional distribution (von Mises). The results in Fig. 2 suggest that the estimate of the reference density is slightly better (smaller bias and variability) for the case I than for the case II. We conclude that the approximation (17) works well in the Takacs-Fiksel method here.

6 Existence of a stationary process with given conditional intensity

In order to be able to use the approximation (18) correctly we need a sufficient condition for the existence of a stationary Gibbs segment process in $\mathbb{R}^{2}$ with prescribed conditional intensity. Various conditions are present in the literature starting with [16]. We shall use a recent work of [8] who deals with the concept of an hereditary energy function, invariant with respect to shifts, satisfying the finite range property. While his results are formulated for point processes in $\mathbb{R}^{d},$ what we need here from [8] is valid also for particle processes in the sense of [17]. This straightforward extension of a part of [8]from Gibbs point process to Gibbs particle process is presented in [14], specially also for segment processes in the plane. The energy function is hereditary if for each $\mathbf{x}\in{\mathbf{N}},\;y\in supp\,\mathbf{x}$ we have

[TABLE]

where $supp$ is the support of a measure.

Definition 1

A function $f:{\mathbf{N}}\rightarrow\mathbb{R}$ is local (on $A$) if there exists $A\subset\mathbb{R}^{d}$ bounded such that for all $\mathbf{x}\in{\mathbf{N}}$ it holds

[TABLE]

where $\mathbf{x}_{A}$ is the restriction of $\mathbf{x}$ onto $A.$

Definition 2

An energy function $H$ has finite range $R>0$ if for all $A\subset\mathbb{R}^{d}$ bounded the local energy

[TABLE]

is a local function on $A\oplus b(0,R).$

Proposition 3

The energy function $E$ from (7) with pair potential (13) has finite range property with $R=r.$

Proof: Let $A\subset\mathbb{R}^{d}$ be bounded, $\mathbf{x}\in{\mathbf{N}}.$ We have to show that

[TABLE]

i.e.

[TABLE]

This is true since to both sides exactly intersections of such pairs of segments contribute, which have either both centres in $A$ or one centre in $A$ and the other in $(A\oplus b(0,r))\cap A^{c}.$ $\Box$

Corollary 1

There exists a stationary segment process in $\mathbb{R}^{2}$ with conditional intensity (14).

Proof: The energy function $E$ from (7) with pair potential (13) is nonnegative, invariant with respect to shifts, hereditary and has finite range. According to [8] and [14] the existence of a stationary segment process is guaranteed. $\Box$

The uniqueness issue is more complex, see [8] for Gibbs point process in $\mathbb{R}^{d}.$ Here the uniqueness is not investigated, we deal with an existing stationary process in order to use the approximation (18).

7 The segment process with reference length distribution

Consider a circle $B\subset{\mathbb{R}}^{2}$ centered in the origin with diameter $e_{a}>0.$ Let $L_{o}=[0,e_{a}]$ be the interval of segment lengths. Then

[TABLE]

is the space of segments $u=(y,r,\varphi)\in Y$ which have centre $y,$ length $r$ and axial direction $\varphi.$

We consider the Poisson segment process $\eta$ with the intensity measure

[TABLE]

on $Y.$ Let the segment process $X$ have a density $p$ w.r.t. $\eta:$

[TABLE]

$b\in\mathbb{R},\,\tau>0,\;c$ is the normalizing constant, $r_{i}$ is the length of $i-$th segment $u_{i},\;f_{1}$ is a reference probability density on $L_{o}$ and

[TABLE]

That means in the general model (6) we have $a=0,$ $g=f_{1}$ depends only on the length variable $r,$ marginal reference directional distribution is uniform. Moreover a factor is added in (21) which forces the segments to lie completely in the window.

For $u\subset B$ the quantity $d(u)\in(0,\frac{1}{2}]$ is the normalized distance of the most distant point of $u$ from the centre of $B.$ For positive, negative values of parameter $b$ more, less distant segments from the origin prevail, respectively, cf. Fig. 3.

The corresponding conditional intensity is

[TABLE]

For the intensity function it holds

[TABLE]

There are no interactions among the segments in the model (21) with statistics $D(\mathbf{x}),$ $X$ is in fact an inhomogeneous Poisson process with unknown reference density $f_{1}$ and a condition $X\subset B.$ The process $X$ is isotropic since both the reference Poisson process and the density $p$ are invariant with respect to rotations around the origin. We denote $Ay$ rotation of $y\in\mathbb{R}^{2},$ where $A$ is an orthonormal matrix $2\times 2,$ analougously $Au$ is rotation of the segment $u$ and finally we use the same symbol $A\varphi$ for the corresponding direction of the segment $u$ after rotation $Au.$ Then also the intensity function is invariant with respect to rotations, i.e.

[TABLE]

for all $u\in Y$ and all rotations $A.$

Let $f_{X}^{(y)}$ be the bivariate (Palm) density of the distribution of length and direction of the segment centered at $y,$ from Proposition 1 we have for $y\in B$ a normalizing constant $C_{y}>0$ such that

[TABLE]

From the isotropy of $X$ it holds

[TABLE]

and therefore

[TABLE]

for all $y\in B$ and all rotations $A.$

8 Semiparametric estimation using the maximum likelihood method

The aim is to estimate parameters $b,\tau$ and density $f_{1}$ from simulated data. We are using again a semiparametric approach so that $f_{1}$ is not parametrized. Because of inhomogeneity of $X$ the solution has to be discretized, but we make use of isotropy. In the parametric part, since the process is Poisson we use maximum likelihood method for parameter estimation. For an observed realization $\mathbf{x}$ the likelihood is defined as

[TABLE]

The logarithmic likelihood

[TABLE]

has to be maximized with respect to $\tau,b.$ We have

[TABLE]

[TABLE]

Using (22) we obtain equations

[TABLE]

[TABLE]

In the estimation procedure we proceed in several steps:

(i) Consider discrete levels of $||y||$ like $0<w_{1}<\dots<w_{k}<e_{a}/2$ and put $y_{j}=(0,w_{j})^{T}\in B.$ For $\Delta_{f}=\frac{e_{a}}{2k}$ let

[TABLE]

[TABLE]

The kernel estimator of the bivariate densities $f_{X}^{(y_{j})},\;j=1,\dots,k,$ is evaluated from the observed data in each class, i.e. from the sample of segments $u_{j}^{(i)}=(y_{j}^{(i)},r_{j}^{(i)},\varphi_{j}^{(i)}),\;i=1,\dots,m_{j}$ centered in $B_{j},\;j=1,\dots,k.$ Because of (23) these segments are first rotated by such $A_{j}^{(i)}$ that

[TABLE]

for some $q>0.$ Then we apply kernel estimation to each sample

[TABLE]

to estimate the length-direction density $f_{X}^{(y_{j})},\;j=1,\dots,k.$ Since the length component has values on a compact $L_{o},$ we use here a system of beta kernels [5]. The angular component is a circular variable, cf. [11], jointly we use a bivariate product kernel.

(ii) Since the length of the longest possible segment centered in $j-$th class is

[TABLE]

we use the step $\Delta_{c}^{(j)}=l_{j}/m$ for numerical integration of (22) with $y=y_{j}$ and a fixed $\varphi$ using Simpson rule, to express unknown constants $C_{y_{j}}=C_{j},\;j=1,\dots,k,$ as functions of $b.$

(iii) Integrals in the equations (24) are evaluated by Monte Carlo simulation of $M$ segments in $Y$ uniformly randomly. Let $n$ of them lie completely in $B,$ denoted $\bar{u}_{j}^{(i)}=(\bar{y}_{j}^{(i)},\bar{r}_{j}^{(i)},\bar{\varphi}_{j}^{(i)}),\;i=1,\dots,n_{j}$ where $n_{j}$ of them are centered in $B_{j},\;\sum_{j=1}^{k}n_{j}=n.$ Then

[TABLE]

is an equation with a single variable $b$ which is solved numerically.

(iv) Having estimated $b$ we obtain $C_{j}$ from step (ii) and then $\tau$ from any of the equations

[TABLE]

[TABLE]

which are Monte Carlo analogues of the equations in (24).

(v) Finally the estimator of the reference length density $f_{1}$ is obtained by plugging these estimators in (22) and renormalizing. We choose $\varphi=0,$ from the levels of $y$ those close to the boundary are mostly variable and therefore not used. We can average results of several levels of $y.$

A numerical study is based on a sample of 60 simulated independent realizations of the segment process with parameters $b=3,\tau=900,\alpha=2,\beta=4,e_{a}=1.$ In step (i) of the estimation procedure $k=6$ classes were considered. In Fig. 4 there are kernel estimates of the observed length density in all six classes separately. We can observe the difference between the estimated observed length distribution and the true reference beta distribution, which increases along classes towards the boundary of the window. The results of estimation of parameters $b,\tau$ (step (iii), (iv)) are in Table 2. We observe still a reasonably small difference between the true and the mean value of estimates of $b$ and $\tau,$ but it is larger than in the Gibbs model with directional distribution in Section 5, also CV is slightly larger. This can be justified by the fact that the present model is inhomogeneous and the estimation procedure is more complex. The semiparametric estimator of the reference length density $f_{1}$ (step (v)) is in Fig. 5, here we do not consider two outer classes where the inhomogeneity is the highest (since there are little admissible directions of longer segments close to the border). We observe a small bias of the estimator of the reference density, but the variability is larger than in the first model in Section 5 from the same reasons as we argued to Table 2.

Acknowledgements

The research was supported by the Czech Science Foundation, project 16-03708S and by Charles University, grant SVV-2016-260334.