# using nls for gamma distribution (a,b,d) Classic List Threaded 3 messages Open this post in threaded view
|

## using nls for gamma distribution (a,b,d)

 Dear all i want to estimated the parameter of the gamma density(a,b,d) f(x) = (1/gamma(b)*(a^b)) * ((x-d)^(b-1)) * exp{-(x-d)/a)} for x>d f(x) = Age specific fertility rate x = age when i run this in R by usling nls() gamma.asfr <- formula(asfr ~ (((age-d)^(b-1))/((gamma(b))*(a^b)))* exp(-((age-d)/a))) gamma.asfr1 <- nls(gamma.asfr, data= asfr.aus, start = list(b = 28, a = 1, d= 0.5), trace = TRUE) error: Error in numericDeriv(form[[3L]], names(ind), env) :  Missing value or an infinity produced when evaluating the model when I use plinear algoritm, and run this gamma.asfr1 <- nls(gamma.asfr, data= asfr.aus, start = list(b = 28, a = 1, d= 0.5), trace = TRUE, algorithm="plinear") error: number of iterations exceeded maximum of 50 then i fix the iteration and minFactor even then its can't work gamma.asfr1 <- nls(gamma.asfr, data= asfr.aus, start = list(b = 28, a = 1, d= 0.5), trace = TRUE, algorithm="plinear", nls.control(maxiter=500, minFactor=0.000001)) error: Missing value or an infinity produced when evaluating the model Can any body tell the problem, where i am doing wrong thanks in advanc .. Muhammad Asif Wazir Ph.D student Institut für Statistik und Decision Support Systems (ISDS). University of Vienna, Austria cell: 00436509092298         [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
 I never had seen someone using a density function for a mean function in nonlinear regression. The convergence problem observed can be due the model derivatives respect to the parameters. How could computationally be d.gamma/d.alpha? digamma function? Does R understand this? If your phenomena exhibit a gamma density form you could use a mean function like: ```fun <- function(x, alpha, beta, theta){ x*(alpha+beta*x)^(-1/theta) } da <- data.frame(x=1:20) da\$y <- fun(da\$x, 1, 2, 0.9)+rnorm(da\$x,0,0.001) plot(y~x, da) curve(fun(x, 1,2,0.9), 1, 100, add=TRUE) n0 <- nls(y~fun(x, alpha, beta, theta), data=da, start=list(alpha=1, beta=2, theta=0.9)) ```Or you can use a gamma density without the integration equal one. So, you can replace the gamma(alpha) for B term: ```gammafun <- function(x, alpha, B, b){ (1/B)*alpha^b*x^(b-1)*exp(-x/alpha) } curve(gammafun(x, 5, 1, 2), 0, 20) da <- data.frame(x=1:20) da\$z <- gammafun(da\$x, 5, 1, 2)+rnorm(da\$x,0,1) plot(z~x, da) curve(gammafun(x, 5, 1, 2), add=TRUE) n1 <- nls(z~gammafun(x, alpha, B, b), data=da, start=list(alpha=5, B=1, b=2))```Sincerely. Walmes. ========================================================================== Walmes Marques Zeviani LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W) Departamento de Estatística - Universidade Federal do Paraná homepage: http://www.leg.ufpr.br/~walmes        twitter: @walmeszeviani ==========================================================================