2
$\begingroup$

I have this $f$ function below.

$$ f(x_1,x_2)\propto \left(\dfrac{x_1}{x_2}\right)\left(\dfrac{\alpha}{x_2}\right)^{x_1-1}exp\left\{-\left(\dfrac{\alpha}{x_2}\right)^{x_1} \right\}I_{R^+}(x) $$ where $\alpha > 0$.

I would like to estimate expected values and variances of $x_1$ and $x_2$ by using Gibbs Sampling or Metropolis-Hastings algorithms.

First of all, I tried to run Metropolis-Hastings algorithm, however, finding proper proposal distributions for both random variables is too hard. I tried some known distributions such as Gamma, Weibull or Exponentials, because they have the support of $x_1$ and $x_2$, I did not get any good results, considering the plots of Markov Chains I got were too bad.

Another solution that I tried was Gibbs sampling. I tried to find full conditionals, i.e., $f(x_1|x_2)$ and $f(x_2|x_1)$ which are $x_1$ given $x_2$ and $x_2$ given $x_1$, respectively.

I found $f(x_1|x_2)$ as follows:

$$ f(x_1|x_2)\propto x_1 \beta^{x_1}\exp\left\{-\left(\beta \right)^{x_1} \right\} $$ where $\beta=\dfrac{\alpha}{x_2}$ and $0<x_1<\infty$.

I tried to derive a standard distribution by using transformation. Thus, I made the change of variable with $\left(\beta \right)^{x_1}=y$. In the end I get the function below.

$$ g(y|x_2)\propto ye^{-y}e^{-e^{-y}} $$

The last two terms, here, resemble to a Gumbel Distribution with shape parameter $\mu=0$ and scale parameter $B=1$. Thus, Gumbell distribution has a support of $(-\infty, +\infty)$, that was not the case because, here $y$ has a support of $(-\infty, 0)$. So, I thought about truncating the Gumbel distribution, and then using acceptance rejection method to get random numbers of $y$'s. However, that did not work out, either because according to this method, the unknown $g(x)$ function, which is $y$ here, must be within $(0,1)$ to compare it with a uniform random number. So, I left it, here and thought maybe I can run MH for $x_1$ within Gibbs sampling. For that, once again, I truncated Gumbel distribution but this time for the support $(0, \infty)$. Because I assumed, since I can derive somewhat truncated Gumbel distribution for $f(x_1|x_2)$, this truncated Gumbel distribution might be a good proposal distribution for $x_1$, overall.

I got the probability distribution function (pdf) of truncated Gumbel distribution as given in below:

$$ f(x)=\frac{\dfrac{1}{B}e^{-{\dfrac{(x-\mu)}{B}+e^{-\dfrac{(x-\mu)}{B}}}}}{1-e^{-e^{\frac{\mu}{B}}}} $$

and its cumulative distribution function (cdf) given as follows:

$$ F(x)=\dfrac{e^{-{e^{-\dfrac{(x-\mu)}{B}}}}-e^{-e^{\frac{\mu}{B}}}}{1-e^{-e^{\frac{\mu}{B}}}} $$

Since I need to generate new random variates for MH algorithm, I need the inverse of this cdf. I will generate a new $z$ from $Uniform(0,1)$ and put it in the inverse of CDF of this distribution to get a random variates from truncated Gumbel distribution. The inverse cdf is given below:

$$ x=-B~ln\left\{-ln\left[z\left( 1-e^{-e^{\frac{\mu}{B}}} \right) + e^{-e^{\frac{\mu}{B}}} \right] \right\}+\mu $$

Anyway, now it's turn for $f(x_2|x_1)$ which is $x_2$ given $x_1$.

$$ f(x_2|x_1)=\dfrac{1}{{x_2}^{x_1}}\exp\left\{-\dfrac{\gamma}{{x_2}^{x_1}} \right\} $$ where $\gamma=\alpha^{x_1}$ and $0<x_1<\infty$. Thus, I made the change of variable with $\dfrac{1}{{x_2}^{x_1}}=t$. In the end I get the function below.

