2
$\begingroup$

Assume X and Y have a bivariate lognormal distribution (x,y>0) that is:

$f_{X,Y}(x,y)$=$$\frac{1}{2p\sqrt{1-r^2}xy\sigma_1\sigma_2}exp\{\frac{-1}{2(1-r^2)}[(\frac{ln(x)-\mu_1}{\sigma_1})^2-2r(\frac{ln(x)-\mu_1}{\sigma_1})(\frac{ln(y)-\mu_2}{\sigma_2})+(\frac{ln(y)-\mu_2}{\sigma_2})^2]\}$$
I want to know the median of |x-y|.
first I calculated the density of |x-y| based on changing the variables, as follows:

I took u=|x-y| and v=y then

\begin{cases} x=u+v, & \text{if x>y which result in: u>0 , v>0} \\ x=v−u, & \text{if x<y which result in v>u>0} \end{cases} The Jacobian is one, so the density function of u=|x-y| is:
$$f_{U}(u)=\int_{0}^\infty f_{X,Y}(u+v,v)dv+\int_{u}^\infty f_{X,Y}(v-u,v)dv$$

my question is: how can I calculate the median of the above density using numerical integration(assume for the simplest case where $\mu=\begin{pmatrix} 1 \\ 1 \\ \end{pmatrix} and \sigma= \begin{pmatrix} 1 & 0.95 \\ 0.95 & 1 \\ \end{pmatrix}$)?
I am quite new in R..I really appreciate any help

$\endgroup$
  • $\begingroup$ If you ask a numerical question (using numerical integration), you need to provide numerical parameter values and the functional form. Also see related stats.stackexchange.com/questions/152850/… and your structure is more complicated via dependence. $\endgroup$ – wolfies Feb 2 '17 at 3:32
  • $\begingroup$ I made some changes regarding your comment. is it clear now? $\endgroup$ – Ella Feb 2 '17 at 4:40
1
$\begingroup$

A simple Monte Carlo approach could work. This method approximates the distribution using a finite set of samples, then computes the median based on the samples. Draw $n$ independent samples of $X$ and $Y$ from the joint distribution: $\{(x_1, y_1), ... (x_n, y_n)\}$. I'll explain how do do this below. This gives samples of $U$: $\{u_1, ..., u_n\}$, where $u_i = |x_i - y_i|$. Take the median of these samples. Accuracy will increase with $n$, which you could easily set in the thousands or millions.

Your problem has special structure that makes it possible to sample $X$ and $Y$ directly from the joint distribution. $\log X$ and $\log Y$ have a bivariate normal distribution from which you can easily sample using a multivariate Gaussian random number generator (using the appropriate mean and covariance matrix). This blog post describes how to do this in R. Once you have samples of $\log X$ and $\log Y$, exponentiate these values to obtain samples of $X$ and $Y$. It's not always possible to sample directly from a given distribution, in which case methods like MCMC algorithms can be used. Luckily, this isn't necessary in your case.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you, I am a little confused. how can I get the $u_i$ samples from MCMC algorithm? Would you mind providing some codes? $\endgroup$ – Ella Feb 4 '17 at 18:16
  • $\begingroup$ There's actually an easy way to do sample directly. See my edited answer. $\endgroup$ – user20160 Feb 6 '17 at 5:16

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.