How can I generate random samples from following probability density function?

76 Views Asked by At

Let $\mathbf{\alpha}=(\alpha_1, \ldots, \alpha_m)$. The posterior density function of $\mathbf{\alpha}$ is given by $$h_0(\mathbf{\alpha}|\mathbf{x})=‎\frac{\prod_{i=1}^{m}\alpha_i^{a_i}}{\left(1+\sum_{i=1}^{m}\alpha_i\right)^{\sum_{i=1}^{m+1}a_i}}‎\prod_{j=1}^{2}\Big[\sum_{i=1}^{m}\alpha_i n_i(x_i-x_{i-1})+d_j\Big]^{-(m_j+c_j)}$$ where $a_i, x_i, d_j, m_j, c_j >0$and are known. I used the Metropolis Hastings algorithm to generate samples from $h_0$. I have taken multivariate normal distribution as the proposal density. That is I have taken a multivariate normal distribution $N_m({\mathbf{\alpha}},H^{-1}(\hat{\mathbf{\alpha}}) )$, where $\hat{\mathbf{\alpha}}$ is the posterior mode and $$H(\hat{\mathbf{\alpha}})=\Big(-\frac{\partial^2 \log h_0(\mathbf{\alpha|x})}{\partial \alpha_i \partial \alpha_j}|_{\mathbf{\alpha}=\hat{\mathbf{\alpha}}} \Big)$$ is the precision matrix . It may be noted that if I generate observation from multivariate normal distribution, I may get negative observations which are not acceptable for $\mathbf{\alpha}.$ Keeping this in mind, the M-H algorithm steps are given below.

Step 1. Set initial values $\mathbf{\alpha}^{(0)}=\hat{\mathbf{\alpha}}.$

Step 2. For $t = 1, . . . , T$ repeat the following steps.

a) Set $\mathbf{\alpha}=\mathbf{\alpha}^{(t-1)}.$

b) Generate new candidate parameter values $\mathbf{\theta}$ from $N_m(\mathbf{log(\mathbf{\alpha})}, H^{-1}(\hat{\mathbf{\alpha}})).$

c) Set $\mathbf{\alpha}^\prime=\exp(\mathbf{\theta})$.

d) Calculate $p=\min\left\{1,\frac{h_0(\mathbf{\alpha}^\prime|\mathbf{x})\prod_{i=1}^{m}\alpha_i^\prime}{h_0(\mathbf{\alpha}|\mathbf{x})\prod_{i=1}^{m}\alpha_i} \right\}.$

e) Update $\mathbf\alpha^{(t)}=\mathbf{\alpha}^\prime$ with probability p, otherwise set $\mathbf\alpha^{(t)}=\mathbf{\alpha}.$

But the above algorithm is not efficient.