Forum: help
Monitor Forum | 
Start New Thread
| RE: Simulation using maxLik().......... [ Reply ] By: Arne Henningsen on 2015-04-07 22:26 | [forum:42140] |
|
Dear Rajibul You could call maxLik() within the try() command to avoid that the simulation terminates at an error message. If you provide a *minimal* reproducible example, it is more likely that somebody reads your code. Best regards, Arne |
|
| Simulation using maxLik().......... [ Reply ] By: Rajibul Mian on 2015-04-07 15:21 | [forum:42139] |
|
Dear All: Hope everybody is fine. I am trying to run the following (given below) small simulation and find some problems. In my simulation maxLik() needs to be run several times. I find out that some times maxLik() can converge sometimes end up showing the error message. Lets say I want to run simulation 10 times. then sometimes I got the results with 10 possible solutions that indicates maxLik() runs well in all 10 occasions. But in many times maxLik() stooped and showing the error message which also stopped the simulation. My concern is to ignore the error messages and go to the next step ....to continue the simulation. For example: say in 6th occasion among the 10 time of simulation maxLik() couldn't converge and show the error message, then I want to ignore that error message and continue to the 7th occasion. With my very little knowledge, I realized that maxLik() is programmed in such a way that it shows the error message in its internal iteration, which is very hard to manage (for me) by using the available error/message suppressing tools.... I am giving my code (tried to give the minimal one but....) here for the consideration. Any kind of suggestions will be accepted with great appreciation. Thanks a lot. Best regards. Rajibul Mian rajibulmian@gmail.com rajibulmian.r@gmail.com PhD candidate, Math Stat, Univ of Windsor ###################################### ############## Code start ############## ###################################### library(gamlss.dist) library(gamlss) library(maxLik) ### ### function for check ### rm(list = ls()) rajib_1 = function(n, k, par){ ### function starts mat = matrix(,ncol = 6, nrow = k, byrow = T) for(i in 1:k){ ### for loop starts # par = c(-2,1,-1,2,2,0.3) ### for checking # n = 100 x_a1 = 0 x_a2 = 0.5 x = rnorm(n,x_a1,x_a2) #x = runif(n, -10,5) x g_a1 = 3 g_a2 = 0.5 g = rnorm(n,g_a1,g_a2) #g = runif(n, -50,50) g z_a1 = -1 z_a2 = 1 z_a3 = 0.5 z_p = exp(z_a1+z_a2*x+z_a3*g) / (1+exp(z_a1+z_a2*x+z_a3*g)) z_size = 1 z = rbinom(n,z_size,round(mean((z_p)),1)) #z = rbinom(n,z_size,mean((z_p))) #z = rbinom(n,z_size,0.5) z a = par theta = a[1] + a[2]*x + a[3]*g + a[4]*z c = a[5] m = (1/c) * ( exp(theta)/(1-exp(theta)) ) ### m needs to be increased around 1 or more than 1 data_1 = rZINBI(n, mu = 1+round(mean(na.omit(m))), sigma = a[5], nu = a[6]) #data_1 = rZINBI(n, mu = abs(m), sigma = a[5], nu = a[6]) y = data_1[data_1!=0] #data = rZINBI(100, mu = 1, sigma = 2, nu = 0.3) #y = data[data!=0] ### func_2 = function(par){ ### function starts a = par theta = a[1] + a[2]*x + a[3]*g + a[4]*z c = a[5] delta = a[6] loglik_2 = sum(( - log(1+delta) + log( delta + exp( (1/c) * log( 1-exp(theta) ) ) ) + y*theta + (1/c)*log( 1-exp(theta) ) + log( gamma(y+1/c) / (gamma(1/c)*factorial(y)) ) ) ) loglik_2 } ### function ends model_1.2 = maxLik(logLik = func_2, start = a, iterlim = 500) model_1.2 #model_1.3 = suppressMessages(model_1.2) #if(is(model_1.3)=="try-error")next model_1.2_est = model_1.2$est # mat[i,] = c(mean(na.omit(m)),model_1.2_est) #mat[i,] = model_1.2_est mat[i,] = summary(model_1.2)$est[,1] } ### for loop ends return(mat) } ### function ends # par = c(-2,1,-1,2,0.3,0.3) #n = 100 res = rajib_1(n = 100, k = 10, par = c(-2,1,-1,2,2,0.3) ) res Thanks a lot ###################################### ############## Code end ############## ###################################### #### Error message look like..... Iteration 23 Parameter: [1] -4.393563954 0.026897608 0.105529277 -0.169426621 0.005932189 [6] 475.534579010 Gradient: [,1] [,2] [,3] [,4] [,5] [,6] [1,] NaN NaN NaN NaN NaN NaN Error in maxNRCompute(fn = logLikAttr, fnOrig = fn, gradOrig = grad, hessOrig = hess, : NA in gradient In addition: There were 50 or more warnings (use warnings() to see the first 50) ########################### |
|