$$ g(t|x_1)\propto \exp\{-\gamma ~ t \}~t^{-\dfrac{1}{x_1}} $$

At first, I thought this was a Gamma distribution with the shape parameter $a=-\dfrac{1}{x_1}+1$ and scale parameter $b=\gamma=\alpha^{x_1}$. It seems, no doubt, a gamma distribution just because here $x_1$ is a constant because it's given/known. But the problem comes when it is applied. Because when I do Metropolis-Hastings within Gibbs, things get a bit complicated. Since I am going to update $x_1$ and put it in conditional distribution $g(t|x_1)$ in order to get $t$'s random variates,so then, I am going to inverse the transformation $\dfrac{1}{{x_2}^{x_1}}=t$ in order to get $x_2$'s. However, the parameters of gamma distribution creates a big problem.

As you know Gamma distribution has two parameters which is greater than zero. Here, $b=\gamma=\alpha^{x_1}$ is not a big deal, because $0<x_1<\infty$ so $b>0$. HOWEVER, in order to get $a>0$, the shape parameter, I need to generate $x_1>1$ which is impossible because $0<x_1<\infty$.

I realized this problem long after I wrote and run the code, which giving me "NaN" for $x_2$ because the problem I talked above. I did what I said. I gave truncated Gumbel distribution as a proposal distribution for $x_1$ and run Gamma distribution for full conditional of $x_2$.

My initial values are 1 for each of $x_1$ and $x_2$. I am calculating acceptance ratio for $x_1$ since I run Metropolis-Hastings on it. I update each $x_1$ for each iteration in the parameters of Gamma distribution for $x_2$.

M=1000;
x=zeros(M+1,2);
t=zeros(M,1); %I create a Mx1 vector for inversing the transformation of x_2
alpha=2; % arbitrarily picked alpha value for my function
mu=0; %the shape parameter of truncated Gumbel distribution
beta=1; %the scale parameter of truncated Gumbel distribution

fx=@(x1,x2)((x1/x2)*((alpha/x2)^(x1-1))*\exp(-1*((alpha/x2)^x1))); %the function I have
gx=@(x1)(((1/beta)*\exp(-(((x1-mu)/beta)+\exp(-((x1-mu)/beta)))))/(1-\exp(-\exp(mu/beta)))); %the pdf of truncated Gumbell distribution

alpha21=@(x1)((-1/x1)+1); %the shape parameter of Gamma distribution for x2
beta21=@(x1)(alpha^(x1)); %the scale parameter of Gamma distribution for x2

acceptance_ratio=0;

t(1,1)=1;
x(1,1)=1;
x(1,2)=1;

for i=2:(M+1)
        old=fx(x(i-1,1),x(i-1,2));
        oldprob=gx(x(i-1,1)); %finding probabilities of x1's from truncated Gumbel distribution
        ex=x(i-1,1);

        z=unifrnd(0,1);
        t=(-beta)*log(-log(z*((1-\exp(-\exp(mu/beta))))+\exp(-\exp(mu/beta))))+mu; %inverse of cdf of truncated Gumbel distribution

        x(i,1)=t;
        new=fx(x(i,1),x(i-1,2));
        newprob=gx(x(i,1));
        acceptance_prob=min(1,new*newprob/(old*oldprob));
        if unifrnd(0,1)>acceptance_prob
            x(i,1)=ex;

        else
            acceptance_ratio=acceptance_ratio+1;
        end

        t(i,1)=gamrnd(alpha21(x(i,1)),beta21(x(i,1))); %updating the parameters of Gamma distribution in each iteration
        x(i,2)= t(i,1)^(-1/x(i,1));
end
mean(x)
var(x)
AR=acceptance_ratio/M

