# How to model a biased coin with time varying bias?

Models of biased coins typically have one parameter $\theta = P(\text{Head} | \theta)$. One way to estimate $\theta$ from a series of draws is to use a beta prior and compute posterior distribution with binomial likelihood.

In my settings, because of some weird physical process, my coin properties are slowly changing and $\theta$ becomes a function of time $t$. My data is a set of ordered draws i.e. $\{H,T,H,H,H,T,...\}$. I can consider that I have only one draw for each $t$ on a discrete and regular time grid.

How would you model this? I'm thinking of something like a Kalman filter adapted to the fact that hidden variable is $\theta$ and keeping the binomial likelihood. What could I use to model $P(\theta(t+1)|\theta(t))$ to keep inference tractable?

Edit following answers (thanks!): I would like to model $\theta(t)$ as a Markov Chain of order 1 like it is done in HMM or Kalman filters. The only assumption I can make is that $\theta(t)$ is smooth. I could write $P(\theta(t+1)|\theta(t)) = \theta(t) + \epsilon$ with $\epsilon$ a small Gaussian noise (Kalman filter idea), but this would break the requirement that $\theta$ must remain in $[0,1]$. Following idea from @J Dav, I could use a probit function to map the real line to $[0,1]$, but I have the intuition that this would give a non-analytical solution. A beta distribution with mean $\theta(t)$ and a wider variance could do the trick.

I'm asking this question since I have the feeling that this problem is so simple that it must have been studied before.

• You can get an estimate if you have a model for how the success proportion chanhes with time. Many different models would work and the estimates could vary a lot based on the assumed model. I do not think tractability is a practical criterion for choosing a model. I would want to understand the process and look for a model thaat demonstrates characteristic that agree with the behavior you expect. Jul 25 '12 at 11:42
• @MichaelChernick : Thanks. The only assumption I can make is that $\theta$ is moving smoothly and slowly. Moreover tractability is an important criteria since I actually want to extend the solution to multivariate case with non-trivial inter-dependencies. An ideal solution would be analytical and give 'online' update of parameter estimates when a new data arrives. Jul 25 '12 at 13:47
• Can you quantify what you mean by "$\theta$ is moving smoothly and slowly?" The integers are discrete, and there are smooth functions which take on arbitrary values on the integers, which means that smoothness gives no constraints. Some notions of "slowly" still don't give any constraints, while some do. Jul 25 '12 at 16:35
• How fast is "slowly", like a change in probability of 0.1 / unit time or 0.001 or... And how long a sequence do you expect to have? Is the range relatively narrow (e.g, 0.2 - 0.4) or does it come close to (0,1)? Jul 25 '12 at 19:58
• @DouglasZare By 'smooth', I wanted to state that E[θ_t+1|θ_t]=θ_t (or very close) and VAR(θ_t+1|θ_t) is small. θ is not jumping around (otherwise nothing could be done really). Jul 26 '12 at 12:19

I doubt you can come up with a model with analytic solution, but the inference can still be made tractable using right tools as the dependency structure of your model is simple. As a machine learning researcher, I would prefer using the following model as the inference can be made pretty efficient using the technique of Expectation Propagation:

Let $X(t)$ be the outcome of $t$-th trial. Let us define the time-varying parameter

$\eta(t+1) \sim \mathcal{N}(\eta(t), \tau^2)$ for $t \geq 0$.

To link $\eta(t)$ with $X(t)$, introduce latent variables

$Y(t) \sim \mathcal{N}(\eta(t), \beta^2)$,

and model $X(t)$ to be

$X(t) = 1$ if $Y(t) \geq 0$, and $X(t) = 0$ otherwise. You can actually ignore $Y(t)$'s and marginalize them out to just say $\mathbb{P}[X(t)=1] = \Phi(\eta(t)/\beta)$, (with $\Phi$ cdf of standard normal) but the introduction of latent variables makes inference easy. Also, note that in your original parametrization $\theta(t) = \eta(t)/\beta$.

