# Mixture of normals, dependent on past

I have the following probability model:

$(X_k|\text{PastHistory}_{k-1}, \theta_0,\theta_1,\theta_2) \sim (\pi\cdot N(\theta_1+\theta_0\cdot X_{k-1},1)+(1-\pi)\cdot N(\theta_2+\theta_0\cdot X_{k-1},1))$

and $X_k=0.5X_{k-1}+\epsilon_k$ with $\epsilon_k \sim N(0,1)$

I want to build a function that gives a sample of size n, from the mixture.

I've tried something like this in Mathematica:

  Sample[n_] := Module[{lst},
lst = {RandomVariate[NormalDistribution[0, 1/(1 - 0.5^2)]]};
element = 1;
While[element <= n,
AppendTo[lst,
RandomVariate[
TransformedDistribution[
p*X1 + (1 - p) X2, {
X1 \[Distributed]
NormalDistribution[ θ1 + θ0*lst[[element]], 1],
X2 \[Distributed]
NormalDistribution[ θ2 + θ0*lst[[element]], 1]
}]]];
element = element + 1;
];
lst
];


Is this correct?

Many thanks for any help. ;)

• I think I'm missing something, because it looks like you are giving two different ways of constructing the series $X_k$ (and they lead to different results): the method before the "and" and the method after it. You cannot do both at the same time! Which one do you want to implement? – whuber Feb 14 '14 at 15:41
• You're right. I want to implement the first one. – An old man in the sea. Feb 14 '14 at 15:51

We can write

$$X_k = \theta_0 X_{k-1} + \theta_1 U_{k-1} + \theta_2 (1 - U_{k-1}) + Z_{k-1}$$

where $U_i$ are iid Bernoulli$(\pi)$ variables and $Z_i$ are iid Normal$(0,1)$ variables. This form of expression allows us to capitalize on vectorized methods of generating large lists of random variates at once, thereby obtaining $(U_i)$ and $(Z_i)$ quickly and from them the vector $Y = \theta_1 U + \theta_2 (1-U) + Z.$ All that remains is to compute the $X_k$ iteratively, which in Mathematica-speak is a "fold" of the bivariate function

$$f(X, Y) = \theta_0 X + Y.$$

That leads directly to code like the following, whose arguments are the length of the simulated process $n$, the vector of parameters $(\theta_0, \theta_1, \theta_2),$ the Bernoulli parameter $\pi$ (here written $p$), and the starting value $x_0$:

simulate[n_, {θ0_, θ1_, θ2_}, p_: 0.5, x0_: 0] :=
With[
{u = RandomVariate[BernoulliDistribution[p], n],
z = RandomVariate[NormalDistribution[0, 1], n]},
FoldList[θ0 #1 + #2 &, x0, θ1 u + θ2 (1 - u) + z]
]


To test it, consider the process with no autocorrelation and widely separated means in the mixture:

n = 10^5;
AbsoluteTiming[sim = simulate[n, {0, -3, 3}, 1/3];]


{0.100006, Null}

That's a million values per second--not too bad. Let's look at some of these values; I plot every hundredth one in order to keep the plot visually recognizable:

ListPlot[sim[[1 ;; n ;; 100]], DataRange -> {1, n}]


Histogram[sim]


As intended, one-third ($p$) of the values are near $\theta_1=-3$ and the remaining two-thirds are near $\theta_2=3.$