# Weird pdf of a quadratic function of a N(0,1) variable: miscoding or big rounding error?

I would like to calculate the pdf of a random variable y defined by :

y=c+b*x+a*x^2

The pdf is a non-central chi-squared distribution. For a>0, it should be equal to zero if y is smaller than d, where d=c-(b^2)/4a.

Strangely enough, when computing it with R, the pdf shoots up at y>d+e, where e is quite large.

Is there an error in my codes (below) or is it a rounding error? In the latter case, how to address it?

set.seed(101)
x <- seq(-3.5,3.5,length.out=1000)
c<-80
b<-30
a<-6
y<-c+b*x+a*(x^2) # g(x)
## min(y)


Graph 1: just to get an idea of the function

plot(x[order(x)],y[order(x)],
type="l",lwd=2, xlim=c(-4,4),
ylab="y",xlab="x",
main="a. y=g(x)and density of x")
par(new=T)
fx<-exp(-0.5*(x^2))/sqrt(2*pi)
fx<-dnorm(x)
plot(x[order(x)],fx[order(x)],yaxt="n",xaxt="n",xlab="",ylab="",type="l",lty=2,col="grey")
axis(4)
mtext(side=4,"Density",line=2)
legend("topleft",c("y", "x density"),
col=c("black","grey"), lty=1:2, lwd=c(1,2), bty="n")


PDF via method change of variables:

g1.c<-(-b+sqrt((b^2)-4*a*(c-y)))/(2*a)
g2.c<-(-b-sqrt((b^2)-4*a*(c-y)))/(2*a)
g1.prime.c<-abs(1/sqrt((b^2)-4*a*(c-y)))
fy<-dnorm(g1.c)*abs(g1.prime.c)+
dnorm(g2.c)*abs(g1.prime.c)
min(y)
d<-c+(-(b^2)/(4*a))
plot(y,fy,type="l",lwd=2,ylab="density of y",xlab="y", ylim=c(0,0.015),
main="y=80+30x+6x^2")
lines(c(44.4,44.4),c(-1,0.01),lty=2)
lines(c(d,d),c(-1,max(fy)),lty=2,col="red")
legend("topright", c("d=42.5","d+e=44.4"),lty=2,col=c("red","black"))


See how it shoots up?? Results are much worse if you choose larger values for the 'a' parameter.

PDF via the CDF method:

d<-c+(-(b^2)/(4*a))
first<- 1/(2*sqrt(a)*sqrt(y-d))
in_a1<-sqrt(y-d)/sqrt(a)
in_a2<--sqrt(y-d)/sqrt(a)
in_b<-b/(2*a)
A<-in_a1-in_b
B<-in_a2-in_b
d
min(y)
fy_cdf<-first*(dnorm(A)+dnorm(B))
plot(y,fy_cdf,type="l",lwd=2,ylab="density of y",xlab="y", ylim=c(0,0.015),
main="y=80+30x+6x^2")
lines(c(44.4,44.4),c(-1,0.01),lty=2)
lines(c(d,d),c(-1,max(fy)),lty=2,col="red")
legend("topright", c("d=42.5","d+e=44.4"),lty=2,col=c("red","black"))


Note that the results are the same whatever methods is used to derive the pdf:

# library("miscTools")
# compPlot(fy_cdf,fy)
diff<-fy_cdf-fy
summary(abs(diff)) # these are minor rounding errors, I have no issue with that.


An expediment would be to compute either the experimental density or the histogram with:

hist(y)
plot(density(y))


but I would rather have the exact density as I need to compute several statistics from it, notably the probability mass and partial lower mean of some target values. And the derivative of all these with respect to other variables I will interact with X. And well, it is also just for the sack of doing in R what I am writing in the paper...

• For those of us who want to help but don't use R, could you post the graphs you're seeing? – barrycarter Nov 8 '14 at 16:44
• The PDF should be infinite at $c - b^2/(4a)$ and undefined for any smaller values. – whuber Nov 8 '14 at 20:56
• Ok, I thought the PDF was undefined only if y is smaller than $d=c-b^2/(4a)$. Anyway, if I redo the plot with only values for y greater than this quantity, the spike is still there as it starts at 44.4 instead of 42.5. – Gravier Nov 8 '14 at 21:05
• Assuming $x$ has a standard Normal distribution, the correct PDF is $$f(t)=\frac{1}{\sqrt{2 \pi d}} e^{-\frac{d}{2 a^2}} \left(e^{\frac{\left(b-\sqrt{d}\right)^2}{8 a^2}}+e^{\frac{\left(b+\sqrt{d}\right)^2}{8 a^2}}\right) e^{-\frac{c-t}{a}}$$ (for all legitimate values of $t$; otherwise it is $0$) where $d = b^2 - 4a(c-t)$. Is this what you have implemented? – whuber Nov 8 '14 at 21:09
• I am going to try it. No, I implemented this: – Gravier Nov 8 '14 at 21:10

I will treat the case $a > 0$. The other cases can be handled with the same methods.

Let $\Phi$ be the CDF of the standard Normal distribution governing the random variable $X$. Because the graph of $y = ax^2 + bx + c$ is a parabola, the chance that $Y = aX^2 + bX + c$ is less than or equal to $t$ is the chance that $x$ lies between the two roots of $ax^2 + bx + c - t$ (which are easily computed with a quadratic formula). Calling them $x_{+}(t)$ and $x_{-}(t)$, this is given by $\Phi(x_{+}(t)) - \Phi(x_{-}(t))$. Differentiation with respect to $t$ gives the PDF:

$$\phi(t) = \frac{1}{\sqrt{2 \pi d}} e^{-\frac{d}{2 a^2}} \left(e^{\frac{\left(b-\sqrt{d}\right)^2}{8 a^2}}+e^{\frac{\left(b+\sqrt{d}\right)^2}{8 a^2}}\right) e^{-\frac{c-t}{a}}$$

where $d=b^2 - 4a(c-t)$ is the discriminant of the quadratic form and $t$ is not less than the minimum of $y$, equal to $c - b^2/4a.$ Note the factor of $\sqrt{d}$ in the denominator: as $d$ approaches $0$ from above (that is, $c-t$ approaches $b^2/4a$ from higher values of $t$), the exponentials stay finite but the entire PDF has to diverge. Necessarily this PDF will have a vertical asymptote at $t = c - b^2/4a$.

Plots of a typical CDF and PDF (for $a=1, b=-2, c=0$} are:

Because these plots coincide with the (obviously correct) plot of a $\chi^2(1)$ distribution when $b=0$--the two plots just get translated one unit to the right--they likely are correct. Let's look at their counterparts for the example $a=6, b=30, c=80$ in the question:

The PDF looks just like that shown in the question, suggesting the R code actually is correct.

The reason the PDF shoots up at the left is made clear by the CDF: the quadratic behavior of the parabola near its vertex at $-30/12$, where it has a height of $85/2$, focuses a great deal of probability just to the right of $t=85/2$. The infinite density reflects the verticality of the CDF at that location. (The verticality is clear in the first plot; in the second it is apparent only when zoomed in closely near its left endpoint.) Values of $Y$ smaller than $85/2$ cannot be achieved at all, whence there is no probability associated with $t\lt 85/2$.

• That's great, I begin to uderstand now. However, it runs against my intuition: does the infinite density imply that the corresponding $y$ value is the most likely one? Intuitively, it should not. And how comes the PDF sum to 1 if for a given interval the probability mass equals infinity? Is there a rule to identify where to truncate the PDF near the vertex, i.e. the upper bound of this interval? – Gravier Nov 9 '14 at 7:43
• By the upper bound of this interval, I mean where the PDF starts to shoot up. I am interested in recovering several statistics from these statistical distributions: central moments, probabiliy mass below a target level, distance between the lower partial mean and a target level, etc. Does the behaviour of $\chi^2_1 (d^2)$ at the vertex of $y=g(x)$ requires such computation to be done on the CDF rather the PDF? – Gravier Nov 9 '14 at 7:51
• The probability mass never exceeds $1$. The probability density does exceed $1$, but it is never infinite, even though it diverges. This phenomenon happens even with some benign, well-known distributions, such as any Beta$(\alpha,\beta)$ distribution where either parameter is less than $1$. The sharp bimodal shape of your distribution indicates that the simple descriptive statistics you are thinking of will not do a good job at characterizing the distribution. But if you neglect the spike by truncating the distribution, you will potentially ignore a great deal of the probability. – whuber Nov 9 '14 at 17:04