If you set range of x in [0, something ] for for a function that goes to infinity for x=0, then what do you expect?
> On Mar 10, 2014, at 7:30 AM, Levi Robinson <robin...@gmail.com> wrote: > > 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. ______________________________________________ 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.
