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() > > > ------------------------------------------------------------------------ > > ______________________________________________ > [email protected] mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html ______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
