mle estimate for model containing hypergeometric function

deeb = function(x, a , q , m, l , log = FALSE ){

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", ""))

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", ""))

loglik = log( a* gamma(q) *gamma(a+0.5) *gamma(q+0.5)) - log( (pi^0.5)* l * m* gamma(q-0.5) * gamma(1+q+a) )+ (0.5*log((x^2))) + log( hyperg_2F1( a+1, q+0.5, q+a+1,(1-((x^2)/(l*m)))))

if(log == FALSE)

density = exp(loglik)

else density = loglik

return (density)

}

theta=maxlogL(x, dist = "deeb", link = list( over = c("a", "q", "m", "l"), fun = c("log_link", "log_link", "log_link", "log_link")))

summary (theta)

I wrote this code to find the estimate of the parameter of my new distribution but it give me the below error

Error in if (any(x < -sqrt(2 * m * l)) | any(x > sqrt(2 * m * l))) stop(paste("x must be greater than -sqrt(2lm) and less than sqrt(2lm) ", :

missing value where TRUE/FALSE needed

In addition: There were 31 warnings (use warnings() to see them)

summary (theta)

Error in h(simpleError(msg, call)) :

error in evaluating the argument 'object' in selecting a method for function 'summary': object 'theta' not found

could u please help me to correct the error. i have here constraint on x because of the hypergeometric function

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

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.