Here are the packages needed, the data, example with the same codes and results and then when i tried with different model i got error
library(stats)
library(stats4)
library(base)
library(survival)
library(lsei)
library(npsurv)
library(splines)
library(VGAM)
library(flexsurv)
library(rmutil)
library(MASS)
library(methods)
library(bbmle)
library(gsl)
library(hypergeo)
library (methods)
library(miscTools)
library(maxLik)
library(EstimationTools)
x=read.table("ten.csv")
> x
1.312
1.997
2.27
2.478
2.642
2.848
3.585
1.314
2.006
2.272
2.49
2.648
2.88
3.585
1.479
2.027
2.274
2.511
2.684
2.954
1.552
2.055
2.301
2.514
2.697
3.012
1.7
2.063
2.301
2.535
2.726
3.067
1.803
2.098
2.359
2.554
2.77
3.084
1.861
2.14
2.382
2.566
2.773
3.09
1.865
2.179
2.382
2.57
2.8
3.096
1.944
2.224
2.426
2.586
2.809
3.128
1.958
2.24
2.434
2.629
2.818
3.233
1.966
2.253
2.435
2.633
2.821
3.433
Here is example with the results:
dpl=function(x, mu, sigma, log= FALSE){
if(any(x<0))
stop(paste("x must be positive", "\n", ""))
if (any(mu<=0))
stop(paste("mu must be positive","\n", ""))
if (any(sigma<=0))
stop(paste("sigma must be positive", "\n", ""))
loglik=log(mu)+2*log(sigma)-log(sigma+1)+log(1+(x^mu))+(mu-1)*log(x)-sigma*(x^mu)
if(log==FALSE)
density=exp(loglik)
else density=loglik
return(density)
}
x=read.table("ten.csv")
x=x$V1
theta=maxlogL(x, dist="dpl",link=list(over=c("mu", "sigma"), fun= c("log_link", "log_link")))
summary(theta)
---------------------------------------------------------------
Optimization routine: nlminb
Standard Error calculation: Hessian from optim
---------------------------------------------------------------
AIC BIC
100.8022 96.8022
---------------------------------------------------------------
Estimate Std. Error
mu 3.877652 0.3175
sigma 0.048756 0.0159
-----
Then I used the same data and the same codes but i changed the distribution and i got error:
deeb = function(x, a , q , m, l , log = FALSE ) {
if (any(a <= 0))
stop(paste("a must be positive", "\n", ""))
if (any(q <= 0.5))
stop(paste("q must be greater than 0.5", "\n", ""))
if (any(m <= 0))
stop(paste("m must be positive", "\n", ""))
if (any(l <= 0))
stop(paste(" l must be positive", "\n", ""))
if (any(x <= -sqrt(2* m* l)) | any(x >= sqrt(2* m* l)))
stop(paste("x must be greater than -sqrt(2*l*m) and less than sqrt(2*l*m) ", "\n", ""))
loglik = log(a)+ log( gamma(q)) + log (gamma(a+0.5))+ log (gamma(q+0.5)) - (0.5*log(pi))-log( l)- log( m)- log ( gamma(q-0.5))- log ( gamma(1+q+a) )+(0.5*log((log(x))^2))-log(x) + log( hyperg_2F1( a+1, q+0.5, q+a+1,(1-((log(x))^2)/(l*m))))
if(log == FALSE)
density = exp(loglik)
else density = loglik
return (density)
}
theta1=maxlogL(x, dist = "deeb", link = list( over = c("a", "q", "m", "l"), fun = c("log_link", "log_link", "log_link", "log_link")))
Error in if (any(a <= 0)) stop(paste("a must be positive", "\n", "")) :
missing value where TRUE/FALSE needed
In addition: There were 28 warnings (use warnings() to see them)
> summary (theta1)
Error in h(simpleError(msg, call)) :
error in evaluating the argument 'object' in selecting a method for function 'summary': object 'theta1' not found