3

My probability distribution function is

f(x;r)=(1+r)/Γ(1/(1+r))exp(-x^(1+r)) , 0<x<Inf, r>-1

I am using following code:

library(stats4)   
x3=c(0.927869895,0.000559193,0.059761785,0.361542986,0.274291999,0.563373589,
     0.565878385,0.126281994,0.46623362,0.796111638,0.809460469,0.28011156,
     0.263443012,0.704551045,0.978333171,0.139962088,0.599336759,0.108009066,
     0.202133384,0.05401611)

nLL <- function(r)-sum(log(1+r)-log(gamma(1/(1+r)))-x3^(1+r))
fit0 <- mle(nLL, start = list(r = 0.1), nobs = NROW(x3))
fit1 <- mle(nLL, start = list(r = 0.1), nobs = NROW(x3),
            method = "Brent", lower = -1, upper = 150)
stopifnot(nobs(fit0) == length(x3))

Now I want to extract value of MLE stored in fit1. I have used following command

fitv = fit1$Coefficients

But I am unable to get the value as it give following error

Error in fit1$Coefficients : $ operator not defined for this S4 class

Now I get the answer 

nLL <- function(r)-sum(log(1+r)-log(gamma(1/(1+r)))-x1^(1+r))
fit0 <- mle(nLL, start = list(r = 0.1), nobs = NROW(x1))
fit1 <- mle(nLL, start = list(r = 0.1), nobs = NROW(x1),
        method = "Brent", lower = -1, upper = 200)
stopifnot(nobs(fit0) == length(x1))
class(fit1)
str(fit1)
fiteg[k] = fit1@coef
Zheyuan Li
  • 71,365
  • 17
  • 180
  • 248
aorstat
  • 31
  • 4
  • What exactly did you try and what didn't work? Are you sure that's a PDF? It doesn't seem to integrate to 1 for `r=.5` for example. Are you trying to solve this numerically or analytically? – MrFlick Jun 05 '15 at 04:03
  • It is a pdf because it is a special case of gamma distribution with c=1,scale parameter, and b=1/(1+r), shape parameter. I am so sorry about typing. I am studying to install MathJax. – aorstat Jun 05 '15 at 05:44
  • Is it really called the "unformal" distribution? – Mike Wise Jun 05 '15 at 06:10
  • I am so sorry. I just prove it that it is a special case of gamma distribution. I will change my tropic. I apologize about my wrongs. – aorstat Jun 05 '15 at 06:29
  • @Sania, Thanks for your help-. – aorstat Jun 05 '15 at 06:44
  • Also, for using S4 classes, see for example [here](http://stackoverflow.com/q/13099780/903061). – Gregor Thomas Jun 05 '15 at 06:45
  • @Gregor, Thank you for your example. Your example is very useful for my question. Now I can get the value. – aorstat Jun 05 '15 at 07:10
  • Have you tried `logLik(fit1)`? Works for me. You can also try examining the structure of the object by `str(fit1)`. First elements (called slots for S4) are accessible via @, lower lever elements through $. – Roman Luštrik Jun 05 '15 at 12:58
  • @Roman Lustrink, Thank you very much. Yes, your answer is correct. – aorstat Jul 01 '15 at 13:00

0 Answers0