我想计算由以下定义的随机变量 y 的 pdf：
y=c+b*x+a*x^2
pdf 是非中心卡方分布。对于 a>0，如果 y 小于 d，则它应该等于 0，其中 d=c-(b^2)/4a。
奇怪的是，当用 R 计算它时，pdf 在 y>d+e 处上升，其中 e 非常大。
我的代码（如下）有错误还是舍入错误？在后一种情况下，如何解决？
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)
图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")
PDF通过变量的方法更改
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"))
看看它是怎么弹起来的？？
通过 CDF 获得 PDF
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"))
无论使用何种方法导出 pdf，结果都是相同的
# library("miscTools")
# compPlot(fy_cdf,fy)
# diff<-fy_cdf-fy
# summary(abs(diff)) # these are minor rounding errors,
# no issue with that.
