Thanks for your suggestion above. I modified the input to the quantile function by restricting it away from the endpoints 0 and 1. (Also now the integral of the binomial density itself only goes from 0.01 to 0.99 unlike (0,1) before.) I am wondering now if I can reformulate the entire model as a Random effects model. Specifically, since the likelihood function is being integrated over all possible values of "y", I am wondering if "y" can be considered as a random effect. In the existing model, the reason for integration in the likelihood function is that the binomial probability is a CONDITIONAL probability, conditional on values taken by y. (Also y really is unbounded- -infty to +infty, but through a change of variable using y = Normcdf(x), I restricted the integral to being over (0,1). Otherwise, the variable of integration would be x = normcdf, and dx = d(normcdf), which is difficult to write an integration formula for. ). Anyway I tried using Proc NLmixed, with the code below, but got this error: NOTE: A finite difference approximation is used for the derivative of the 'PDF' function. ERROR: Quadrature accuracy of 0.000100 could not be achieved with 31 points. The achieved accuracy was 0.503703. Here is the nlmixed code. proc nlmixed data=bin1 gconv=1e-10 maxiter=200000; parms rho=0.5; bounds rho >= 0.01, rho <= 0.99; /** gamma= 0.3 =const **/ /** Id is the observation number **/ w= (gamma - sqrt(rho)*w) / (1-sqrt(rho) ); p = cdf("Normal",w); /** Likelihood function **/ L= pdf("Binomial",D,p, N); model d ~ general(L); random y ~ Normal(0,1) subject=id;; run; Thanks so much for your all your guidance and help so far!
... View more