Skip to content

trying to plot function using curve

3 messages · Faheem Mitha, Uwe Ligges, Thomas W Blackwell

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)
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)
      }
  }
Uwe Ligges
#
Faheem Mitha wrote:
[Code snipped]

Some comments:
- Indeed, I'd use curve(), but of course you can do it "manually"
   with plot(.., type="l") as well.
- In your usage of curve() you rely heavily on R's scoping rules
   (should work, from my point of view, but always looks a bit
   dangerous):
       densityfn <- function(x) density(x,theta,pos,len) # No defaults!
       # Now curve() grabs theta, pos, len from any environment before:
       curve(densityfn)
- I would not use that many functions as in your code,
   but that's a matter of taste.
- It's not that amusing to debug all that code - try it yourself,
   R has many nice debugging tools (debug(), browser(), etc.).

Uwe Ligges
Thomas W Blackwell
# 
Faheem  -

I observe that function  unnorm1() explicitly returns a vector of
length one, constructed from the first and second elements of its
argument.  Possibly, it was intended to use the first and second
columns of its argument, in which case the subscripts should read
 y[ ,1]  and  y[ ,2].  Uwe Ligges' comment applies.

HTH  -  tom blackwell  -  u michigan medical school  -  ann arbor  -
On Tue, 15 Apr 2003, Faheem Mitha wrote: