---
title: Hidden Markov Models and State Estimation
author: 36-467/36-667
date: 19 November 2020 (Lecture 23)
output: slidy_presentation
bibliography: locusts.bib
---

\[
\newcommand{\Prob}[1]{\mathbb{P}\left( #1 \right)}
\]

```{r, include=FALSE}
set.seed(2018-12-05)
library(knitr)
opts_chunk$set(size="small", background="white",
               cache=TRUE, autodep=TRUE,
               tidy=TRUE, tidy.opts=list(comment=FALSE),
               warning=FALSE, message=FALSE,
               echo=FALSE)
```


## Housekeeping

- Class is canceled for Tuesday
    + I'll post notes for what would have been the optional lecture
- Homework 12 will be due on 3 December (as scheduled)


## Previously

- Markov processes are ones where conditioning on the whole past is equivalent to conditioning on the current observation:
\[
\Prob{X(t+1)|X(1:t)} = \Prob{X(t+1)|X(t)}
\]
- This makes $X(t)$ the **state** of the process at time $t$, the variable which fixes the distribution of future events
- The Markov property makes it easy to:
    + Do probability calculations
	+ Do prediction
	+ Do inference
- The Markov property commits us to $X(t+1) \perp X(1:t-1) | X(t)$

## Markov Property vs. Measurement Noise

- What if we don't observe the state perfectly?
- $S(t)$ is a Markov process
- We only observe $X(t) = S(t) + \epsilon(t)$, where $\epsilon(t)$ is IID noise
- Will $X(t)$ also be Markov?


## A Concrete Example

- $S(t)$ is a Markov chain with states $-1$ and $+1$
- $\Prob{S(t+1)=+1|S(t)=+1} = \Prob{S(t+1)=-1|S(t)=-1} \equiv q > 0.5$
- We observe not $S(t)$ but $X(t) = S(t) + \epsilon(t)$
- Noise $\epsilon(t) \sim \mathcal{N}(0,1)$ and is IID
- Is this $X(t)$ a Markov process?


```{r}
# From lecture 18: generic, abstract Markov simulator
rmarkov <- function(n, rinitial, rtransition) {
  x <- vector(length=n)
  x[1] <- rinitial()
  for (t in 2:n) {
    x[t] <- rtransition(x[t-1])
  }
  return(x)
}

# Simulate a chain that's +1/-1, with persistence in state
  # Pr(X(t+1)=+1|X(t)=+1) = Pr(X(t+1)=-1|X(t)=-1) = q, assumed > 0.5
  # Exercise: Prove that the invariant distribution is (1/2, 1/2) for all q
rpmchain <- function(n, q) {
    # Take one step of the chain
    transition.with.persistence <- function (x) {
        # Stay in the same state with probability q
        stay <- sample(c(TRUE, FALSE), size=1, prob=c(q, 1-q))
        if (stay) { return(x) }
        # Or flip to the opposite state
        if (!stay) { return(-x) }
    }        
    pm.seq <- rmarkov(n, rinitial=function() { sample(c(-1, +1),
                                                          size=1) },
                      rtransition=transition.with.persistence)
    return(pm.seq)
}

# Simulate a +1/-1 chain, but observe it through noise
noisy.rpmchain <- function(n, q, sd=1) {
    pm.seq <- rpmchain(n, q)
    x <- rnorm(n=n, mean=pm.seq, sd=sd)
    return(list(x=x, s=pm.seq))
}
```

```{r}
sim <- noisy.rpmchain(n=100, q=0.75)
plot(sim$x, xlab="t", ylab="X(t)", xlim=c(0, 120))
lines(1:100, sim$s, col="grey")
legend("right", legend=c("Hidden state S(t)", "Observable X(t)"),
       pch=c(NA, 1), lty=c("solid", "blank"), col=c("grey", "black"),
       cex=0.5)
```

## The Current State "Screens Off" the Past from the Future

$\Prob{X(t+1)=x(t+1)|S(t)=s, X(1:t)=x(1:t)}$
\begin{eqnarray}
& = & \sum_{r}{\Prob{X(t+1)=x(t+1)|S(t+1)=r, S(t)=s, X(1:t)=x(1:t)}\Prob{S(t+1)=r|S(t)=s,X(1:t)=x(1:t)}}\\
& = & \sum_{r}{\Prob{X(t+1)|S(t+1)=r}\Prob{S(t+1)=r|S(t)=s}}
\end{eqnarray}

## Relating $X(t)$ to Its Past through the State

$\Prob{X(t+1)=x(t+1)|X(1:t)=x(1:t)}$
\begin{eqnarray}
& = & \sum_{s}{\Prob{X(t+1)=x(t+1)|S(t)=s,X(1:t)=x(1:t)}\Prob{S(t)=s|X(1:t)=x(1:t)}}\\
& = & \sum_{s}{\left(\sum_{r}{\Prob{X(t+1)|S(t+1)=r}\Prob{S(t+1)=r|S(t)=s}}\right) \Prob{S(t)=s|X(1:t)=x(1:t)}}
\end{eqnarray}

$\Rightarrow$ Given $X(t)$, the only way $X(t-1)$ can matter for $X(t+1)$ is by changing the conditional distribution of $S(t)$


## Relating $X(t)$ to $S(t)$

- Any value of $X(t)$ _could_ have come from _either_ value of $S(t)$
```{r}
curve(0.5*dnorm(x,-1,1)+0.5*dnorm(x,1,1), xlab="x", ylab="p(x)",
      main="Marginal and Conditional Distributions of X(t)", xlim=c(-5,5), ylim=c(0, 0.5), col="purple", lty="dotted")
curve(dnorm(x,-1,1), add=TRUE, col="blue")
curve(dnorm(x,+1,1), add=TRUE, col="red")
legend("topright", legend=c("X(t)", "X(t)|S(t)=-1", "X(t)|S(t)=+1"),
       lty=c("dotted", "solid", "solid"), col=c("purple", "blue", "red"))
```
- $X(t) \gg 0$ is _unlikely_ if $S(t)=-1$ but not impossible
- $X(t) \approx 0$ is about equally like for either value of $S(t)$


## History Matters for $X(t)$

- If $X(t) = `r signif(max(sim$x), 2)`$ then it's really likely that $S(t)=+1$ and conditioning on earlier $X$s can't really change that
- Earlier values of $X$ give information about $S(t)$ when $X(t)$ isn't too extreme, _and so_ they give information about $X(t+1)$
- In a Markov process, $X(t-1) \perp X(t+1)|X(t)$
- That isn't true here
- $\therefore$ the observed process isn't Markov



## Hidden Markov Models / State-Space Models

- **Hidden** or **latent** (= unobserved) **state** $S(t)$, which follows a Markov process
- Transition probabilities:
\[
q(s,r) \equiv \Prob{S(t+1)=r|S(t)=s}
\]
    + Need a distribution for initial state as well
	+ A.k.a. transition model, process model, state model, plant model
- **Observable** or **manifest** $X(t)$ is a noisy function of the current state:
\[
X(t) = f(S(t), \epsilon(t))
\]
for some function $f$ and IID noise $\epsilon$
    + Function $f$ needn't be additive and might be very nonlinear
- Leads to **observation probabilities** or **observation model**
\[
g(s,x) \equiv \Prob{X(t)=x|S(t)=s}
\]

(for continuous variables, replace PMFs by pdfs as needed)



## Some Basic Calculations for HMMs {.smaller}

- Probability of a state sequence from Markov property:
\begin{eqnarray}
\Prob{S(1:n) = s(1:n)} & = & \Prob{S(1)=s(1)} \prod_{t=2}^{n}{S(t)=s(t)|S(1:t-1)=s(1:t-1)}\\
& = &  \Prob{S(1)=s(1)} \prod_{t=2}^{n}{S(t)=s(t)|S(t-1)=s(t-1)}\\
& = & \Prob{S(1)=s(1)} \prod_{t=2}^{n}{q(s(t-1), s(t))}
\end{eqnarray}
- Probability of an observed sequence, given a state sequence:
\begin{eqnarray}
\Prob{X(1:n)=x(1:n)|S(1:n)=s(1:n)} & = & \prod_{t=1}^{n}{\Prob{X(t)=x(t)|S(t)=s(t)}}\\
& = & \prod_{t=1}^{n}{g(s(t), x(t))}
\end{eqnarray}
- Joint probability of states and observations:
\begin{eqnarray}
\Prob{X(1:n)=x(1:n), S(1:n)=s(1:n)} & = & \Prob{S(1:n)=s(1:n)}\Prob{X(1:n)=x(1:n)|S(1:n)=s(1:n)}\\
& = & \Prob{S(1)=s(1)} \prod_{t=2}^{n}{q(s(t-1), s(t))}\prod_{t=1}^{n}{g(s(t), x(t))}
\end{eqnarray}
- Probability of an observation sequence:
\begin{eqnarray}
\Prob{X(1:n)=x(1:n)} & = & \sum_{s(1:n)}{\Prob{X(1:n)=x(1:n), S(1:n)=s(1:n)}}\\
& = & \sum_{s(1:n)}{ \Prob{S(1)=s(1)} \prod_{t=2}^{n}{q(s(t-1), s(t))}\prod_{t=1}^{n}{g(s(t), x(t))}}
\end{eqnarray}


## Why We Want the Probability of an Observation Sequence

\[
\Prob{X(1:n)=x(1:n)}  =  \sum_{s(1:n)}{ \Prob{S(1)=s(1)} \prod_{t=2}^{n}{q(s(t-1), s(t))}\prod_{t=1}^{n}{g(s(t), x(t))}}
\]

- When we predict, we want
\[
\Prob{X(t)=x|X(1:t-1)=x(1:t-1)} = \frac{\Prob{X(1:t)=x(1:t)}}{\Prob{X(1:t-1)=x(1:t-1)}}
\]
- $\therefore$ Predicting comes down to find the probabilities of sequences and dividing
- $\Prob{X(1:n)=x(1:n);\theta}=$ the likelihood of the data with parameter value $\theta$ so we want this for inference as well
- But our equation is a sum over state histories, and that's bad
    + Even 2 states means $2^n$ possible values of $s(1:n)$
	+ Never mind continuous states


## A Recursion Trick (a.k.a. "The Forward Algorithm") [1] {.smaller}

- If we knew $S(1)$, it'd be easy to predict $X(1)$:
\[
\Prob{X(1)=x|S(1)=s} = g(s, x)
\]
- Leads to the _marginal_ distribution of $X(1)$, what we should predict for the first observation:
\[
\Prob{X(1)=x} = \sum_{s}{\Prob{X(1)=x|S(1)=s}\Prob{S(1)=s}} = \sum_{s}{g(s, x) \Prob{S(1)=s}}
\]
- Now estimate the state at time 1 _after_ seeing that $X(1)$ is actually $x(1)$:
\begin{eqnarray}
\Prob{S(1)=s|X(1)=x(1)} & = & \frac{\Prob{S(1)=s, X(1)=x(1)}}{\Prob{X(1)=x(1)}}\\
& = & \frac{\Prob{S(1)=s}\Prob{X(1)=x(1)|S(1)=s}}{\Prob{X(1)=x(1)}}\\
& = & \frac{\Prob{S(1)=s} g(s, x(1))}{\sum_{s^{\prime}}{\Prob{S(1)=s^{\prime}} g(s^{\prime}, x(1))}} \equiv F_1(s)
\end{eqnarray}

## A Recursion Trick (a.k.a. "The Forward Algorithm") [2] {.smaller}

- We have the conditional distribution of $S(1)$, extrapolate it forward:
\begin{eqnarray}
\Prob{S(2)=r|X(1)=x(1)} & = & \sum_{s}{\Prob{S(2)=r, S(1)=s|X(1)=x(1)}}\\
& = & \sum_{s}{\Prob{S(2)=r|S(1)=s, X(1)=x(1)}\Prob{S(1)=s|X(1)=x(1)}}\\
& = & \sum_{s}{\Prob{S(2)=r|S(1)=s}\Prob{S(1)=s|X(1)=x(1)}}\\
& = & \sum_{s}{q(s,r) F_1(s)}
\end{eqnarray}
- The extrapolated state lets us predict $X(2)$:
\begin{eqnarray}
\Prob{X(2)=x|X(1)=x(1)} & = & \sum_{s}{\Prob{X(2)=x, S(2)=s|X(1)=x(1)}}\\
& = & \sum_{s}{\Prob{X(2)=x|S(2)=s, X(1)=x(1)}\Prob{S(2)=s|X(1)=x(1)}}\\
& = & \sum_{s}{\Prob{X(2)=x|S(2)=s}\Prob{S(2)=s|X(1)=x(1)}}\\
& = & \sum_{s}{g(s, x)\Prob{S(2)=s|X(1)=x(1)}}
\end{eqnarray}
- After we see $X(2)$, use both it _and_ $X(1)$ to update $S(2)$:
\begin{eqnarray}
F_2(s) & \equiv & \Prob{S(2)=s|X(2)=x(2), X(1)=x(1)}\\
& = & \frac{\Prob{S(2)=s, X(2)=x(2)|X(1)=x(1)}}{\Prob{X(2)=x(2)|X(1)=x(1)}}\\
& = & \frac{\Prob{X(2)=x(2)|S(2)=s, X(1)=x(1)}\Prob{S(2)=s|X(1)=x(1)}}{\Prob{X(2)=x(2)|X(1)=x(1)}}\\
& = & \frac{\Prob{X(2)=x(2)|S(2)=s}\Prob{S(2)=s|X(1)=x(1)}}{\Prob{X(2)=x(2)|X(1)=x(1)}}\\
& = & \frac{g(s, x(2))\Prob{S(2)=s|X(1)=x(1)}}{\Prob{X(2)=x(2)|X(1)=x(1)}}
\end{eqnarray}

## A Recursion Trick (a.k.a. "The Forward Algorithm") [3] {.smaller}


#### The Forward Algorithm

- Begin with $\Prob{S(t)=s|X(1:t)=x(1:t)} \equiv F_t(s)$
- Extrapolate the state forward:
\[
\Prob{S(t+1)=s|X(1:t)=x(1:t)} = \sum_{s^{\prime}}{q(s^{\prime}, s)\Prob{S(t)=s^{\prime}|X(1:t)=x(1:t)}} = \sum_{s^{\prime}}{q(s^{\prime}, s) F_t(s^{\prime})}
\]
- Predict the next observation:
\[
\Prob{X(t+1)=x|X(1:t)=x(1:t)} = \sum_{s}{g(s,x) \Prob{S(t+1)=s|X(1:t)=x(1:t)}}
\]
- Condition $S(t+1)$ on $X(t+1)$:
\[
\Prob{S(t+1)=s|X(1:t+1)=x(1:t+1)} = \frac{g(s,x(t+1)) \Prob{S(t+1)=s|X(1:t)=x(1:t)}}{\Prob{X(t+1)=x|X(1:t)=x(1:t)}}
\]
- End with $F_{t+1}(s) \equiv \Prob{S(t+1)=s|X(1:t+1)=x(1:t+1)}$, or begin the cycle over again with $t+2$.

## Back to the Likelihood

- We've gotten predictions along the way with the forward algorithm
- The likelihood comes by multiplying:
\[
\Prob{X(1:n)=x(1:n)} = \Prob{X(1)=x(1)}\prod_{t=2}^{n}{\Prob{X(t)=x(t)|X(1:t-1)=x(1:t-1)}}
\]

## Actually Doing the Calculations

- We can do everything exactly when $S(t)$ and $X(t)$ are both discrete
    + The original "forward algorithm"
- We can do everything exactly when $S(t)$ and $X(t)$ are both Gaussian, and everything's linear
    + "Kalman filter"
- Otherwise we need numerical approximations


## The Particle Filter: Words {.smaller}

- The key thing is to get $F_t(s) = \Prob{S(t)=s|X(1:t)=x(1:t)}$
- In the **Monte Carlo filter**, or **particle filter**, we approximate $F_t(s)$ by a large number of random samples from the right distribution

#### The Particle Filter

1. At $t=1$, make $m$ independent draws  $S^*_1, S^*_2, \ldots S^*_m$ from $\Prob{S(1)}$
    + These are the **particles**
2. For each $S^*_i$, calculate $g(S^*_i, x(t)) \equiv w_{i}$
    + $\Prob{X(t)=x(t)|X(1:t-1)=x(1:t-1)} \approx m^{-1} \sum_{i=1}^{m}{w_i}$
3. Resample the particles: draw $m$ of them, with replacement, from $S^*$, with probabilities $\propto w_i$; call these $\tilde{S}_1, \tilde{S}_2, \ldots \tilde{S}_m$
4. Set $\hat{F}_t$ to the sample distribution of the $\tilde{S}$
    + The distribution of the resampled particles approximates $\Prob{S(t)|X(1:t)=x(1:t)}$
5. Increment $t$ by 1, and **evolve** the particles forward in time: for each $\tilde{S}_i$, draw $S^*_i$ from $q(\tilde{S}_i, S^*_i)$
    + Distribution of the new $S^*$ approximates $\Prob{S(t+1)|X(1:t)=x(1:t)}$.
6. Go to (2)


## The Particle Filter: Code

```{r, echo=TRUE, tidy=TRUE, tidy.opts=list(comment=FALSE)}
# Particle filter for an HMM
# Inputs: time series of length n (x)
  # function generating draws for initial state (rinitial)
  # function generating draws for next state (rtransition)
  # function giving observational probabilities (dobs)
  # number of particles (m)
# Output: n*m array of particles, approximating Pr(S(t)|X(1:t))
# Presumes: states can be elements in a vector
  # observations can be elements in a vector
  # rinitial is a no-argument function (for compatibility with rmarkov())
  # rtransition is a one-argument function (ditto)
  # dobs is a two-argument function
particle.filter <- function(x, rinitial, rtransition, dobs, m) {
    # How long is the time series?
    n <- length(x)
    # Our initial guesses at the state, before seeing anything
    particles <- replicate(m, rinitial())
    # A storage area for our guesses
    particle.history <- matrix(0, nrow=m, ncol=n)
    for (t in 1:n) {
        # Calculate observationa likelihoods/weights for current observation
        # given state=particles
        weights <- dobs(state=particles, observation=x[t])
        # Resample the particles, with probabilities (proportional to)
        # the weights
        particles <- sample(particles, size=m, replace=TRUE, prob=weights)
        # Store the current set of guesses
        particle.history[,t] <- particles
        # Evolve the particles forward using the Markov transitions for the
        # next time step
        particles <-sapply(particles, rtransition)
    }
    return(particle.history)
}
```


## The Particle Filter: For Our Little Demo

```{r, echo=TRUE, tidy=FALSE}
sim.pf <- particle.filter(x=sim$x,
                          rinitial=function() { sample(-1,+1, size=1) },
                          rtransition=function (x) {
                              stay <- sample(c(TRUE, FALSE), size=1, prob=c(0.75, 0.25))
                              if (stay) { return(x) }
                              if (!stay) { return(-x) }
                          },
                          dobs=function(state, observation) {
                              dnorm(x=observation, mean=state, sd=1)
                          },
                          m=100)
```



## The Particle Filter: Run on Our Little Demo

```{r}
plot(sim$x, xlab="t", ylab="X(t)", xlim=c(0, 120))
lines(1:100, sim$s, col="grey", lwd=2)
pf.limits <- apply(sim.pf, 2, quantile, prob=c(0.05, 0.95))
pf.means <- apply(sim.pf, 2, mean)
lines(1:100, pf.limits[1,], col="blue", lty="dashed")
lines(1:100, pf.limits[2,], col="blue", lty="dashed")
lines(1:100, pf.means, col="blue")
legend("bottomright", legend=c("Hidden state S(t)", "Observable X(t)", "Particle filter (mean)", "Particle filter (90% Interval)"),
       pch=c(NA, 1, NA), lty=c("solid", "blank", "solid"), "dashed", col=c("grey", "black", "blue", "blue"),
       cex=0.5)
```


## Issues with the Particle Filter

1. Monte Carlo error: really want $m$ to be as large as possible
    + But it costs computing time
	+ Error is $O(1/\sqrt{m})$ so diminishing returns
2. "Particle collapse": if one $w_i \gg w_j$ for all the other $j$, then (almost) all new particles are copies of $S^*_i$ and the approximation is bad
    + Lots of hacks around this
3. What if we don't know transitions $q$ or observation probabilities $g$?



## Parameter Estimation

- Be explicit about parameters in the notation: it's $q(s,s^{\prime};\theta)$ and $g(s,x;\theta)$; how do we find $\theta$?
- We could use the likelihood: use forward algorithm to get $\prod_{t}{\Prob{X(t)=x(t)|X(1:t-1) = x(1:t-1);\theta}}$, and maximize
    + Or, as usual, maximize the log
	+ It's a complicated function of $\theta$ but it is just a function of $\theta$
- Alternate between estimate states and estimating parameters
    + Start with a guess about $q$ and $g$, estimate $S(1:t)$
	+ Using $S(1:t)$, estimate $q$ and $g$
	+ Iterate to convergence
- Calculating the likelihood is still computationally demanding, but simulating is easy $\Rightarrow$ Use all our simulation-based-inference tricks

## What About Space?  What About Time?

- Hidden Markov random field: $S(r)$ is a Markov random field, $X(r) = f(S(r), \epsilon(r))$
    + $X(r)$ isn't Markov
	+ Simulation: Gibbs for the $S$ field, then draw from $X(r)|S(r)$ independently
	+ Likelihood: There are recursive approaches, though rather more complicated than the forward algorithm
- Spatio-temporal hidden Markov random fields: $S(r,t)$ is a spatio-temporal random field, $X(r,t) = f(S(r,t), \epsilon(r,t))$
- OR: dynamic factor models
    + $S(t)$ is a low-dimensional vector
	+ $S(t)$ evolves according to a Markov process (often assumed linear)
	+ $X(r,t) = f_r(S(t), \epsilon(r,t))$, so $X(r,t) \perp X(q,t) | S(t)$
	    * This measurement model is also often assumed to be linear


## Summing Up

- If we observe a Markov process with noise, the observables aren't Markov
- We often _postulate_ such hidden Markov models, because they let us fit non-Markov processes but still have some of the nice structure of Markov models
- Estimating the distribution of the hidden state $S(t)$ given the history of observations $X(1:t)$ is crucial
    + The forward algorithm lets us do this recursively
	+ The particle filter lets us do this recursively by simulation
- Inference then proceeds either by likelihood (found using the forward algorithm) or maybe by simulation-based inference
- The spatial and spatio-temporal variants are easy to formulate, but harder to work with

Many more details in the handout on the website


## Backup: Applying Indirect Inference to Our Demo

-  We need:
    1. Simulator for the generative model
	2. Estimator for the auxiliary model
	3. Code to optimize the generative model so auxiliaries fit to simulations match the auxiliary fit to the data
- We'll use AR(2) as the auxiliary
    + One more auxiliary parameter than we really need to estimate $q = \Prob{S(t+1)=+1|S(t)=+1} = \Prob{S(t+1)=-1|S(t)=-1}$


## Backup: Finding the Auxiliary Parameter Estimates

```{r, echo=TRUE}
# Estimate the auxiliary model (here, an AR(2)) for indirect inference
# Inputs: a time series (x)
# Output: Point estimate of the auxiiary model (here, a vector of length 2)
aux.estimator <- function(x) {
    ar2.fit <- ar.ols(x, aic=FALSE, order=2)
    return(ar2.fit$ar)
}

# Repeatedly simulate the generative model, fit the auxiliary, and average
# Inputs: parameter of generative model (q)
  # number of simulations to run (s)
  # length of each simulation (n)
# Output: vector of averaged auxiliary parameter estimates
aux.from.sim <- function(q, s, n=length(sim$x)) {
    estimates <- replicate(s, aux.estimator(noisy.rpmchain(n=n, q=q)$x))
    return(rowMeans(estimates))
}    
```


## Backup: Can We Go From Auxiliary Parameters to Generative Parameters?

```{r}
plottable.aux.from.sim <- Vectorize(aux.from.sim, vectorize.args="q")
ar2.vs.q <- sapply(seq(from=0.01, to=0.99, by=0.01),
                   plottable.aux.from.sim, s=500)
```

```{r}
plot(ar2.vs.q[1,], ar2.vs.q[2,], type="l", xlab="Coefficient on X(t-1)",
     ylab="Coefficient on X(t-2)")
text(ar2.vs.q[1,1], ar2.vs.q[2,1], labels="q=0.01")
text(ar2.vs.q[1,99], ar2.vs.q[2,99], labels="q=0.99")
```



## Backup: Optimizing the Generative Model so Simulation-Auxiliaries Match Data-Auxiliaries


```{r, echo=TRUE}
toy.indirect.inference <- function(x, starting.q, s) {
    n <- length(x)
    # What's the auxiliary model on the actual data?
    aux.data <- aux.estimator(x)
    # Calculate the discrepancy between the auxiliary on the data and
    # the auxiliary we get from simulating at a particular parameter value
    discrepency <- function(param) {
        q <- param
        # Get the average of auxiliary estimates over s simulation runs
        aux.mean <- aux.from.sim(q=q, s=s, n=n)
        # Mean squared difference over auxiliary parameters
        return(sum((aux.data-aux.mean)^2))
    }
    # Now minimize the discrepancy
    fit <- optim(par=starting.q, fn=discrepency, method="Brent", lower=0,
                 upper=1)
    # Return just the best-fitting value
    return(fit$par)
}    
```

- Run:

```{r, echo=TRUE}
(toy.ii <- toy.indirect.inference(x=sim$x, starting.q=0.5, s=500))
```
- Remember we set the true $q$ to $0.75$

- This isn't a fluke:
```{r, echo=TRUE}
summary(replicate(30, toy.indirect.inference(x=sim$x, starting.q=0.5, s=500)))
```

    + All the variation here is from simulation to simulation



