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) } }
Faheem Mitha wrote:> 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. >[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
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:> 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. > . . . > > 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. > > *********************************************************************** > . . . > 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))) > } > . . . >