# Dividing a uniform by a normal random variable: What's the distribution?

I have two random variables $s\sim \mathcal{N}(\nu,\sigma)$ and $a\sim \mathcal{U}(0,A)$, $0<A<1$, calculate a third r.v. $t=(1-a)/s$ and want to find its distribution $p(t|\nu,\sigma,A)$.

My reasoning is, that for any value of $a$, there is exactly one $s=(1-a)/t$ such that the $(a,s)$ pair will produce $t$ and I will therefore have to integrate over the probabilities for these values:

$$p(t|\nu,\sigma,A) = \int_0^A p_a(x)p_s((1-x)/t)\,\mathrm dx\\$$

which can be expressed as the sum of two error functions.

But already when checking this step numerically, I get a discrepancy between estimated (by sampling) and calculated distribution:

Why do the histogram and the calculated distribution do not match?

This is the code I used:

n=10000000

A=0.7
nu=.7
sigma=.1

# sampling from the target distribution
s=rnorm( n, mean=nu, sd=sigma )
a=runif( n, min=0, max=A )
t=(1-a)/s
hist(t,200, freq=F, xlim=c(0,5), ylim=c(0,2.5))

# plot the analytical result
analytic <- function(t, A, nu, sigma ){
tmp <- function(x, A, t, nu, sigma){
return (1/A*dnorm( (1-x)/t, mean=nu, sd=sigma) )
}
return( integrate( tmp, 0, A, A=A, t=t, nu=nu, sigma=sigma)$value) } x=seq(0,5,by=.01) y=rep(0,length(x)) for( i in seq(1,length(x)) ){ y[i]=analytic(x[i], A, nu, sigma) } lines( x, y, col="red", type="l")  - could someone include the image into the question? Thanks – thias Jul 19 '12 at 13:06 You have sufficient reputation that you should be able to include the image yourself. Are you encountering difficulties with that? – cardinal Jul 19 '12 at 13:09 ah, sorry - I read somewhere that you need 2k reputation but that may be only in SO :-) – thias Jul 19 '12 at 13:11 The ratio distribution can be calculated using this. – user10525 Jul 19 '12 at 13:15 By the way, in the problem title, you mention dividing a uniform by a normal, but in the problem statement it looks like you're dividing a normal by a uniform. – Macro Jul 19 '12 at 13:37 show 6 more comments ## 1 Answer Let$X$and$Y$be two independent random variables and define$Z=\frac{X}{Y}$. Consider the change of variable$(X,Y)\leftrightarrow (Z,Y)$, then the corresponding inverse transformation is$(ZY,Y)$. Therefore, the Jacobian matrix is given by (see this and this links for a reference of this procedure) \begin{eqnarray} J=\left( \begin{array}{cc} \dfrac{\partial X}{\partial Z}& \dfrac{\partial Y}{\partial Z}\\ \dfrac{\partial X}{\partial Y}& \dfrac{\partial Y}{\partial Y} \\ \end{array} \right)=\left( \begin{array}{cc} Y& 0\\ Z& 1 \\ \end{array} \right). \end{eqnarray} The absolute value of the determinant of$J$is$\vert \mbox{det}(J)\vert=\vert Y\vert$. Then the density of$Z$is given by $$\newcommand{\rd}{\,\mathrm d}f_Z(z)=\int_{-\infty}^{\infty}f_{X,Y}(zy,y)\vert y \vert \rd y. \tag{\star}$$ Using independence you have that$f_{X,Y}=f_X f_Y$. Now, note that in your case$1-a\sim U(1-A,1)$. Then, by replacing the corresponding densities in$(\star)$you obtain $$f_Z(z)=\int_{-\infty}^{\infty}\vert y \vert\dfrac{I_{(1-A,1)}(zy)}{A} \varphi(y;\nu,\sigma) \rd y.$$ Here,$I_{(1-A,1)}(zy)$gives you the integration domain of$y$which depends on the sign of$z$. These are$\left(\dfrac{1-A}{z},\dfrac{1}{z}\right)$for$z>0$and$\left(\dfrac{1}{z},\dfrac{1-A}{z}\right)$for$z<0$. The density is$0$for$z=0$. Therefore, if you replace analytic = function(z){ if(z>0) return( integrate(function(x) abs(x)*dnorm(x,nu,sigma), (1-A)/z, 1/z)$value/A)
else return( integrate(function(x) abs(x)*dnorm(x,nu,sigma), 1/z, (1-A)/z )\$value/A)
}


and y[i]=analytic(x[i]) in the loop, you will observe the desired fit.

NB Your argument did not work because it is not a formal change of variable, in particular it does not include the Jacobian.

I hope this helps.

-
thank you for the great answer! Given your previous hints, I managed to work through it myself as well... –  thias Jul 19 '12 at 15:42