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

Question

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")

Answer 1

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.