[R] threshold distribution
Dear ALL I have a list of data below 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 all of them are ninty. Nowaday, i try to find a distribution to fit those data. Firstly, I normalize the data, i.e.. (x-mean(X))/(sd(X)) i utilize the SAS to fit my data. Then i obtain the result below ##- Parameters for Lognormal Distribution Parameter Symbol Estimate Threshold Theta -1.51062 Scale Zeta 0.237283 Shape Sigma 0.593481 Mean 0.001321 Std Dev 0.982435 ##--- however, i confuse about the threshold parameter.. How to get it? Does it be able to be calculated by R? Thanks a lot.. __ 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.
Re: [R] threshold distribution
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') meanlogsdlog -0.134806360.19118861 ( 0.03061468) ( 0.02164785) > lines(dlnorm(x, -.1348, .1911), col='red') On Sat, Apr 4, 2009 at 8:30 AM, Abelian wrote: > Dear ALL > I have a list of data below > 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 > all of them are ninty. > Nowaday, i try to find a distribution to fit those data. > Firstly, I normalize the data, i.e.. (x-mean(X))/(sd(X)) > i utilize the SAS to fit my data. Then i obtain the result below > ##- > Parameters for Lognormal > Distribution > > Parameter Symbol > Estimate > > Threshold Theta > -1.51062 > Scale > Zeta 0.237283 > Shape Sigma > 0.593481 > > Mean > 0.001321 > Std > Dev 0.982435 > ##--- > however, i confuse about the threshold parameter.. > How to get it? Does it be able to be calculated by R? > Thanks a lot.. > > __ > 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. > -- Jim Holtman Cincinnati, OH +1 513 646 9390 What is the problem that you are trying to solve? __ 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.
Re: [R] threshold distribution
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')
> meanlogsdlog
> -0.134806360.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')
meanlogsdlog
-0.134806360.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()
adj.pdf
Description: Adobe PDF document
__
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.
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
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.
Re: [R] 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.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
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,
Re: [R] threshold distribution
In my haste I forgot to take the logs of the likelihoods. How embarrassing. The conclusion is unchanged - the data support a non-zero threshold. Here's the corrected code: x.threshold <- 0 loglikelihood <- log(1/(x-x.threshold) * dnorm(log(x - x.threshold), mean(log(x - x.threshold)), sd(log(x - x.threshold sum.loglikelihood <- sum(loglikelihood) print(sum.loglikelihood) x.threshold <- 0.63 loglikelihood.63 <- log(1/(x-x.threshold) * dnorm(log(x - x.threshold), mean(log(x - x.threshold)), sd(log(x - x.threshold sum.loglikelihood.63 <- sum(loglikelihood.63) print(sum.loglikelihood.63) x.minus.x.threshold <- x - x.threshold plot(loglikelihood.63 ~ x.minus.x.threshold, xlim = c(0, 2), ylim =c(-10, 2)) points(x, loglikelihood, col="red") 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 Charles Annis, P.E. Sent: Sunday, April 05, 2009 10:39 AM To: 'Mike Lawrence'; t...@centroin.com.br Cc: 'r-help' Subject: Re: [R] 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.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 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 sd
Re: [R] threshold distribution
Abelian wrote:
Dear ALL
I have a list of data below
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
all of them are ninty.
Nowaday, i try to find a distribution to fit those data.
Firstly, I normalize the data, i.e.. (x-mean(X))/(sd(X))
i utilize the SAS to fit my data. Then i obtain the result below
##-
Parameters for Lognormal
Distribution
Parameter Symbol
Estimate
Threshold Theta
-1.51062
Scale
Zeta 0.237283
Shape Sigma
0.593481
Mean
0.001321
Std
Dev 0.982435
##---
however, i confuse about the threshold parameter..
How to get it? Does it be able to be calculated by R?
Function pelln3 in package lmom will estimate all 3 parameters
of the 3-parameter lognormal distribution, including the threshold ...
> x <- scan(textConnection("
+ 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
+ "))
Read 39 items
>
> y <- (x-mean(x))/sd(x)
>
> library(lmom)
> pelln3(samlmu(y))
zeta mu sigma
-1.5362134 0.2554631 0.5896735
J. R. M. Hosking
__
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.
