R help - Nov 2005 - help with legacy R code

Hi there,


Could somebody help me disect this legacy R script I inherited at work, I
have two questions:
1. I've tried to upgrade our R version from 1.6.2 (yeah, I know), to R 2.0,
but some of the lines in this script are not compatible with R 2.0, could
someone help me figure out where the problem is?
2. the jpeg generated (attached) seems to be off on some of the data, is
there a better way of doing this.

Thanks very much in advance!

David


 library(MASS)
 jpeg(filename = "diswrong.jpg", width = 800, height = 600, pointsize
= 12,
quality = 75, bg = "white")

 myfunc <- function(x, mean, sd, nfalse, ntotal, shape, rate) {
 (nfalse*dgamma(x,shape,rate)+(ntotal-nfalse)*dnorm(x,mean,sd))/ntotal
 }

 wrong <- scan("wrongrawdata.txt", list(x=0))
 wrongfit <- fitdistr(wrong$x, "gamma")
 wrongmean <- mean(wrong$x)
 wrongshape <- wrongfit[[1]][1]
 wrongrate <- wrongfit[[1]][2]

 good <- scan("rawdata.txt", list(x=0))
 xmin = 0
 newx = good$x
 xmean = mean(newx)


 xmax = max(newx)+0.15
 goodhist <- hist(newx, br=seq(from=0,to=xmax,by=0.15), probability=T,
col="lightyellow")

 initmean <- (min(newx)+max(newx))/2
 totalx <- length(newx)

 wrongmeanshift <- wrongmean + 0.2
 wrongper <- pgamma(wrongmeanshift, wrongshape, wrongrate)
 nfalseundermean <-
which(abs(newx-wrongmeanshift)==min(abs(newx-wrongmeanshift)))
 initnfalse <- nfalseundermean / wrongper

 fitmean <- -1
 fitsd <- 0
 fitnfalse <- initnfalse
 fitshape <- wrongshape
 fitrate <- wrongrate

 curve((fitnfalse*dgamma(x,fitshape,fitrate))/totalx, add=T,
col="red",
lwd=2)

 breaksllength <- length(goodhist$breaks)
 endi = breaksllength - 1
 binprob = c(1)
 for (i in 1:endi) {
 expnegative <- fitnfalse * (pgamma(goodhist$breaks[i+1],wrongshape,
wrongrate)-pgamma(goodhist$breaks[i],wrongshape, wrongrate))
 if (goodhist$counts[i] == 0)
 binprob[i] = 0
 else
 binprob[i] = (goodhist$counts[i] - expnegative) / goodhist$counts[i]
 }

 result = data.frame(newx)
 prob = c(1)
 for (i in 1:totalx) {
 bini = which ((goodhist$breaks < newx[i]) & (goodhist$breaks >
newx[i]-0.15
))
 if ((binprob[bini] < 0.8) | (newx[i] < wrongmean))
 prob[i] = -1
 else
 prob[i] = binprob[bini]*100
 }

 result = data.frame(result, prob)
 write.table(result, file="probwrong.txt", sep=" ",
row.name=F, col.name=F)
 fitpars = c(fitmean, fitsd, fitnfalse, fitshape, fitrate, totalx)
 result = data.frame(fitpars)
 write.table(result,file="parwrong.txt", sep=" ",
row.name=F, col.name=F)
 dev.off()

David Zhao wrote:
> Hi there,
> 
> 
> Could somebody help me disect this legacy R script I inherited at work, I
> have two questions:
> 1. I've tried to upgrade our R version from 1.6.2 (yeah, I know), to R
2.0,
> but some of the lines in this script are not compatible with R 2.0, could
> someone help me figure out where the problem is?
> 2. the jpeg generated (attached) seems to be off on some of the data, is
> there a better way of doing this.
1a. R 2.0 must be a software I am not familar with, since for the R I 
know such a version has never been released.
1b. We are unable to reproduce the stuff given below. Not even an error 
message is given.
1c. Do you expect anybody has the time to make your own homework, in 
particular on an unreproducible example? There are very convenient 
debugging tools made available for you in R.
2. We do not see where your jpeg produced with your data is "off".


Uwe Ligges

PS: Let me quote
 > PLEASE do read the posting guide! 
http://www.R-project.org/posting-guide.html


> Thanks very much in advance!
> 
> David
> 
> 
>  library(MASS)
>  jpeg(filename = "diswrong.jpg", width = 800, height = 600,
pointsize = 12,
> quality = 75, bg = "white")
> 
>  myfunc <- function(x, mean, sd, nfalse, ntotal, shape, rate) {
>  (nfalse*dgamma(x,shape,rate)+(ntotal-nfalse)*dnorm(x,mean,sd))/ntotal
>  }
> 
>  wrong <- scan("wrongrawdata.txt", list(x=0))
>  wrongfit <- fitdistr(wrong$x, "gamma")
>  wrongmean <- mean(wrong$x)
>  wrongshape <- wrongfit[[1]][1]
>  wrongrate <- wrongfit[[1]][2]
> 
>  good <- scan("rawdata.txt", list(x=0))
>  xmin = 0
>  newx = good$x
>  xmean = mean(newx)
> 
> 
>  xmax = max(newx)+0.15
>  goodhist <- hist(newx, br=seq(from=0,to=xmax,by=0.15), probability=T,
> col="lightyellow")
> 
>  initmean <- (min(newx)+max(newx))/2
>  totalx <- length(newx)
> 
>  wrongmeanshift <- wrongmean + 0.2
>  wrongper <- pgamma(wrongmeanshift, wrongshape, wrongrate)
>  nfalseundermean <-
> which(abs(newx-wrongmeanshift)==min(abs(newx-wrongmeanshift)))
>  initnfalse <- nfalseundermean / wrongper
> 
>  fitmean <- -1
>  fitsd <- 0
>  fitnfalse <- initnfalse
>  fitshape <- wrongshape
>  fitrate <- wrongrate
> 
>  curve((fitnfalse*dgamma(x,fitshape,fitrate))/totalx, add=T,
col="red",
> lwd=2)
> 
>  breaksllength <- length(goodhist$breaks)
>  endi = breaksllength - 1
>  binprob = c(1)
>  for (i in 1:endi) {
>  expnegative <- fitnfalse * (pgamma(goodhist$breaks[i+1],wrongshape,
> wrongrate)-pgamma(goodhist$breaks[i],wrongshape, wrongrate))
>  if (goodhist$counts[i] == 0)
>  binprob[i] = 0
>  else
>  binprob[i] = (goodhist$counts[i] - expnegative) / goodhist$counts[i]
>  }
> 
>  result = data.frame(newx)
>  prob = c(1)
>  for (i in 1:totalx) {
>  bini = which ((goodhist$breaks < newx[i]) & (goodhist$breaks >
newx[i]-0.15
> ))
>  if ((binprob[bini] < 0.8) | (newx[i] < wrongmean))
>  prob[i] = -1
>  else
>  prob[i] = binprob[bini]*100
>  }
> 
>  result = data.frame(result, prob)
>  write.table(result, file="probwrong.txt", sep=" ",
row.name=F, col.name=F)
>  fitpars = c(fitmean, fitsd, fitnfalse, fitshape, fitrate, totalx)
>  result = data.frame(fitpars)
>  write.table(result,file="parwrong.txt", sep=" ",
row.name=F, col.name=F)
>  dev.off()
> 
> 
> ------------------------------------------------------------------------
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html

1) There is no 'R 2.0'.  What version did you mean?

2) We cannot reproduce your script (no data files), and JPEGs are not 
allowed on R-help: see http://www.r-project.org/mail.html.  (PNGs are, 
though).

3) You have given us no indication of what the problems are nor in which 
lines.

Please give us more usable information.


On Wed, 9 Nov 2005, David Zhao wrote:
> Could somebody help me disect this legacy R script I inherited at work, I
> have two questions:
> 1. I've tried to upgrade our R version from 1.6.2 (yeah, I know), to R
2.0,
> but some of the lines in this script are not compatible with R 2.0, could
> someone help me figure out where the problem is?
> 2. the jpeg generated (attached) seems to be off on some of the data, is
> there a better way of doing this.
> library(MASS)
> jpeg(filename = "diswrong.jpg", width = 800, height = 600,
pointsize = 12,
> quality = 75, bg = "white")
>
> myfunc <- function(x, mean, sd, nfalse, ntotal, shape, rate) {
> (nfalse*dgamma(x,shape,rate)+(ntotal-nfalse)*dnorm(x,mean,sd))/ntotal
> }
>
> wrong <- scan("wrongrawdata.txt", list(x=0))
...


-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595