R help - Mar 2014 - Help resolving issues with generating a chi-squared density plot from scratch

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]]

I am far too lazy to go through all your rather complicated code, but my 
basic impression is that you are re-inventing quite a few wheels.

Just to get a plot of the density function you can simply do, e.g.:

	plot(function(x){dchisq(x,1)},0,qchisq(0.999,1))

I can't see why you are fooling about with the ylim values; just let the 
plot() function choose ylim.

As to why you can't get things to work when df=1 ... well don't try to 
set the upper y limit equal to infinity! You have a finite plotting 
region, after all.

I have no idea what you mean by "The y-axis is numbered only 
relatively"; this makes no sense at all.  What *do* you want?  The 
y-axis labelling that you get will be/is the appropriate labelling.

The arrows will go where you tell them to go, so if they are "going 
diagonal" you are telling them to go diagonal.

cheers,

Rolf Turner

On 10/03/14 13:30, Levi Robinson 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")
>    }
>
> }

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.