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?

Thanks in advance!

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.