9.6 Precision of parameter estimates
9.6.1 Error measures
Let ˆθbe an estimator of the true parameter θ∈Rand (ˆθi)i=1,...,n be a sample of realizations of ˆθ. To quantify the estimation precision, we introduce (or recall) two error measures. The absolute error (AE) is defined as
AE(ˆθi) := AEθ(ˆθi) :=|θˆi−θ|, and theabsolute relative error (RE) is given by
RE(ˆθi) := REθ(ˆθi) := AE(ˆθi) θ . 9.6.2 Continuous presentation
In order to investigate the estimation precision for small data sets, we applied parametric bootstrap. For each parameter combination (µi, σi)i=1,...,61 estimated from the data set Schmack et al. (2015) we simulated 1000 response patterns with lengthT = 240sas in the original data. We then compared the estimators (ˆµi,σˆi)i=1,...,61 with the true parameter values using the RE. The median relative errors for the 61 parameter constellations are shown in Figure 9.12 A. Out of these, 54 (89%) showed estimation errors with median REs less than 0.25 (across the two parameters µand σ, black). The remaining simulations (dark orange) showed only few percept changes,n <20, as well as large coefficient of variation (σ/µ, Figure 9.12 B).
Simulation results for the Gamma distribution are comparable (data not shown).
Precision of parameter estimates including censored dominance times
The estimation precision of the ML estimators including the censored last dominance times can be evaluated using the same parametric bootstrap method as described above. The results are comparable, i.e., the effect of including the censored last dominance time is negligible.
This is due to the relatively large number of dominance times during continuous presentation (the median number of dominance times in Schmack et al. (2015) is 27) such that the effect of
including an additional (censored) dominance time is small.
9. A Hidden Markov Model
Figure 9.12: Precision of parameter estimates in the one-state HMM. For each of 61 parameter constellations in continuous presentation, 1000 simulations were performed with sample sizes as in Schmack et al. (2015). (A) Median of the relative error (RE) for each parameter and (B) a scatterplot of CV and n for the parameter estimates. Black points indicate constellations with mean RE across the parameters smaller than 0.25.
9.6.3 Intermittent presentation
9.6.3.1 Precision of parameter estimates using the BWA
In order to investigate the estimation precision of the Baum-Welch algorithm, we again apply parametric bootstrap to the 61 parameter combinations estimated from the response patterns to intermittent presentation in the sample data set. Figure 9.13 shows the median errors obtained in 1000 simulations for every parameter constellation forT1 = 1200 s (like in Schmack et al. (2013, 2015)) andT2= 3600 s whereπstart,S is not analyzed as it is estimated as function ofpSS andpU U (9.12). For pSS and pU U the absolute errors are presented due to the small values of the two parameters. For the time horizon of the data, T1 (panel A), 50 of the 61 parameter combinations yielded average errors (i.e., mean median errors across all parameters) smaller than 0.25 (black). The remaining cases (dark orange) showed a large CVU =σU/µU, i.e., less distinguishable stable and unstable distributions, or small sample sizen≈10 (panel C). For the larger time horizonT2 (panel B), almost all parameter combinations showed errors smaller than 0.25.
The HMM with Gamma-distributed dominance times yields comparable results concerning the precision of parameter estimation. Moreover, subjects yielding large estimation errors for the IG-model often also yield comparatively large errors for the model with Gamma-distributed dominance times (data not shown).
9.6.3.2 Precision of parameter estimates using UMVU inspired estimates We discussed in Section 9.4.1 (page 100) that UMVU inspired estimators could improve the estimation results of the estimators of the standard deviations, especially concerning bias.
Figure 9.14 now shows the simulated median bias (defined as the difference between the empirical median of the estimator and the true value) of ˆσS (panel A) and ˆσU (B) depending on the estimation procedure (ML or UMVU inspired) for the 61 subjects of Schmack et al.
(2015) and usingT1= 1200 s.
104
9. A Hidden Markov Model
Figure 9.13: Precision of parameter estimates in the two-state HMM. For each of 61 parameter constellations in intermittent presentation 1000 simulations were performed.
Log(median) of RE (for µS, σS, µU, σU) or of AE (for pSS, pU U) for T1 = 1200 (A) and T2 = 3600 (B). Black lines indicate constellations with mean error across the parameters
<0.25. (C) and (D): Scatterplot of CVU and the sample sizenwhere parameter combinations yielding large mean errors ≥0.25 are printed dark orange forT1= 1200 (C) andT2 = 3600 (D).
Figure 9.14: Comparison of ML and UMVU inspired estimators. The simulated median bias (defined as the difference between the empirical median of the estimator and the true value) is shown for ML and UMVU inspired estimators for σS (A) andσU (B). Each point represents the median bias of the ML estimator (x-axis) and the UMVU inspired estimator (y-axis) for the response pattern of one subject in Schmack et al. (2015). To derive the median bias 1000 simulations were performed. Additionally, the main diagonal is drawn.
9. A Hidden Markov Model
For almost all subjects shown in Figure 9.14 A the median bias of the UMVU inspired estimator ofσS is smaller than the median bias of the ML estimator. ForσU this effect is visible much less clearly (panel B) as there are more unstable dominance times which decrease the bias of the ML estimator and make it comparable to the UMVU inspired estimator in terms of bias. The relative errors for the two types of estimators are comparable (data not shown).
However, recall that the UMVU inspired estimators are just intuitive estimators which lack further detailed mathematical considerations.
9.6.3.3 Fitting the HMM to the data set Schmack et al. (2015) via Direct Nu-merical Maximization
We now estimate the HMM parameters of the 61 response patterns in the data set Schmack et al. (2015) via direct numerical maximization (DNM) explained in Section 9.4.2. Moreover, we contrast the results with the estimation results obtained by the Baum-Welch algorithm.
For 37 of the 61 subjects in Schmack et al. (2015) the estimated likelihoods of BWA and DNM are identical, for ten the DNM likelihood is less than one percent worse and for another ten subjects the likelihood of the DNM is even more than one percent smaller than the BWA likelihood. For four subjects the DNM approach yields slightly larger likelihoods than the BWA. In addition, simulations with the DNM estimated parameters in some cases with smaller likelihoods do not yield such convincing results as simulations with the parameters estimated using the BWA. Therefore, we recommend to apply the BWA to estimate the HMM parameters.
Precision of parameter estimates using the DNM approach
Comparing the estimation precision of the DNM approach with the BWA approach using the response patterns of the 61 subjects in Schmack et al. (2015) yields better results for the BWA. For the recording lengthT = 1200 s the DNM yields 27 subjects with a mean median relative error larger than 0.25, where the BWA only yields eleven subjects. ForT = 3600 s there are for the BWA three subjects with a mean median relative error larger than 0.25 and for the DNM twelve subjects.
Precision of parameter estimates including censored dominance times
The estimation precision when including the censored dominance times via the DNM method can be evaluated using the data set of (Schmack et al., 2015). It should be compared to the DNM method in the paragraph above. The results are comparable, i.e., the effect of including the censored last dominance time is small.
106
Chapter 10
The HMM: Theoretical properties
In this chapter we aim to investigate the mathematical properties of the sequences of dominance times generated by the Hidden Markov Model more in detail thereby also deriving quantities important for the comparison of clinical groups. There are several possibilities to interpret these sequences. First, the points in time where the perception changes can be understood as point process on the (non-negative) real line. Second, the perceptual reversals may be interpreted as renewal points and the whole process in the case of continuous stimulation as renewal process. When incorporating the stable and unstable phases occurring in intermittent presentation, the process can be connected to an alternating renewal process (e.g., Medhi, 2009) introduced in Chapter 8.2.3.1.
In this chapter we always assume HMMs with inverse Gaussian distributed emissions. In Section 10.1 we derive results on the distribution and expectation of the number of perceptual reversals during continuous presentation. The central theme of Section 10.2 is the theoretical investigation of the point process induced by the HMM for intermittent presentation (HMMi).
First, we explain the connection to semi-Markov processes. Then, the number of changes is discussed as well as first passage times, stationary distributions and renewal equations. Note that first passage times are required for the derivation of the steady-state distribution, which is very important to analyze differences between clinical groups as an increased time spent in the unstable state is an indicator for a less stable perception. Thus, Corollary 10.8 about the steady-state distribution is the result in this chapter being most important for application.
Moreover, we investigate the expectation and distribution of the residual time. A knowledge about the residual time enables us to make a prognosis about the next perceptual reversal.
We understand ΞHMMc and ΞHMMi as the point processes on the non-negative line generated by the points in time (t0, t1, t2, . . . , tn) of the perceptual reversals of the HMM for continuous (HMMc) and for intermittent presentation, respectively. There, we uset0= 0. Formally, for
the set of realized dominance times (d1, d2, . . . , dn) of a HMMc ΞHMMc :={0} ∪ {t∈R|
k
X
i=1
di =t, k= 1,2, . . . , n}={t0, t1, t2, . . . , tn} and equivalently for the set of dominance times (d1, d2, . . . , dn) of a HMMi
ΞHMMi:={0} ∪ {t∈R|
k
X
i=1
di=t, k= 1,2, . . . , n}={t0, t1, t2, . . . , tn}.
10. The HMM: Theoretical properties
We define ˜Yt:= ( ˜Yt)t≥0 as the hidden state at time tgoverning the point process ΞHMMi (in contrast toYi which is the hidden state of thei-th dominance time). In precise terms,
Y˜t:=j∈ {S, U}|
k
X
i=1
di ≤t <
k+1
X
i=1
di, Yk+1 =j, k= 0,1,2, . . . , n−1
=Yi on ti−1 ≤t < ti, withP0
i=1di := 0. The difference betweenYiand ˜Ytis illustrated in Figure 10.1 for a simulated example of a HMMi. Panel C shows ˜Yt as variable defined on a continuous time space, and panel D showsYi as variable defined on a discrete time space. ˜Yt may be described using an alternating renewal process with states S and U or a regenerative process (compare Sections 8.2.3.1 and 8.2.3.2).
0 T
L R
perc.
t
A
0 T
U S
Y~ t
t
C
●●
●●●●●●●●●
●
0 4 8 12
5 180di
i
B
●●
●●●●●●●●●
●
0 4 8 12
U S Yi
i
D
Figure 10.1: Overview of different variables in the HMM. (A) The perception at time t∈[0, T]. (B) The resulting dominance times (di)i≥1. (C) The hidden state Y˜t on the continuous axis [0, T]. (D) The hidden state Yi of each dominance time di.
Note that ΞHMMi is not a continuous time Markov chain as the probabilities
P( ˜Yt=S|Y˜s=S),0≤s < tdepend on how long the current state isS and thus also depend on the state of the chain at points in times1, s2, . . . < s. As a simple example, imagine that the hidden state changes at time t = 100 from unstable to stable. Then, we observe that P( ˜Y200 =S|Y˜150 =S)6=P( ˜Y200=S|Y˜150=S,Y˜101 =S) as the inverse Gaussian distribution is not memoryless.
Often we speak of ”phases”: A stable phase comprises all stable dominance times between a change from the unstable to the stable state and the change back to the unstable state. An unstable phase is defined analogously.
10.1 Continuous presentation
In this section the number of changes in the HMM as well as the residual time are analyzed.