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.
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 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.
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 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
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 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
[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 abelian1...@gmail.com 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.
