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()
______________________________________________ R-help@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