# Big rounding error in R or miscoding? Pdf of a quadratic of function of N(0,1)

#### Gravier1Billion

##### New Member
Dear all,
I am trying to plot 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<d, where d=c-(b^2)/4a (see attachment for the details).
Strangely enough, when computing it with R, the pdf shoots up at y>d+e, where is quite large.
Is there an error in my codes (attached) or is it a rounding error?
In the latter case, how to address it?
Xav
ps. not sure the attachment went through, so please find below the code:

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

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

# method by CDF
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"))

# # 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)) # likely to to rounding error, but I have
# no issue with that.

#### Gravier1Billion

##### New Member
With the attachement this time, little internet access problem 