[R] help with legacy R code

[David Zhao]

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()