[R] help with legacy R code

[Uwe Ligges]

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