Difference of two squared normal dependent variables I need to find the distribution of the random variable $Z$
$Z = \frac{(X - \mu_0)^2}{\sigma_0^2} - \frac{(X - \mu_1)^2}{\sigma_1^2}$, where $X  \sim \mathcal{N}(\mu_0, \sigma_0)$.
We can find the distribution of each variable. The first variable $\frac{(X - \mu_0)^2}{\sigma_0^2} \sim \chi_1$, and the second variable $\frac{(X - \mu_1)^2}{\sigma_1^2} \sim \frac{\sigma_0^2}{\sigma_1^2}\chi_1(\lambda = \frac{(\mu_0 - \mu_1)^2}{\sigma_0^2})$. Thus, the distribution of $Z$ is weighted sum of non-central $\textbf{dependent}$ chi square distribution. I find a package in R (sadists) for $\textbf{independent}$ weighted sum to do inference. How can I find the distribution of $Z$ or at least a package to do some inference?
 A: \begin{align}
Z&= \frac{(X-\mu_0)^2}{\sigma_0^2}-\frac{(X-\mu_1)^2}{\sigma_1^2}\\
&= \left(\frac{1}{\sigma_0^2} - \frac{1}{\sigma_1^2}\right)X^2 
-2\left(\frac{\mu_0}{\sigma_0^2} - \frac{\mu_1}{\sigma_1^2}\right)X
+\left(\frac{\mu_0^2}{\sigma_0^2} - \frac{\mu_1^2}{\sigma_1^2}\right)\\
&= \left(\frac{1}{\sigma_0^2} - \frac{1}{\sigma_1^2}\right)
\left(X-A\right)^2 + B
\end{align}
where the last step is obtained via "completing the square" and the 
values of $A$ and $B$ are left for you to determine 
(it is a tedious exercise
at the level of middle-school algebra). So $Z$ is a displaced and 
scaled noncentral
$\chi_1^2$ random variable. Note that the scale factor is a negative
number when $\sigma_0^2 > \sigma_1^2$, and that if 
$\sigma_0^2 =\sigma_1^2$, then $Z$ is not a $\chi_1^2$ random variable
at all but a normal random variable! 
You can work with this instead of calling
on the sadists (?) who run R for help.
A: The distr package for R lets you create distributions based on functions of known distributions.  So you could use it to create a difference of 2 Chi-squares and do various computations based on the new distribution.
edit
Here is a basic example with values for the means and standard deviations chosen:
> library(distr)
> 
> x <- Norm(1, 2)
> z <- ((x-1)/2)^2 - ((x-3)/4)^2
> 
> plot(z)
> distr::q(z)(c(0.025,0.5,0.975))
[1] -1.738  0.121  4.612
> p(z)(0)
[1] 0.408
> mean(r(z)(10000))
[1] 0.503