I think I am right at what I am doing. I know my problem. But I am stuck at it. I cannot find any solution. Would you care to help me?

$\endgroup$
7
$\begingroup$

Are you sure the joint density$$f(x_1,x_2)=\left(\dfrac{x_1}{x_2}\right)\left(\dfrac{\alpha}{x_2}\right)^{x_1-1}\exp\left\{-\left(\dfrac{\alpha}{x_2}\right)^{x_1} > \right\}\mathbb{I}_{\mathbb{R}^*_+}(x_1,x_2)$$ is integrable?


When I consider the conditional$$f(x_2|x_1)=\dfrac{1}{{x_2}^{x_1}}\exp\left\{-\dfrac{\gamma}{{x_2}^{x_1}} \right\}$$ it should be a proper probability density if the joint above is a proper density. However, when considering the case $x_1=1$, it simplifies to $$f(x_2|x_1=1)=\dfrac{1}{{x_2}}\exp\left\{-\dfrac{\gamma}{{x_2}} \right\}$$ which does not integrate over $(0,\infty)$ since the change of variable $\eta=1/x_2$ leads to the density $$f(\eta|x_1=1)=\eta\exp\left\{-\gamma\eta\right\}\eta^{-2}.$$ Which is equivalent to $\eta^{-1}$ when $\eta\approx 0$. Your change of variable $t=x_2^{-x_1}$ shows the same issue occurs when $x_1\le 1$.

Note: When you operate the change of variable $$x_1\longrightarrow\left(\beta \right)^{x_1}=y$$ the new variate $y$ is either supported by $(0,1)$ or $(1,\infty)$ depending on whether or not $\beta<1$. But this is irrelevant relative to the main issue.

$\endgroup$
  • $\begingroup$ Professor, I have been told that this is a real density but without a normalizer constant. My mistake is that I forgot to put proportionality symbol. $\endgroup$ – ARAT Jan 17 '15 at 20:19
  • 2
    $\begingroup$ If all there is missing is only a normalising constant, this does not modify integrability or lack thereof. So I am afraid the density as given is not integrable. $\endgroup$ – Xi'an Jan 17 '15 at 20:28
  • 1
    $\begingroup$ this function is not integrable so this mean it is not bounded, right? $\endgroup$ – ARAT Jan 18 '15 at 8:54
  • $\begingroup$ warning: non-integrable functions may be bounded, take for instance $f(x)=1/(1+|x|)$. $\endgroup$ – Xi'an Jan 18 '15 at 10:30
  • $\begingroup$ Professor, I found out that this function is the density of Weibull distribution, whose parameters are random variables. Thus I think I found my solution. I can use two full conditionals for $x_1$ and $x_2$ which they are Gamma distribution and Logarithmically concave distribution. $\endgroup$ – ARAT Jan 20 '15 at 8:30
2
$\begingroup$

$f(x_1,x_2)$ is not a density function since it is not integrable.

Note that

$$\int_0^{\infty}\int_0^{\infty} f(x,y)dydx \geq \int_1^{\infty}\int_0^{\infty} f(x,y)dydx = \int_1^{\infty}\Gamma\left(\dfrac{x-1}{x}\right)dx = \infty,$$

where $\Gamma$ is the Gamma function. The first inequality follows by using that you are integrating a positive function over a smaller domain. The following equality is obtained after some algebra (you can also use Mathematica if you don't believe in magic). The integrand is lower bounded by $1$, which implies the infite integral.

Perhaps there is a typo in your expression since you mentioned that it is trivial to show it is a density? It would be interesting to see the original source.

$\endgroup$
  • $\begingroup$ Well you are right. It seems to be improper. However I think by using a proper prior you can get a proper posterior, Am I wrong? Here's a solution to bayesian analysis of weibull distribution goo.gl/YNfwVy $\endgroup$ – ARAT Feb 5 '15 at 9:13

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.