Multivariate geostatistical mapping of radioactive contamination in the Maddalena Archipelago

18 Pages·2011·1.03 MB·English

Working Paper Series
Department of Economics
University of Verona
Multivariate geostatistical mapping of radioactive contamination
in the Maddalena Archipelago (Sardinia, Italy)
Marco Minozzo
WP Number: 21 December 2011
ISSN: 2036-2919 (paper), 2036-4679 (online)
Advances in Statistical Analysis manuscript No.
(will be inserted by the editor)
Multivariate Geostatistical Mapping of Radioactive
Contamination in the Maddalena Archipelago
(Sardinia, Italy) – SPATIAL SPECIAL ISSUE
Marco Minozzo · Clarissa Ferrari
Received:date/Accepted:date
Abstract To improve the quality of prediction of radioactive contamination,
geostatistical methods, and in particular multivariate geostatistical models,
are increasingly being used. These methods, however, are optimal only in the
case in which the data may be assumed Gaussian and do not properly cope
withdatameasurementsthatarediscrete,nonnegativeorshowsomedegreeof
skewness. To deal with these situations, here we consider a hierarchical model
in which non-Gaussian variables of different kind are handled simultaneously.
Weshowthatwhenobservationsareassumedtobeconditionallydistributedas
PoissonandGamma,variogramsandcross-variogramshaveconvenientsimple
forms, and estimation of the parameters of the model can be carried out by
Monte Carlo EM. This work has been inspired by radioactive contamination
data from the Maddalena Archipelago (Sardinia, Italy).
Keywords Generalizedlinearmixedmodel,linearmodelofcoregionalization,
Markov chain Monte Carlo, Monte Carlo EM, spatial factor model
Mathematics Subject Classification (2010) MSC 62M30 · MSC 62H11
1 Introduction
Theanalysisofspatiallycorrelateddatahasreceivedinrecentyearsincreasing
attention both from an applied as well as from a methodological perspective,
MarcoMinozzo
DepartmentofEconomics,UniversityofVerona,
Viadell’Artigliere,19
Tel.:[phone]
E-mail:[email protected]
ClarissaFerrari
DepartmentofEconomics,UniversityofVerona,
Viadell’Artigliere,19
Tel.:[phone]
E-mail:[email protected]
2 MarcoMinozzo,ClarissaFerrari
and environmental studies have played in this respect a leading role with the
assessmentandmonitoringofair,soilandwatercontamination,inwhichmea-
sures are collected over specific geographic areas. In these studies, to improve
theaccuracyofprediction,geostatisticalmethods,andinparticularmultivari-
ate geostatistical models and techniques in which the additional information
embodied in correlated variables can be exploited, have often been used. Tra-
ditional multivariate geostatistical methods are provided by cokriging and by
thelinearmodelofcoregionalization(withitssimplercounterpart,thepropor-
tional covariance model, otherwise known as intrinsic correlation model) to-
getherwiththerelated‘factorialkriginganalysis’(Matheron1982;Goovaerts,
Sonnet and Navarre 1993; Wackernagel 1995; Chil`es and Delfiner 1999). For
instance, recent works concerned with the use of the variogram to study the
spatial autocorrelation present in radioactive contamination data are Dubois
and De Cort (2001), Dunois and Bossew (2003), and Dewiere et al. (2004).
Works dealing with kriging prediction are instead Simon, Bouville and Beck
(2004), and O’Dea and Dowdall (1999) which consider kriging prediction of
logarithmictransformeddata.Othergeostatisticalapproachesfortheanalysis
ofmultivariatespatialdatahavebeenproposedbyVerHoefandBarry(1998)
andChristensenandAmemiya(2001,2002),whereasaBayesianversionsofthe
linear model of coregionalization have been considered by Gelfand, Schmidt
and Sirmans (2002), Schmidt and Gelfand (2003), and Banerjee, Carlin and
Gelfand (2004). The widespread use of the classical linear model of coregion-
alization for the modelling of multivariate spatial data is mainly due to the
convenient characterization of the multivariate spatial autocorrelation struc-
ture, as well as to the availability of widely accepted estimation procedures
based on the (eigenvalue) principal component decomposition and on the al-
gorithm of Goulard and Voltz (1992). An advantage of this model is that no
hypotheses are made about the distribution of the data and only variograms
and cross-variograms are needed. The drawback is that the commonly used
estimation and prediction procedures may be considered adequate only in the
case in which the data can be assumed to be multivariate Gaussian. In pres-
ence of count or skewed measurements, as in contamination data, the use of
this model can lead to misleading predictions and to erroneous conclusions
about the underlying factors representing pollution levels.
In the case of a single spatial measurement, a preferable approach to krig-
ing in presence of non-Gaussian data, in particular contamination count data,
has been suggested by Diggle, Moyeed and Tawn (1998) which propose to
model data of this kind following a hierarchical approach (see also Alexander,
Moyeed and Stander 2000, Cressie 2000, and Christensen and Waagepetersen
2002). On the other hand, here, to cope with situations in which more spatial
variables are gathered simultaneously, we consider, following the suggestions
put forward in the univariate case, a hierarchical modeling framework (see
also Minozzo 2002; Minozzo and Fruttini 2004) based on a generalization of
the classical geostatistical proportional covariance model. In this modeling
framework, whereas the underlying (unobservable) factors are assumed to be
multivariate Gaussian, the actual data may be assumed to be (conditionally)
MultivariateGeostatisticalMappingofRadioactiveContamination 3
distributed as, for example, Poisson, Negative Binomial, Skew-Normal, etc.,
allowing non-Gaussian variables of different kind to be dealt with simultane-
ously.
Interestingly, assuming for the observed data, Poisson and Gamma dis-
tributions, to cope with discrete and nonnegative measurements, in Section
2 we show that variograms and cross-variograms, even among variables of
different kind, have convenient simple forms. Adopting a non-Bayesian infer-
entialframework,inSection3weshowhowestimatesoftheparametersofthe
model can be obtained by employing the method of moments and standard
factor analysis techniques. Moreover, we also show how likelihood inference
can be performed by using stochastic versions of the EM algorithm, like the
stochastic EM (StEM) algorithm or the Monte Carlo EM (MCEM) algorithm
(see, for instance, Nielsen 2000, and Fort and Moulines 2003). Prediction of
the underlying factors can then be carried out by Markov chain Monte Carlo
(MCMC) algorithms. In Section 4 an application of this modeling framework
is considered in the analysis of radioactive contamination data from the Mad-
dalena Archipelago which originally motivated our study. Lastly, in Section 5
we report some conclusions.
2 The hierarchical spatial factor model
Let y (x ), i = 1,...,m, k = 1,...,K, be a set of geo-referenced data mea-
i k
surementsrelativetomregionalizedvariables,gatheredatK spatiallocations
x . These m regionalized variables are seen as a partial realization of a set of
k
m random functions Y (x), i=1,...,m, x∈D ⊆R2. For these functions we
i
assume, for any x, and for i6=j,
Y (x)⊥⊥Y (x) | Z (x), and Y (x)⊥⊥Z (x) | Z (x),
i j i i j i
′ ′′
and, for x 6=x , and j =1,...,m,
′ ′′ ′ ′ ′′ ′
Y (x)⊥⊥Y (x ) | Z (x), and Y (x)⊥⊥Z (x ) | Z (x),
i j i i j i
where Z (x), i=1,...,m, are mean zero joint stationary Gaussian processes.
i
Moreover, we assume that, for any given i and x, conditionally on Z (x), the
i
random variables Y (x) have conditional distributions f (y;M (x)), specified
i . i i
by the values of the conditional expectations M (x) = E Y (x) Z (x) , and
i i i
that hi(Mi(x)) = βi +Zi(x), for some parameters βi and(cid:2)some(cid:12)(cid:12)know(cid:3)n link
functions h (·).
i
For the latent part of the model, we adopt the following structure. For the
m joint stationary Gaussian processes Z (x), let us assume the linear factor
i
model
P
Z (x)= a F (x)+ξ (x), (1)
i ip p i
Xp=1
where a are m × P coefficients, F (x), p = 1,...,P, are P ≤ m non-
ip p
observable spatial components (common factors) responsible for the cross
4 MarcoMinozzo,ClarissaFerrari
correlation between the variables Z (x), and ξ (x) are non-observable spa-
i i
tial components (unique factors) responsible for the residual autocorrelation
in the Z (x) unexplained by the common factors. We assume that F (x) and
i p
ξ (x) are mean zero stationary Gaussian processes with covariance functions
i
Cov F (x),F (x+h) =ρ(h), and Cov ξ (x),ξ (x+h) =ψ ρ(h),
p p i i i
(cid:2) (cid:3) (cid:2) (cid:3)
where ρ(h) is a real spatial autocorrelation function common to all factors
such that ρ(0) = 1 and ρ(h) → 0, as ||h|| → ∞, and ψ are non-negative
i
real parameters. We also assume that the processes F (x) and ξ (x) have all
p i
cross-covariances identically equal to zero.
As far as the latent processes Z (x), i=1,...,m, are concerned, we have
i
that
Cov[Z (x),Z (x+h)]=E[Z (x)·Z (x+h)]−E[Z (x)]·E[Z (x+h)]
i j i j i j
P P
=E a F (x)+ξ (x) · a F (x+h)+ξ (x+h)
ip p i jp p j
(cid:20)(cid:18) (cid:19) (cid:18) (cid:19)(cid:21)
Xp=1 Xp=1
P
=E a a F (x)F (x+h)+ξ (x)ξ (x+h)
ip jp p p i j
(cid:20) (cid:21)
Xp=1
P
= a a E[F (x)F (x+h)]+E[ξ (x)ξ (x+h)],
ip jp p p i j
Xp=1
that is, we have (taking i = j) the covariance functions Cov[Z (x),Z (x+
i i
h)] = P a2 ρ(h)+ψ ρ(h), and (taking i 6= j) the cross-covariance func-
p=1 ip i
tions CPov[Z (x),Z (x+h)] = P a a ρ(h). Of course, taking h = 0, we
i j p=1 ip jp
also have that Var[Z (x)] = PP a2 + ψ , and that Cov[Z (x),Z (x)] =
i p=1 ip i i j
P a a . P
p=1 ip jp
P
2.1 Poisson distribution
To completely specify the model we need to detail the conditional distribu-
tions f (y;M (x)) and the link functions h (·). For a fixed i, that is, for a
i i i
fixed random function Y (x), let us now assume that f (y;M (x)) are Poisson
i i i
distributions with means M (x), that is, that
i
f (y;M (x))=exp(−M (x))·(M (x))y/y!, y =0,1,2,..., (2)
i i i i
and that h (·) = ln(·). Then it is possible to show that Y (x) is second order
stationary.iInfact,definingς2 =. Var[Z (x)],τ∗ =. exp{β +i ς2/2}andτ2(x)=.
i i i i i i
Var[Y (x)|Z (x)], and considering that, for every x ∈ D, the expectations
i i
MultivariateGeostatisticalMappingofRadioactiveContamination 5
E[exp{Z (x)}] = exp{ς2/2} and E[exp{2Z (x)}] = exp{2ς2}, and also that,
i i i i
for the Poisson distribution, τ2(x)=M (x), with some algebra we obtain:
i i
. ∗
µ =E[Y (x)]=E[E[Y (x)|Z (x)]]=E[exp{β +Z (x)}]=τ ;
i i i i i i i
σ2 =. Var[Y (x)]=E[τ2(x)]+Var[M (x)]=τ∗+(τ∗)2[exp{ς2}−1];
i i i i i i i
C (h)=. Cov[Y (x),Y (x+h)]=(τ∗)2[exp{ς2ρ(h)}−1], h6=0. (3)
ii i i i i
Itiseasytocheckthat,inthecaseinwhichρ(h)=ρ(−h),wehaveC (−h)=
. ii
C (h). Moreover, C (∞) = 0 and C (0) = σ2 6= C (0+), that is, the covari-
ii ii ii i ii
ance function C (h) is discontinuous at the origin.
ii
As for the variogram of the observable Y (x) we have, for h6=0,
i
. 1 1
γ (h)= Var Y (x+h)−Y (x) = E (Y (x+h)−Y (x))2 (4)
ii i i i i
2 2
(cid:2) (cid:3) (cid:2) (cid:3)
1 1
= Var Y (x) + Var Y (x+h) −Cov Y (x), Y (x+h)
i i i i
2 2
(cid:2) (cid:3) (cid:2) (cid:3) (cid:2) (cid:3)
1 1 1 1
= E τ2(x) + E τ2(x+h) + Var M (x) + Var M (x+h)
2 i 2 i 2 i 2 i
(cid:2) (cid:3) (cid:2) (cid:3) (cid:2) (cid:3) (cid:2) (cid:3)
−E M (x)·M (x+h) +E M (x) ·E M (x+h) (5)
i i i i
=τ∗+(cid:2) (τ∗)2 exp ς2 −(cid:3)exp (cid:2)ς2ρ(h)(cid:3) . (cid:2) (cid:3) (6)
i i i i
(cid:2) (cid:8) (cid:9) (cid:8) (cid:9)(cid:3)
Let us note that in this case the variogram of Y (x) is always equal to zero
i
at the origin with nugget γ (0+) = τ∗ that equals the mean µ and it is, in
ii i i
general, strictly greater than 0. It is also easy to check that the sill of the
.
variogram satisfies γ (∞)=C (0)=σ2.
ii ii i
Considering two random functions Y (x) and Y (x), with i 6= j, for which
i j
f (y;M (x))andf (y;M (x))arePoissondistributionswithmeansM (x)and
i i j j i
M (x), and h (·) = h (·) = ln(·), and considering that E[exp{Z (x)+Z (x+
j i j . i j
h)}]=exp{ς2/2+ς2/2+ς ρ(h)},whereς =Cov[Z (x),Z (x)],wehavethat
i j ij ij i j
. ∗ ∗
C (h)=Cov[Y (x),Y (x+h)]=τ τ [exp{ς ρ(h)}−1]. (7)
ij i j i j ij
Note that C (h) = C (h) and that, if ρ(h) = ρ(−h),we have C (−h) =
ij ji ij
C (h). From (7) it follows that the random functions Y (x) and Y (x) are
ij i j
jointly second order stationary.
6 MarcoMinozzo,ClarissaFerrari
Fig.1 Theoreticalvariogramγii(h)andcross-variogramγij(h)oftheobservableYi(x)and
Yj(x)assumingaconditionalPoissondistribution,alogarithmiclinkfunction,andautocor-
relationfunctionρ(h)=exp{−((α||h||)δ)}.Variousshapescanbeobtainedbyappropriately
choosingvaluesfortheparametersβi,βj,ςi2,ςj2,ςij,andfortheparametersα,δ.
Moreover, the cross-variogram between Y (x) and Y (x) is given by
i j
. 1
γ (h)= Cov(Y(x+h)−Y(x)),(Y(x+h)−Y(x))
ij i i j j
2
(cid:2) (cid:3)
1
= E(Y(x+h)−Y(x))(Y(x+h)−Y(x))
i i j j
2
(cid:2) (cid:3)
1 1
= CovM (x+h),M (x+h) − CovM (x+h),M (x)
i j i j
2 2
(cid:2) (cid:3) (cid:2) (cid:3)
1 1
− CovM (x),M (x+h) + CovM (x),M (x)
i j i j
2 2
(cid:2) (cid:3) (cid:2) (cid:3)
1 1 1
= EM(x)M(x) − EM(x)EM(x) + EM(x+h)M(x+h)
i j i j i j
2 2 2
(cid:2) (cid:3) (cid:2) (cid:3) (cid:2) (cid:3) (cid:2) (cid:3)
1 1
− EM(x+h)EM(x+h) − EM(x)M(x+h)
i j i j
2 2
(cid:2) (cid:3) (cid:2) (cid:3) (cid:2) (cid:3)
1 1 1
+ EM(x)EM(x+h) − EM(x+h)M(x) + EM(x+h)EM(x)
i j i j i j
2 2 2
∗ ∗(cid:2) (cid:3) (cid:2) (cid:3) (cid:2) (cid:3) (cid:2) (cid:3) (cid:2) (cid:3)
=τ τ exp ς −exp ς ρ(h) . (8)
i j ij ij
(cid:2) (cid:8) (cid:9) (cid:8) (cid:9)(cid:3)
Noticethatforthecross-variogramwehaveγ (0)=γ (0+)=0andγ (∞)=
ij ij ij
∗ ∗
τ τ [exp{ς }−1] = C (0), which can be smaller or greater than zero, de-
i j ij ij
pendingonthesignofς .Thetheoreticalvariogramγ (h)andthetheoretical
ij ii
cross-variogram γ (h) are depicted in Fig. 1.
ij
MultivariateGeostatisticalMappingofRadioactiveContamination 7
2.2 Gamma distribution
Alternatively to the Poisson distribution, let us now assume that, for a given
i and for any given x, conditionally on Z (x), the random variables Y (x) are
i i
Gamma distributed with conditional expectations
.
M (x)=E Y (x) Z (x) =exp β +Z (x) =νb,
i i i i i
(cid:2) (cid:12) (cid:3) (cid:8) (cid:9)
(here h (·)=ln(·)) and conditi(cid:12)onal variances
i
Var Y (x) Z (x) =νb2 =ν−1exp 2β +2Z (x) =ν−1(M (x))2,
i i i i i
(cid:2) (cid:12) (cid:3) (cid:8) (cid:9)
where ν >0 and(cid:12) b>0 are parameters, that is, let us assume that, for y >0,
yν−1 yν ν ν
f (y;M (x))= exp − . (9)
i i
Γ(ν) (cid:26) M (x)(cid:27)(cid:18)M (x)(cid:19)
i i
Here the ‘shape’ parameter ν is constant for x ∈ D, whereas the ‘scale’ pa-
rameterbvariesoverDdependingontheconditionalexpectationM (x).Con-
i
sidering ν constant, f (y;M (x)) belongs to the (one parameter) exponential
i i
family with natural parameter −1/b.
EveninthiscaseitispossibletoshowthatY (x)issecondorderstationary.
. ∗ i
Indeed, µ =E[Y (x)]=τ ,
i i i
σ2 =. Var[Y (x)]= 1e−2βi(τ∗)4+(τ∗)2[exp{ς2}−1],
i i ν i i i
.
and C (h)=Cov[Y (x),Y (x+h)], for h6=0, is the same as in (3).
ii i i
Then, considering (5), the variogram of the observable Y (x) is given by,
i
for h6=0,
1
γ (h)= E (Y (x+h)−Y (x))2
ii i i
2
(cid:2) (cid:3)
1 1 1 1
= E e2βi+2Zi(x) + E e2βi+2Zi(x+h)
2 (cid:20)ν (cid:21) 2 (cid:20)ν (cid:21)
1 1
+ Var eβi+Zi(x) + Var eβi+Zi(x+h)
2 2
h i h i
−E eβi+Zi(x)·eβi+Zi(x+h) +E eβi+Zi(x) ·E eβi+Zi(x+h)
h i h i h i
1
= e−2βi(τ∗)4+(τ∗)2 exp ς2 −exp ς2ρ(h) . (10)
ν i i i i
(cid:2) (cid:8) (cid:9) (cid:8) (cid:9)(cid:3)
It is immediate to see that whereas the nugget is equal to ν−1e−2βi(τ∗)4, the
i
second addendum is the same as in (6) in the Poisson case.
For the cross-covariance C (h) and for the cross-variogram γ (h), when
ij ij
Y (x) and Y (x) are both conditionally Gamma distributed, or when one is
i j
Gamma and the other is Poisson, we obtain the same forms (7) and (8) as
in the case in which Y (x) and Y (x) are both Poisson, and thus Y (x) and
i j i
Y (x)areagainjointlysecondorderstationary.Indeed,theformsofthecross-
j
covariance and of the cross-variogram depend on the link function which we
assumed to be ln(·) both in the Poisson and in the Gamma case.
8 MarcoMinozzo,ClarissaFerrari
3 Estimation and prediction
Adopting a non-Bayesian inferential framework, a simple procedure to obtain
a set of preliminary estimates for the parameters of the model can be carried
out by implementing some version of the classical method of moments, or
by using some iterative least squares algorithm, similarly to the procedures
considered, for instance, in Wackernagel (1995) (see also Goulard and Voltz
1992) for the linear model of coregionalization. A computationally convenient
estimation procedure based on the method of moments is as follows. Let us
assume, for instance, that m=2 and that Y (x), and Y (x) are conditionally
1 2
distributed as Poisson and Gamma, respectively, with h (·) = ln(·). For a
i
specific choice of the ‘shape’ parameter ν, let us first consider the estimation
ofthem+m+m(m−1)/2parametersβ ,ς2 andς .SinceY (x),i=1,2,are
i i ij i
jointly second order stationary, we are allowed to consider the sample means
and sample variances and covariances of the observable Y (x )
i k
K K
y¯ =. 1 y (x ), s2 =. 1 y (x )−y¯ 2, i=1,2,
i K i k i K i k i
kX=1 kX=1(cid:0) (cid:1)
K
. 1
s = y (x )−y¯ y (x )−y¯ , i=1, j =2,
ij i k i j k j
K
kX=1(cid:0) (cid:1)(cid:0) (cid:1)
and estimates of the quantities µ , C (0) and C (0) can simply be obtained
i ii ij
as µˆ = y¯, Cˆ (0) = s2 and Cˆ (0) = s . Then, by inverting the expressions
i i ii i ij ij
givingµ ,C (0)andC (0)intermsofparametersβ ,ς2 andς ,estimatesof
i ii ij i i ij
the latters can be obtained by βˆ =ln(µˆ )−(1/2)ςˆ2, i=1,2, where
i i i
Cˆ (0)−µˆ Cˆ (0)+(µˆ )2
ςˆ2 =ln 1+ 11 1 , ςˆ2 =ln 22 2 ,
1 (cid:18) (µˆ )2 (cid:19) 2 (cid:18)(µˆ )2 (1+1/ν)(cid:19)
1 2
and
Cˆ (0)
ij
ςˆ =ln 1+ , i=1, j =2.
ij
(cid:18) µˆ µˆ (cid:19)
i j
Let us note that estimates ςˆ2 and ςˆ do not guarantee that the symmetric
i ij
matrix Σˆ (with elements ςˆ , where ςˆ stay for ςˆ2) is positive semidefinite. A
Z ij ii i
positive semidefinite matrix Σ˜ minimizing
Z
2 2 1/2
kΣ˜ −Σˆ k= (ς˜ −ςˆ )2 ,
Z Z ij ij
(cid:18) (cid:19)
Xi=1Xj=1
can be obtained by considering Σ˜ = QΛ0QT, where Q is the 2×2 matrix
Z
whose columns are all the eingenvectors of Σˆ , and Λ0 is the 2×2 diagonal
. Z
matrixΛ0 =diag(λ ,λ ),whereλ andλ areeitherthepositiveeingenvalues
1 2 1 2
of Σˆ or zero (Rao 1973).
Z
The estimated (positive semidefinite) variance-covariance matrix of the
latent Z (x) can then be decomposed to find estimates of the parameters a
i ip
MultivariateGeostatisticalMappingofRadioactiveContamination 9
and ψ , by following standard factor analysis techniques (see, for instance,
i
Mardia, Kent and Bibby 1979). These might involve the choice of the number
P of factors to extract, as well as of which rotation to use (varimax, etc.).
The choice of the ‘shape’ parameter ν and of the autocorrelation function
ρ(h) can be based on the empirical variograms and cross-variograms γ (h)
ij
of the observable Y (x), using, for instance, the classical unbiased estimator
i
of the variogram (see, for instance, Cressie 1991). The function ρ(h) can be
selected by a direct comparison of empirical variograms γ (h) and estimated
ij
theoretical variograms γˆ (h), by a process of selection and trial. An estimate
ij
of ν can be obtained by inverting the formula giving the (theoretical) nugget
of γ (h):
22
Cˆ (0)+(µˆ )2
22 2
νˆ= −1,
γˆ (0+)
22
where γˆ (0+) is a nugget estimate.
22
3.1 Estimation with the stochastic EM algorithm
Although the above estimation procedures may easily be implemented, they
lack of sufficient optimality properties and more efficient procedures, such as
likelihood based procedures, should be considered (Minozzo 2004). Assuming
that the number of factors P, the ‘shape’ parameter ν and the spatial auto-
correlationfunctionρ(h)havealreadybeenchosen,likelihoodinferenceonthe
parameters β , ψ and a would require the maximization of the likelihood
i i ip
based on the marginal density function of the observations y (x ). However,
i k
since this marginal density is not available, and since the integration required
in the E step of the EM algorithm (Dempster, Laird and Rubin 1977) would
notbeeasy,wehavetoresorttosomestochasticversionsoftheEMalgorithm,
like the StEM algorithm or the MCEM algorithm (see, for instance, Wei and
Tanner 1990; Celeux and Diebolt 1985; Nielsen 2000; Zhang 2002; Kuhn and
Lavielle 2004).
To illustrate this estimation procedure, let us assume that m = 2 and
that Y (x) and Y (x) are conditionally distributed as Poisson and Gamma,
1 2
respectively, with h (·) = ln(·), and assume that P = 1, that is, that there is
i .
onlyonecommonfactor.Moreover,foragivensetofdatay =(y ,...,y )T,
. i i1 iK .
i=1,2, where y =y (x ), let us define the vectors of random variables ξ =
ik i k . . i
(ξ ,...,ξ )T, i = 1,2, and F = (F ,...,F )T, where ξ = ξ (x ) and
i1 . iK 1 11 1K ik . i k
F =F (x ),fori=1,2,andk =1,...,K.Then,bydefiningθ =(β,A,ψ),
1k 1 . k . .
where β = (β ,β )T, A = (a ,a )T, and ψ = (ψ ,ψ )T, the full density
1 2 11 21 1 2
function of all observed and unobserved random variables is given by
f(y ,y ,ξ ,ξ ,F )=f(y ,y |ξ ,ξ ,F )·f(ξ ,ξ ,F )
1 2 1 2 1 1 2 1 2 1 1 2 1
K K
= f (y ;M )· f (y ;M )·f(ξ )·f(ξ )·f(F),
(cid:20) 1 1k 1k(cid:21)(cid:20) 2 2k 2k(cid:21) 1 2 1
kY=1 kY=1