[R] threshold distribution

Sun Apr 5 16:38:58 CEST 2009

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.