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,
           p*X1 + (1 - p) X2, {
           X1 \[Distributed] 
           NormalDistribution[ θ1 + θ0*lst[[element]], 1], 
           X2 \[Distributed] 
           NormalDistribution[ θ2 + θ0*lst[[element]], 1]
     element = element + 1;

Is this correct?

Many thanks for any help. ;)

  • $\begingroup$ 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? $\endgroup$ – whuber Feb 14 '14 at 15:41
  • $\begingroup$ You're right. I want to implement the first one. $\endgroup$ – 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] := 
  {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}]



Figure 2

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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.