[R] 2 plots 1 figure
Dermot MacSweeney
dsweeney at nmrc.ucc.ie
Wed Nov 15 11:16:21 CET 2000
Hi,
I am wondering would it not be easier to fit curves to the histograms and then
plot these curves. These are not as accurate but it is easier to compare the
results. For instance the following function has as an input two measured data
sets, e.g.,
plothist(rnorm(500),rnorm(500))
The function takes these data sets and fits a normal function to each of them
and also scales the function to the dataset size. It then plots the fitted
curves (It is a little messy as I had only intended to use it to fit one data
set at a time but I am sure that someone should should be able to modify it
easily).
function (input1,input2)
{
# Calculate Stats for input1
sigma1<-sd(input1)
mu1<-mean(input1)
xlow1<-(min(input1)) - 2*sigma1
xhigh1<-(max(input1)) + 2*sigma1
# Calculate Stats for input2
sigma2<-sd(input2)
mu2<-mean(input2)
xlow2<-(min(input2)) - 2*sigma2
xhigh2<-(max(input2)) + 2*sigma2
# Calculate scaling for fit1
hist.obj1<-hist(input1,xlim=range(xlow1,xhigh1),plot=F)
k1<-hist.obj1$counts[1]/hist.obj1$intensities[1]
# Calculate scaling for fit2
hist.obj2<-hist(input2,xlim=range(xlow2,xhigh2),plot=F)
k2<-hist.obj2$counts[1]/hist.obj2$intensities[1]
# Calculate fit1
normal.fit1 <- function(x)
k1/(sqrt(2*3.142)*sigma1)*exp(0.5*(-(x-mu1)**2)*sigma1**-2)
# Calculate fit2
normal.fit2 <- function(x)
k2/(sqrt(2*3.142)*sigma2)*exp(0.5*(-(x-mu2)**2)*sigma2**-2)
# Plot measured and fitted data
hist(input1,xlim=range(xlow1,xhigh1),freq=T)
#hist(input2,xlim=range(xlow1,xhigh1),freq=T,col="blue")
curve(normal.fit1,add=T,xlow1,xhigh1,col="red")
curve(normal.fit2,add=T,xlow2,xhigh2,col="blue")
box()
}
**************************************************************
Dermot MacSweeney NMRC,
Email: dsweeney at nmrc.ucc.ie Lee Maltings,
Tel: +353 21 904178 Prospect Row,
Fax: +353 21 270271 Cork,
WWW: http://nmrc.ucc.ie Ireland
**************************************************************
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list