# Fiducial distribution and sequential monte-carlo algorithm

I have spent a little time reading "Generalized fiducial inference for normal mixed models" by Cisewski and Hannig. First of all I am interested in understanding how to simulate the fiducial distribution of the parameters. But I don't understand this point.

Consider for instance the one-way random effects model $y_{ij}=\mu + \alpha_i + \epsilon_{ij}$ with $\alpha_i \sim {\cal N}(0,\sigma_a)$ independent of $\epsilon_{ij} \sim {\cal N}(0,\sigma)$.

• First of all the authors write $\alpha_i= \sigma_a u_{i}$ and $\epsilon_{ij}=\sigma v_{ij}$ with $u_{i}$ and $v_{ij}$ independent standard normal random variables.

• In fact they consider the data are rounded. For instance when $y_{ij}=1.32$ then actually the $(i,j)$-datum is the interval $[a_{ij}, b_{ij}[=[1.315, 1.325[$ because the "real" observation is rounded with two decimals. Thus the model is $$a_{ij} \leq \mu + \sigma_a u_i+\sigma v_{ij} < b_{ij}.$$.

• Denote by $Z$ the random vector made of all the random components of the model: all the $u_i$'s and all the $v_{ij}'s$.

• Define the set $Q\bigl(a,b, u, v\bigr) = \bigl\{ (\mu,\sigma_a,\sigma) \mid a \leq \mu + \sigma_a u + \sigma v < b\bigr\}$ and then define the random set $$Q(\text{data}, Z) = \bigcap_{i,j}Q\bigl(a_{ij},b_{ij},u_i,v_{ij}\bigr)$$ (with randomness contained in $Z$, and the data are fixed). When this set is nonempty we denote by $V\left(Z\right)$ an element of this set selected according to a specified rule (at least in a measurable way). The (generalized) fiducial distribution of $(\mu,\sigma_a,\sigma_b)$ is the conditional distribution of $V(Z)$ given $Q(\text{data}, Z) \neq \varnothing$.

• Then the authors propose a sequential Monte-Carlo (SMC) generator of the fiducial distribution, and this is the point that confuses me. It starts on page 9 of the paper.

• The authors claim that generating the fiducial distribution is equivalent to generating $Z$ conditional on $Q(\text{data}, Z) \neq \varnothing$. I agree with this part but they do not propose a measurable rule $V$, or I have missed this point while reading the paper (or is it an obvious point ?).

• For $t=1, \ldots, n$ where $n$ is the sample size, denote by $Z_{1:t}$ the random components corresponding to the first $t$ data values (my notation is not appropriate here). Then the authors propose a sequential algorithm whose first step at time $t$ requires them to simulate the distribution $\pi_t$ of $Z_{1:t}$.

• There are some undefined notations in the preprinted version whose link is given above: $Q_t$ and $Q_t^{(J)}$ and the definition of the set $C_t$ is not well-written. It is easy to guess the definition of $Q_t$: for $t=1,\ldots,n$ the random set $Q_t$ is constructed similarly to $Q(\text{data}, Z)$ but only with the first $t$ available data points and their associated random components in the model (my notation is not appropriate for defining $Q_t$), and $C_t$ are the possible values of the random components for which $Q_t$ is not empty.

• Then $\pi_t$ is the independent product of the standard normal distributions truncated to $C_t$. I am totally lost here. In particular, does the algorithm require one to derive $C_t$ explicitly ?

I would be glad if you someone could clear up the points that I find confusing, or give me a link to a similar algorithm which is more detailed.

-
Stephane I have edited your question to improve the English. i hope you don't mind. I think your English is good and the meaning of what you write is mostly very clear. I did the for two reasons: (a) I am a compulsive editor and (b) I think that maybe these edits can help you to improve your English a little. If you don't like it I won't do it again. Please just in case I misrepresented something that you were trying to say. –  Michael Chernick Sep 20 '12 at 21:30
+1 Very good question by the way! –  Michael Chernick Sep 20 '12 at 21:33
@MichaelChernick Thanks Michael. Of course I trust you. –  Stéphane Laurent Sep 21 '12 at 5:09
I am glad and I hope this is helpful to you. If it is it is worth it for me to make the effort. –  Michael Chernick Sep 21 '12 at 10:20
@MichaelChernick Yes this is very helpful to reread my English after it has been improved :) –  Stéphane Laurent Sep 21 '12 at 11:58

The whole goal was to generate standard normal conditional on the $C$ (using all of the data). That proved to be a rather difficult problem and therefore we have chosen to do this sequentially. You have correctly understood the meaning of $C_t$ and $Q_t$.

By the way, $C_t$ is formally defined in the middle of the page 9. You are correct the definition of $Q_t$ has slipped through the cracks of several revisions but you did guess its meaning correctly. Thank you for pointing this out.

The main idea is that it is easy to move from $C_t$ to $C_{t+1}$—that is all the business with $m_t$ and $M_t$ that follows. The reason why we choose sequential MC is that we can easily generate from the conditional distribution $Z_{t+1} I_{C_{t+1}}$ given $Z_t$ [using $m_t$ and $M_t$]. However there is a need to re-weight due to the fact that marginal distribution of $Z_{t+1} I_{C_{t+1}}$ marginalized to the time $t$ is different than the distribution of $Z_t I_{C_t}$.

The details of how much to re-weight is on the page 10 with more details in the Appendix C.

-
(+1) Welcome to the site, Jan. This site supports $\LaTeX$; I've tried to typeset the math. Please feel free to make edits if I've introduced errors. –  cardinal Sep 21 '12 at 0:25
@JanHannig Out of curiousity, did you come to the site because someone alerted you about Stephane's question. In any case we are privileged to have you come and answer. We welcome people with your expertise to register on the site and visit when you can to contribute. –  Michael Chernick Sep 21 '12 at 2:02
Thank you Jan, I will carefully read your answer today. Still I'm not sure $C_t$ is well defined: the definition you wrote depends on $\beta_i$ and $\sigma$ which are not quantified. Michael, that's me who sent an e-mail to Pr Hannig. –  Stéphane Laurent Sep 21 '12 at 5:14
Ok now I have read your answer. Now the main idea is clearer but still I do not understand the details (and I don'y know what are $J$ and $Q_t^{J}$. As you suggested me in your email I will ask Prof. Cisewski to send me her thesis. You also pointed your Matlab code out, I will study it unc.edu/~hannig/r_download.html I'm not a Matlab user but some colleagues could help. –  Stéphane Laurent Sep 21 '12 at 7:50

Thank you for your questions. I just want to elaborate on a couple points you made above:

1. I just want to clarify your description of the interval, or "fat", data: each observation is considered an interval rather than an exact value, but the intervals are fixed and defined based on the data collection/data storing process. The clarification is for your example where you say if your observation is 1.32, then the interval would be [1.315, 1.325). In this illustration, you are making the observation the midpoint of the interval and allowing more significant digits than the original datum, which is not what we use. Let's suppose your observation 1.32 is really 1.32 meters, and your measuring devise is a ruler than has tick marks down to millimeters. Then the observation is recorded as 1.32 meters, but really, it could be 1.312 m or 1.318 m (this is supposing that the practice is to round to the next highest mm). Hence, the observation "1.32" would be placed in the interval [1.31, 1.32). The point is that the intervals are not arbitrarily defined around each point, but are set in advance as a fixed grid (in this case, the fixed grid is determined by the ruler).

2. I also want to note that there is no assumption about the distribution across an interval. This leads into your other question about V(). We maintain the intervals throughout the algorithm, and at t = n we have a sample of weighted particles that, geometrically, are polyhedrons on the parameter space. One can then select a value within each weighted polyhedron according to some rule, V(). In the simulation study in our paper, we randomly select either the upper or lower endpoint of the marginal intervals for each dimension of the parameter space. There are other options here, and Hannig (2009) "On Generalized Fiducial Inference" describes some alternatives.

-
Thanks for these clarifications! About the first point, I was just supposing the practice is to round to the closest cm. About the 2nd point I still need to read your dissertation to understand the algorithm and this notion of "weighted particles" which I have never encountered before. –  Stéphane Laurent Sep 23 '12 at 16:23