Conditional distribution of the weight of a mixture gaussian with data augmentation using gibbs sampling This question is relate to Differenciate between two distributions using gibbs sampling .
For $t=1,\,\dots,\,n$, let's $r_t\sim\mathcal{N}(0,\,\sigma_t^2)$ and $$\sigma_t^2=\left\{\begin{array}{lcl}
\sigma^2 & \text{with probability} & p\\
1 & \text{with probability} & 1-p
\end{array}\right.$$ The sample is $r_t, t=1,\,\dots, n$, we suppose we the value of $\sigma$ is known and our aim is to find sample from the posterior distribution of the weight $p$ using Gibbs sampling.
To do so, I follow the data augmentation algorithm as proposed by Diebolt and Robert (1994) at their section 3.1.1 (Example 2). The algorithm work fine if $\sigma$ is far from $1$ but it does not work anymore for $\sigma$ close to $1$. The following graph show what I found for $\sigma=3$ and $\sigma=0.9$ :

Problem :
This algorithm work fine for $\sigma^2$ far from $1$. But, when $\sigma^2$ is close to $1$, the algorithm doesn't not work.
My questions :

*

*Is there another algorithm to better the sampling of the posterior of $p$?

*Theoretically, is it link to the sample size of the number of iteration of my algorithm? In other words, do you think that If I raise the sample size $n$ or the number of iteration $n_{iter}$ I will find the better the sampling of the posterior of $p$?

*What will be the difference between $1$ and $\sigma$ such that the posterior of $p$ can be found properly? Using my empirical trials, $\sigma^2=0.1$, it works, but when $\sigma^2=0.9$ or $\sigma^2=0.8$ it doesn't work.

All kind of contribution will be very helpful. Thanks
Here is the code I've made using R software:
f.sample<-function(r_t,niter,sigma){
  #prior
  alpha_0 <- c(1,1)
  p <- rdirichlet(1,alpha_0)[1]
  
  n<-length(r_t)
  z_sim <- numeric(n)
  p_sim <- numeric(niter)
  
  for(iter in 1:niter){
    for(t in 1:n){
      w_1 <- p*dnorm(r_t[t],0,sigma)/(p*dnorm(r_t[t],0,sigma) + 
              (1-p)*dnorm(r_t[t],0,1))
      z_sim[t] <- rbinom(1,1,w_1)
    }
    n_iter <- sum(z_sim==1)
    
    alpha_post   <- alpha_0+c(n_iter,n-n_iter)
    p            <- rdirichlet(1,alpha_post)[1]
    p_sim[iter]  <- p
  }
  
  return(mean(p_sim[(floor(niter/2)+1):niter]))
  #plot(p_sim,type="l")
}


#Monte Carlo
tru.p <- 0.10
n<- 1000
f.boot<-function(nboot,niter,nseed,sigma){
  set.seed(nseed)
  seeds <- sample.int(n=10^8,size=nboot)
  p_sim <- numeric(nboot)
  for(boot in 1:nboot){
    set.seed(seeds[boot])
    z_t   <- rbinom(n,1,tru.p)
    r_t   <- rnorm(n,mean=0,sd=c(1,sigma)[z_t+1])
    p_sim[boot] <- f.sample(r_t,niter,sigma)
    if(boot%%50==0)  print(p_sim[boot])
  }
  
  return(p_sim)
}

p_3<-f.boot(nboot = 500, niter=50, nseed=3216, sigma=3)
p_0.9<-f.boot(nboot = 500, niter=50, nseed=3216, sigma=0.9)

As you can see, I've made a bootstrap in other to insure the mean of the sampling is not only relate to the sample.

Here are a few runs with different values of $p$ and $\sigma$ (in the title of the plots) and $n=100$ realisations for each:

The histogram is based on $10^4$ Gibbs steps and the curve is the exact likelihood, renormalised to fit the height of the histogram. Whatever the values of $(p,\sigma)$ the fit is quite acceptable and indicate that the convergence properties of the Gibbs sampler are satisfactory. A random walk Metropolis-Hastings algorithm could be used instead but there is no argument to believe it would do better.
Concerning the other questions,


*The number of MCMC iterations is not directly connected with the sample size $n$. As $n$ grows, the likelihood gets peakier (around the true value of the parameter) which means it may be be harder to reach, but also that it is smoother (asymptotically Normal) and thus easier to explore.


*The posterior (equivalent to the likelihood in this case) is found properly, as shown by the pictures.
