Advertisement
Guest User

Untitled

a guest
Feb 20th, 2019
80
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 2.86 KB | None | 0 0
  1. qqnorm(log(y1))
  2. qqnorm(log(y2))
  3. qex <- function(x) qexp((rank(x)-.375)/(length(x)+.25))
  4. plot(qex(y1),log(y1))
  5. plot(qex(y2),log(y2))
  6.  
  7. rm(list=ls())
  8.  
  9. set.seed(123)
  10. x = rlnorm(100,0,1)
  11.  
  12. hist(x)
  13.  
  14. # Loglikelihood and AIC for lognormal model
  15.  
  16. ll1 = function(param){
  17. if(param[2]>0) return(-sum(dlnorm(x,param[1],param[2],log=T)))
  18. else return(Inf)
  19. }
  20.  
  21. AIC1 = 2*optim(c(0,1),ll1)$value + 2*2
  22.  
  23. # Loglikelihood and AIC for Pareto model
  24.  
  25. dpareto=function(x, shape=1, location=1) shape * location^shape / x^(shape + 1)
  26.  
  27. ll2 = function(param){
  28. if(param[1]>0 & min(x)> param[2]) return(-sum(log(dpareto(x,param[1],param[2]))))
  29. else return(Inf)
  30. }
  31.  
  32.  
  33. AIC2 = 2*optim(c(1,0.01),ll2)$value + 2*2
  34.  
  35. # Comparison using AIC, which in this case favours the lognormal model.
  36.  
  37. c(AIC1,AIC2)
  38.  
  39. require(MASS)
  40. hist(x, freq=F)
  41. fit<-fitdistr(x,"log-normal")$estimate
  42. lines(dlnorm(0:max(x),fit[1],fit[2]), lwd=3)
  43.  
  44.  
  45. > fit
  46. meanlog sdlog
  47. 3.8181643 0.1871289
  48.  
  49.  
  50.  
  51. > dput(x)
  52. c(52.6866903145324, 39.7511298620398, 50.0577071855833, 33.8671245370402,
  53. 51.6325665911116, 41.1745418750494, 48.4259060939127, 67.0893697776377,
  54. 35.5355051232044, 44.6197404834786, 40.5620805256951, 39.4265590077884,
  55. 36.0718655240496, 56.0205581625823, 52.8039852992611, 46.2069383488226,
  56. 36.7324212941395, 44.7998046213554, 47.9727885542368, 36.3400338997286,
  57. 32.7514839453244, 50.6878893947656, 53.3756089181472, 39.4769689441593,
  58. 38.5432770167907, 62.350999487007, 44.5140171935881, 47.4026606915147,
  59. 57.3723511479393, 64.4041641945078, 51.2286815562554, 60.4921839777139,
  60. 71.6127652225805, 40.6395409719693, 48.681036613906, 52.3489622656967,
  61. 46.6219563536878, 55.6136160469819, 62.3003761050482, 42.7865905767138,
  62. 50.2413659137295, 45.6327941365187, 46.5621907725798, 48.9734785224035,
  63. 40.4828649022511, 59.4982559591637, 42.9450436744074, 66.8393386407167,
  64. 40.7248473206552, 45.9114242834839, 34.2671010054407, 45.7569869970351,
  65. 50.4358523486278, 44.7445606782492, 44.4173298921541, 41.7506552050873,
  66. 34.5657344132409, 47.7099864540652, 38.1680974794929, 42.2126680994737,
  67. 35.690599714042, 37.6748157160789, 35.0840798650981, 41.4775827114607,
  68. 36.6503753230464, 42.7539062488003, 39.2210050689652, 45.9364763482558,
  69. 35.3687017955285, 62.8299659875044, 38.1532612008011, 39.9183076516292,
  70. 59.0662388169057, 47.9032427690417, 42.4419580084314, 45.785859495192,
  71. 59.5254284342724, 47.9161476636566, 32.6868959277799, 30.1039453246766,
  72. 37.7606323857655, 35.754797368422, 35.5239777126187, 43.7874313667592,
  73. 53.0328404605954, 37.4550326357314, 42.7226751172495, 44.898430515261,
  74. 59.7229655935187, 41.0701258705001, 42.1672231656919, 60.9632847841197,
  75. 60.3690132883734, 45.6469334940722, 39.8300067022836, 51.8185235060234,
  76. 44.908828102875, 50.8200011497451, 53.7945569828737, 65.0432670527801,
  77. 49.0306734716282, 35.9442821219144, 46.8133296904456, 43.7514416949611,
  78. 43.7348972849838, 57.592040060118, 48.7913517211383, 38.5555058596449
  79. )
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement