3
$\begingroup$

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

enter image description here

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

enter image description here

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...

$\endgroup$
  • $\begingroup$ For those of us who want to help but don't use R, could you post the graphs you're seeing? $\endgroup$ – barrycarter Nov 8 '14 at 16:44
  • $\begingroup$ The PDF should be infinite at $c - b^2/(4a)$ and undefined for any smaller values. $\endgroup$ – whuber Nov 8 '14 at 20:56
  • $\begingroup$ 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. $\endgroup$ – Gravier Nov 8 '14 at 21:05
  • $\begingroup$ 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? $\endgroup$ – whuber Nov 8 '14 at 21:09
  • $\begingroup$ I am going to try it. No, I implemented this: $\endgroup$ – Gravier Nov 8 '14 at 21:10
3
$\begingroup$

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:

Figure

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:

Figure 2

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$.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Gravier Nov 9 '14 at 7:43
  • $\begingroup$ 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? $\endgroup$ – Gravier Nov 9 '14 at 7:51
  • 1
    $\begingroup$ 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. $\endgroup$ – whuber Nov 9 '14 at 17:04

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.