[R] help with legacy R code
David Zhao
wzhao6898 at gmail.com
Thu Nov 10 00:11:12 CET 2005
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()
More information about the R-help
mailing list