If you are interested in implementing the inference algorithm, take a look at this paper. They use a very similar model so you can easily adapt the algorithm. To understand EP the following page may found useful. If you are interested in pursuing this approach let me know; I can provide more detailed advice on how to implement the inference algorithm.

To elaborate on my comment a model such as p(t)=p$_0$ exp(-t) is a model that is simple and allows estimation of p(t) by estimating p$_0$ using maximum likelihood estimation. But does the probability really decay exponentially. This model would be clearly wrong if you observe time periods with high frequency of success than you observed at earlier and later times. Oscillatory behavior could be modelled as p(t)=p$_0$ |sint|. Both models are very tractable and can be solved by maximum likelihood but they give very different solutions.

• It appears that the OP is looking to model the success probability at time $t$, $\theta(t)$, as a markovian process, not to specify some functional form for $\theta(t)$. Jul 25 '12 at 12:15
• @macro is right, I'm not able to provide a parametric form for $theta(t)$, and this is not desirable since this function could be anything smooth. I want a Markov model of order 1 similar to a Hidden Markov Model or a Kalman filter, but with a hidden variable that takes real values between 0 and 1, and with a Bernouilli likelihood. Jul 25 '12 at 13:25
• @pierre Okay prior to the edit it appeared that you were looking to estimate the time varying p and were just suggesting the HMM as a possible approach. I was not recommending a functional form for the way it changes with t. I was making a point that without further information many models of various types could be constructed and my two examples were to show that without further information model choices could give very different answers. Why would you insist on an HMM? If one worked and fit your data why reject it because it is "non-analytical. Jul 25 '12 at 14:03
• I am suggesting that finding convenient solutions is not the way to solve practical statistical problems! Jul 25 '12 at 14:03
• @MichaelChernick Lastly: I would like to find an analytical solution since I hope this is a well-known problem and people proposed flexible enough analytical solution. But I agree with our suggestion that modelling the 'real dynamics' is more important than the computational cost in general. Sadly this is for big data and a slow algo will be useless :-( Jul 25 '12 at 14:34

Your probability changes with $t$ but as Michael said, you don't know how. linearly or not ? It looks like a model selection problem where your probablity $p$ :

$p=\Phi(g(t,\theta))$ may depend on a highly non linear $g(t,\theta)$ function. $\Phi$ is just a bounding function that guarantees between 0 and 1 probabilities.

A simple exploratory approach would be to try several probits for $\Phi$ with different non linear $g()$ and to perform a $g()$ model selection based on standard Information Criterias.

As you said using probit would imply numerical solutions only but you may use a logistic function instead :

Logistic function: $P[\theta(t+1)] = \frac{1}{1+\exp{(\theta(t)+\epsilon)}}$

Linearized by : $\log{\frac{P}{1-P}} = \theta(t)+\epsilon$

I'm not sure how this can work under Kalman filter approach, but still believe that a non linear specification like $\theta(t+1)=a t^3 +bt^2+ct + d$ or many others without a random term will do the job. As you can see this function is "smoth" in the sense that it's continous and differentiable. Unfortunately adding $\epsilon$ would generate jumps of the resulting probability which is something you don't want so my advice would be to take out $\epsilon$.

Logit probablity: $P[Coin_{t+1}=H | t] = \frac{1}{1+\exp{(\theta(t))}}$

You already have randomnes in the bernoulli event (Markov Chain) and you are adding an additional source of it due to $\epsilon$. Thus, your problem could be solved as a Probit or Logit estimated by Maximum likelihood with $t$ as explanatory variable. I suppose you agree that that parsimony is very important. Unless your main objective is to apply a given method (HMM and Kalman Filter) and not to give the simplest valid solution to your problem.

• If you use a probit, a multivariate extension is straightforward as a multivariate probit can be estimated. Dependencies would be implicit by the covariance matrix of the implied multivariate normal distribution.
– JDav
Jul 27 '12 at 19:48