Rick, thanks for taking time to reply and look at it.  Your comments are greatly appreciated.  I took some time to consider your comments. First off, here is my code again, and I attach the data now in a data step: data d; input w2 ; cards; 0.086206897 11 1 2 4 1 0.5 0.153846154 0.246575342 0.0125 0.5 0.15625 0.033333333 0.4 0.192982456 11 1 2 13 5.166666667 1 0.333333333 1 9 0.153846154 0.037037037 1 0.333333333 0.181818182 1.5 0.133333333 2 1 2 3 37 0.076923077 0.666666667 14 1 1 1.666666667 0.142857143 1.666666667 20 0.272727273 1.5 0.222222222 0.7 0.25 0.5 1 4 0.454545455 1 1 0.578947368 3 3 10 1 5 2 1 1.5 0.2 2 0.242424242 0.25 1.066666667 8 4 2 0.166666667 0.454545455 0.046511628 0.5 0.2 0.2 0.2 2 3 34 0.25 ; run; proc iml; use d; read all into w2; n = nrow(w2); thres = 20; start logl(par) global(w2,thres,n); eksi = par[1]; theta1 = par[2]; theta2 = par[3]; dens = j(n,1,0); /*eksisq = eksi**2; nover2 = n/2; poissonparameter1 = theta1/2; poissonparameter2 = theta2/2; */ do index = 1 to n; val = 0; w2val = w2[index,1];   do i = 0 to thres;    do k1 = 0 to thres;    do k2 = 0 to thres;      negb = (1-eksi**2)**(n/2) * eksi**(2*i) * gamma(n/2 + i) / (gamma(n/2)*fact(i));      *negb = pdf("negbinomial",i,eksisq,nover2);       poiss = exp(-theta1/2) * (theta1/2)**k1 / fact(k1) * exp(-theta2/2) * (theta2/2)**k2 / fact(k2)  ;      *poiss = pdf("poisson",k2,poissonparameter1) * pdf("poisson",k2,poissonparameter2);      betav = 1/beta(n/2+i+k1+k2,n/2+i+k1+k2);      variab = w2val**(n/2+i+k1+k2-1) * (1 + w2val)**(-1*(n+2*i+k1+k2));         val = val + negb*betav*variab*poiss;    end; *k2;    end; *k1;   end; *i;   dens[index,1] = val; end; *index; v = dens[#,]; return(v); finish logl; *Calculate ML estimates; par={0.2 5 5}; z = logl(par); optn = {1 2}; call nlpnms(rc,xr,"logl",par,optn); z = logl(par); print z; xr = n // xr`; row = {"n" "eksi hat" "theta1 hat" "theta2 hat"}; if rc>0 then print xr[r=row label='ML Estimates for model:']; else print 'Algorithm does not converge'; run; quit; I also added another Poisson factor to my density. I spoke to a colleague regarding my use of NLPNRA and he suggested I also give NLPNMS a try since the latter does not rely so heavily on the second derivatives (i.e. Hessian matrix) being continuous.  Due to the nature of the model and the "brute force" attack, I tried it, and (sometimes...) it is successful.  Like you mentioned, the overflow is indeed reached quickly (which is mentioned in the log, and then the call does not execute further), and then I adjusted the "threshold" to a value such that the model will still give results, even though it is not very reliable.  I am mostly now aiming to just make sure my IML approach is correct and will consider the reliability of the estimates at a later stage. I can run z=logl(par), although it returns a value of 0 - why, I am not sure. As per your suggestion, I used built-in functions of the PDFs of SAS in this new attempt.  I added it here and is commented out for easy use when required.  A possible challenge in this regard is that the call seems to optimize the parameter "eksi" sometimes to be larger than 1, that results in the pdf function of the negative binomial to be undefined.  I kept the written out density in the code just in case. If you have any additional comments or suggestions I would greatly appreciate it. Johan ... View more