I wrote the code to graph a chi-squared density function, shade the
percentile, and point to the CV, but it has a few issues I can't seem to
resolve
1. It won't work at all for DF = 1 due to ylim going to infinity, but I
haven't been able to resolve this still after hours of trying.
2) The y-axis is numbered only relatively; I'd prefer they were the actual
prob densities, but again, I fiddled with a few things, but it just
wouldn't get me what I wanted.
3) For low degrees of freedom and higher percentiles, the arrow pointing to
CV seems to end up going diagonal instead of straight down
Here's the code below and here's the URL for R fiddle for the code which
might make it easier to fix:
http://www.r-fiddle.org/#/fiddle?id=ChFi0dyJ&version=4
chi.dens = function(x = NULL, df = NULL, cv = NULL) {
# x = percentile/quantile
# df = degrees of freedom
# cv = critical value
if(x>1 ||x<0) stop("Percentile must be between 0 and 1") #
Error-handling
qend = qchisq(x, df)
perc = x
qt0=qchisq(0.5, df) # Defining for arrows later
dy0=dchisq(0.45, df) # Defining for arrows later
xrange = qchisq (0.999, df)
x = seq(0, xrange)
y = dchisq(x, df)
yheight = max(dchisq(x, df))
# Creating plot
plot(x,y,type = "l", ylim=c(0,yheight),axes=FALSE)
title( "Chi-squared Density Curve with")
mtext(bquote(paste("DF = ", .(df), " and Percentile = ", .(perc))),
side = 3, line = 0) # Input information
# Shading left tail-region
qt = signif(qend,5)
x0=seq(0, qt)
y0=dchisq(x0, df)
polygon(c(0, x0, qt), c(0, y0, 0), col = "lightblue")
xtks=signif(seq(0,xrange,length=10),3)
axis(1, pos=0, at=xtks, labels=xtks)
y.unit=max(dchisq(x, df))/5
y.pos=seq(0,5*y.unit, length=5)
y.lab=c(0, 1, 2, 3, 4) # Y axis numbers only reflect relative
densities to each other.
axis(2,pos=0, at=y.pos, labels=y.lab) # set up y-axis
# "Normal" CV less than the 99.9 Percentile:
if(qt <= xrange){
if(length(cv)==0) text(0.75*xrange,3*y.unit, bquote(paste("CV = ",
.(qend))), cex=1.2, col = "red",adj=0.5)
#
if(length(perc)==1) text(0.35*xrange,3*y.unit, paste("Area = ",
perc), cex=1.2, col="darkblue", adj=0.5)
if(perc>0.5){
arrows(0.35*xrange, 2.5*y.unit, qt0, 0.3*dy0, length=0.1,
angle=20) # pointing to the shaded region
}
if(perc<=0.5){
arrows(0.35*xrange, 2.5*y.unit, qt/2, 0.3*dy0, length=0.1,
angle=20) # pointing to the shaded region
}
arrows(0.75*xrange, 2.5*y.unit, qt, 0, length=0.1, angle=20)
points(qt,0,pch=17, col="red")
}
# When CV is greater than the 99.9 Percentile:
if(qt > xrange){
print("CV may be too far to the right to be shown on graph due to the
high percentile")
if(length(cv)==0) text(0.75*xrange,3*y.unit, bquote(paste("CV = ",
.(qend))), cex=1.2, col = "red",adj=0.5)
if(length(perc)==1) text(0.35*xrange,3*y.unit, paste("Area = ",
perc), cex=1.2, col="darkblue", adj=0.5)
if(perc>0.5){
arrows(0.35*xrange, 2.5*y.unit, qt0, 0.3*dy0, length=0.1,
angle=20) # pointing to the shaded region
}
if(perc<=0.5){
arrows(0.35*xrange, 2.5*y.unit, qt/2, 0.3*dy0, length=0.1,
angle=20) # pointing to the shaded region
}
arrows(0.75*xrange, 2.5*y.unit, qt, 0, length=0.1, angle=20)
points(qt,0,pch=17, col="red")
}
}
