[R] Help resolving issues with generating a chi-squared density plot from scratch

[David Winsemius]

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 <robinsle at 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 at 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.