[R] trying to plot function using curve

[Faheem Mitha]

Dear People,

I hope someone can help me with this. I have a function (density) which I
am trying to plot using curve. I am calling mg.hist(3,2,1), and getting
the following errors.

Error in xy.coords(x, y, xlabel, ylabel, log) :
	x and y lengths differ
In addition: There were 50 or more warnings (use warnings() to see the
first 50)
> warnings()
1: longer object length
	is not a multiple of shorter object length in: x * y

and more of the same.

I'm not sure where the problem is. I assume it has something to do with me
incorrectly vectorising my functions(s). However, it seems to me that
density() is vectorised.

By the way, if anyone would like to suggest a better way to plot density
plots, please let me know. curve() was suggested to me, which is why I am
using it. Actually, I want to plot a density curve on top of a histogram.
I was thinking of trying to use trellis graphics for this, but I'm not
sure it has any advantages for this purpose.

                                                     Faheem.

***********************************************************************

mg.hist <- function(len,theta,pos)
{
  postscript(file="plot.ps", horizontal = FALSE, onefile = FALSE, paper
             = "special", width=6, height=4)

  densityfn <- function(x)
    {
      density(x,theta,pos,len)
    }
  curve(densityfn)
  dev.off()
 }

density <- function(x,theta,pos,len)
  {
    if((len !=2) && (len !=3))
      stop("length must be 2 or 3")
    if((pos < 1) || (pos > len))
      stop("pos must be between 0 and len-1")

    if(len==2)
      {
        if((pos==1) || (pos==2) )
          {
            unnorm <- function(x)
              {
                ifelse(x==0,2*theta,(exp(theta*x)-exp(-theta*x))/x)
              }
            return(unnorm(x)/
                   integrate(unnorm,lower = -theta, upper = theta,
                             subdivisions=1000)$value)
          }

      }

    if(len==3)
      {
        if((pos==1) || (pos==3))
          {
            unnorm1 <- function(y) ##here y is 2-dim
              {
                ifelse(y[2] == 0,2*theta,
                       (1/y[2])*(exp(y[1]*y[2] + y[2]*theta)
                                 - exp(y[1]*y[2] - y[2]*theta)))
              }
            unnorm2 <- function(y)
              {
                ## return(unnorm1(c(x,y)))  ##here y is 1-dim
                ifelse(y == 0,2*theta,
                       (1/y)*(exp(x*y + y*theta)
                              - exp(x*y - y*theta)))
              }
          }
        return(integrate(unnorm2, lower=-theta, upper=theta,

subdivisions=1000)$value/adapt(ndim=2,lower=c(-theta, -theta),
                           upper=c(theta, theta), functn=unnorm1)$value)
      }

    if(pos==2)
      {
        unnorm <- function(x)
          {
            ifelse(x==0,4,(exp(2*theta*x)+exp(-2*theta*x) - 2)/(x^2))
          }
        return(integrate(unnorm,lower = -theta, upper = x,
                         subdivisions=1000)$value/
               integrate(unnorm,lower = -theta, upper = theta,
                         subdivisions=1000)$value)
      }
  }