threshold distribution
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.Annis at StatisticalEngineering.com phone: 561-352-9699 eFax: 614-455-3265 http://www.StatisticalEngineering.com -----Original Message----- From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Mike Lawrence Sent: Sunday, April 05, 2009 8:13 AM To: tura at 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
<tura at 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 at 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 at 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.