Dear Dr Rick Thank you for your prompt response. I have to doubly integrate the bivariate normal density function of x and y with an unknown polyserial correlation, for eventually finding the maximum likelihood estimate of the polyserial correlation for complex survey. The cdf can’t be used as it’s not a full -ve infinity to +ve infinity limit for the integration according to the R manual for wcorr package for polyserial correlation. The limits are determined based on summing survey weights wrt to different cut points. I want to develop the program for sas although it’s available in R wcorr package. Looking forward to your help. Regards Haider