# Simulation and Stochastic Processes

Supposing we want to take a sample from the distribution $p(x)=cp^*(x)$ where $c$ is the normalization constant and $p^*(x)$ is given by

$$p^*(x)=0.5\exp(-(x-\mu_1)^2)+0.5\exp(-(x-\mu_2)^2).$$

Construct a rejection algorithm, that uses a function from a normal distribution with mean $μ$ and variance $σ^2$. Identify the μ and σ, for given values of $μ_1$ and $μ_2$. Check your results comparing them with a more simple algorithm that takes a sample from $p(x)$ using the method of synthesis.

What I have done so far

-Supposing $g=N(\mu,\sigma^2)$ . -I have found the $\frac{p^*(x)}{g(x)}$.

-Then I found the derivative and calculate $\frac d{dx}\frac{p^*(x)}{g(x)}=0$ in order to found the $x$ that verifies my equation. But here is my basic problem: if my calculations are right I find two solutions $x_1$ and $x_2$. I don't know how I should proceed then, in order to find the $M$ (the upper bound of $\frac{p^*(x)}{g(x)}$, because i have two $x$... should i use both of them?

-I have found also, $\mu=\sigma^2(\mu_1+\mu_2)$.

I would appreciate any help/tip.

-
Do you want to construct a rejection algorithm on purpose or do you just want to sample from your distribution? Since it is a mixture of Gaussians, there is an easier way to do that. –  fabee Nov 18 '11 at 15:25
@fabee unfortunately i have to construct a rejection algorithm because i also have to find the μ,σ.And after this i have to compare this results with the ones that'll take from the algorithm that i'll construct with the method of synthesis –  Panagiotis Nov 18 '11 at 18:15
Of course you should find two solutions: the problem is symmetric around $x = (\mu_1+\mu_2)/2$. To see this more clearly, change from the variable $x$ to a new variable $y - (\mu_1+\mu_2)/2$. This immediately implies $\mu=0$ (in the new coordinates), reducing the problem to finding $\sigma$. (To determine this, start by looking at what happens for very large $x$: that places a lower bound on $\sigma$.) –  whuber Nov 18 '11 at 23:29
So the two solutions that i found are probably correct.Then i equated these solutions and i found the μ.Usually when constructing a rejection algorithm the next step is to place the solution of x that you found in the original p*(x)/g(x),and from there you can find the upper bound(M).But in this case where i have two solutions i don't know how to proceed.I'not sure which x to choose.And also when i place one of the solutions the equation that comes of, is complicated and i can't solve it. –  Panagiotis Nov 19 '11 at 15:40
@fabee p(x) as you said is a mixture of Gaussians, how can i take a sample from them? –  Panagiotis Nov 19 '11 at 16:45