Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- qqnorm(log(y1))
- qqnorm(log(y2))
- qex <- function(x) qexp((rank(x)-.375)/(length(x)+.25))
- plot(qex(y1),log(y1))
- plot(qex(y2),log(y2))
- rm(list=ls())
- set.seed(123)
- x = rlnorm(100,0,1)
- hist(x)
- # Loglikelihood and AIC for lognormal model
- ll1 = function(param){
- if(param[2]>0) return(-sum(dlnorm(x,param[1],param[2],log=T)))
- else return(Inf)
- }
- AIC1 = 2*optim(c(0,1),ll1)$value + 2*2
- # Loglikelihood and AIC for Pareto model
- dpareto=function(x, shape=1, location=1) shape * location^shape / x^(shape + 1)
- ll2 = function(param){
- if(param[1]>0 & min(x)> param[2]) return(-sum(log(dpareto(x,param[1],param[2]))))
- else return(Inf)
- }
- AIC2 = 2*optim(c(1,0.01),ll2)$value + 2*2
- # Comparison using AIC, which in this case favours the lognormal model.
- c(AIC1,AIC2)
- require(MASS)
- hist(x, freq=F)
- fit<-fitdistr(x,"log-normal")$estimate
- lines(dlnorm(0:max(x),fit[1],fit[2]), lwd=3)
- > fit
- meanlog sdlog
- 3.8181643 0.1871289
- > dput(x)
- c(52.6866903145324, 39.7511298620398, 50.0577071855833, 33.8671245370402,
- 51.6325665911116, 41.1745418750494, 48.4259060939127, 67.0893697776377,
- 35.5355051232044, 44.6197404834786, 40.5620805256951, 39.4265590077884,
- 36.0718655240496, 56.0205581625823, 52.8039852992611, 46.2069383488226,
- 36.7324212941395, 44.7998046213554, 47.9727885542368, 36.3400338997286,
- 32.7514839453244, 50.6878893947656, 53.3756089181472, 39.4769689441593,
- 38.5432770167907, 62.350999487007, 44.5140171935881, 47.4026606915147,
- 57.3723511479393, 64.4041641945078, 51.2286815562554, 60.4921839777139,
- 71.6127652225805, 40.6395409719693, 48.681036613906, 52.3489622656967,
- 46.6219563536878, 55.6136160469819, 62.3003761050482, 42.7865905767138,
- 50.2413659137295, 45.6327941365187, 46.5621907725798, 48.9734785224035,
- 40.4828649022511, 59.4982559591637, 42.9450436744074, 66.8393386407167,
- 40.7248473206552, 45.9114242834839, 34.2671010054407, 45.7569869970351,
- 50.4358523486278, 44.7445606782492, 44.4173298921541, 41.7506552050873,
- 34.5657344132409, 47.7099864540652, 38.1680974794929, 42.2126680994737,
- 35.690599714042, 37.6748157160789, 35.0840798650981, 41.4775827114607,
- 36.6503753230464, 42.7539062488003, 39.2210050689652, 45.9364763482558,
- 35.3687017955285, 62.8299659875044, 38.1532612008011, 39.9183076516292,
- 59.0662388169057, 47.9032427690417, 42.4419580084314, 45.785859495192,
- 59.5254284342724, 47.9161476636566, 32.6868959277799, 30.1039453246766,
- 37.7606323857655, 35.754797368422, 35.5239777126187, 43.7874313667592,
- 53.0328404605954, 37.4550326357314, 42.7226751172495, 44.898430515261,
- 59.7229655935187, 41.0701258705001, 42.1672231656919, 60.9632847841197,
- 60.3690132883734, 45.6469334940722, 39.8300067022836, 51.8185235060234,
- 44.908828102875, 50.8200011497451, 53.7945569828737, 65.0432670527801,
- 49.0306734716282, 35.9442821219144, 46.8133296904456, 43.7514416949611,
- 43.7348972849838, 57.592040060118, 48.7913517211383, 38.5555058596449
- )
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement