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") } } [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
