The data suggest a lognormal threshold to me (and perhaps to the originator, based on his title).
###### x <- c(0.80010, 0.72299, 0.69893, 0.99597, 0.89200, 0.69312, 0.73613, 1.13559, 0.85009, 0.85804, 0.73324, 1.04826, 0.84002, 1.76330, 0.71980, 0.89416, 0.89450, 0.98670, 0.83571, 0.73833, 0.66549, 0.93641, 0.80418, 0.95285, 0.76876, 0.82588, 1.09394, 1.00195, 1.14976, 0.80008, 1.11947, 1.09484, 0.81494, 0.68696, 0.82364, 0.84390, 0.71402, 0.80293, 1.02873) # plot the original density x.threshold <- 0 qqnorm(log(x - x.threshold), datax = TRUE) qqline(log(x - x.threshold), datax = TRUE) # plot the lognormal density with a threshold x.threshold <- 0.63 qqnorm(log(x - x.threshold), datax = TRUE) qqline(log(x - x.threshold), datax = TRUE) # compute and plot the associated log likelihoods x.threshold <- 0 loglikelihood <- 1/(x-x.threshold) * dnorm(log(x - x.threshold), mean(log(x - x.threshold)), sd(log(x - x.threshold))) x.threshold <- 0.63 loglikelihood.63 <- 1/(x-x.threshold) * dnorm(log(x - x.threshold), mean(log(x - x.threshold)), sd(log(x - x.threshold))) x.minus.x.threshold <- x - x.threshold plot(x.minus.x.threshold, loglikelihood.63, xlim = c(0, 2), ylim =c(0, 3.5)) points(x, loglikelihood, col="red") # compare loglikelihood ratio with chi-square sum.loglikelihood <- sum(loglikelihood) print(sum.loglikelihood) sum.loglikelihood.63 <- sum(loglikelihood.63) print(sum.loglikelihood.63) log.likelihiood.ratio <- sum.loglikelihood.63 - sum.loglikelihood print(log.likelihiood.ratio) significant.difference <- qchisq(p = 0.95, df = 1)/2 print(significant.difference) ##### Charles Annis, P.E. charles.an...@statisticalengineering.com phone: 561-352-9699 eFax: 614-455-3265 http://www.StatisticalEngineering.com -----Original Message----- From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On Behalf Of Mike Lawrence Sent: Sunday, April 05, 2009 8:13 AM To: t...@centroin.com.br Cc: r-help Subject: Re: [R] threshold distribution You didn't properly specify the x-axis coordinates in your call to lines(): plot(density(x)) lines(x,dlnorm(x, -.1348, .1911), col='red') points(x,dlnorm(x, -.1348, .1911), col='blue') curve(dlnorm(x, -.1348, .1911), col='green',add=T) On Sun, Apr 5, 2009 at 6:40 AM, Bernardo Rangel Tura <t...@centroin.com.br> wrote: > On Sun, 2009-04-05 at 01:13 -0400, jim holtman wrote: >> Here is what I get from using 'fitdistr' in R to fit to a lognormal. >> The resulting density plot from the distribution seems to be a reason >> match to the data. >> >> > x <- scan() >> 1: 0.80010 0.72299 0.69893 0.99597 0.89200 0.69312 0.73613 1.13559 >> 9: 0.85009 0.85804 0.73324 1.04826 0.84002 >> 14: 1.76330 0.71980 0.89416 0.89450 0.98670 0.83571 0.73833 0.66549 >> 22: 0.93641 0.80418 0.95285 0.76876 0.82588 >> 27: 1.09394 1.00195 1.14976 0.80008 1.11947 1.09484 0.81494 0.68696 >> 35: 0.82364 0.84390 0.71402 0.80293 1.02873 >> 40: >> Read 39 items >> > plot(density(x)) >> > library(MASS) >> > fitdistr(x, 'lognormal') >> meanlog sdlog >> -0.13480636 0.19118861 >> ( 0.03061468) ( 0.02164785) >> > lines(dlnorm(x, -.1348, .1911), col='red') > > Hi Jim > > I agree with your solution but my plot result not fine. > I obtain same result. >> fitdistr(x, 'lognormal') > meanlog sdlog > -0.13480636 0.19118861 > ( 0.03061468) ( 0.02164785) > > In plot when I use points (blue) and curve (green) the fit o lognormal > and density(data) is fine but when I use lines (red) the fit is bad (in > attach I send a PDF of my output) > > Do you know why this happen? > > Thanks in advance > > -- > Bernardo Rangel Tura, M.D,MPH,Ph.D > National Institute of Cardiology > Brazil > > P.S. my script is > > x <- scan() > 0.80010 0.72299 0.69893 0.99597 0.89200 0.69312 0.73613 1.13559 > 0.85009 0.85804 0.73324 1.04826 0.84002 > 1.76330 0.71980 0.89416 0.89450 0.98670 0.83571 0.73833 0.66549 > 0.93641 0.80418 0.95285 0.76876 0.82588 > 1.09394 1.00195 1.14976 0.80008 1.11947 1.09484 0.81494 0.68696 > 0.82364 0.84390 0.71402 0.80293 1.02873 > > require(MASS) > fitdistr(x, 'lognormal') > pdf("adj.pdf") > plot(density(x)) > lines(dlnorm(x, -.1348, .1911), col='red') > points(x,dlnorm(x, -.1348, .1911), col='blue') > curve(dlnorm(x, -.1348, .1911), col='green',add=T) > dev.off() > > > > ______________________________________________ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > > -- Mike Lawrence Graduate Student Department of Psychology Dalhousie University Looking to arrange a meeting? Check my public calendar: http://tinyurl.com/mikes-public-calendar ~ Certainty is folly... I think. ~ ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
