single and double integration of a multivariable function in SAS

[Rick_SAS]

Please post the program you have written. That will enable us to suggest modifications.

[Rick_SAS]

>  it produced a constant value for the integrand. however, i need the integrand as a function of rho

No, the program I posted is a function of rho. You can see how the function changes by plotting the double integral versus rho:

rhoVals = do(-0.9, 0.6, 0.05);
IntegVals = j(1, ncol(rhoVals), .);
do i = 1 to ncol(rhoVals);
   rho = rhoVals[i];
   call quad(v, "outer", xinterval);
   IntegVals[i] = v;
end;

title "Double Integral as a Function of rho";
call series(rhoVals, IntegVals) grid={x y};

As I have previously explained, this function is unbounded as rho -> +1. Therefore, it has no maximum on the interval (-1, 1).  

Do you have a reference for the MLE equations for weighted polychoric correlation?  Perhaps you copied an equation incorrectly?

[hrmannan]

Dear Dr Rick

the methodology for estimating weighted polychoric correlation taken from the doc

written by Paul Bailey, Ahmad Emad, Ting Zhang, Qingshu Xie,Source: vignettes/wCorrFormulas.Rmd.

the link for the doc is-

https://american-institutes-for-research.github.io/wCorr/articles/wCorrFormulas.html

 i couldnt paste the formulas as the format is not readable.  kindly see where i went wrong with the formulas . First the double integral of a bivariate normal density over finite limits (both are proportions so below 1) has to be solved and then the loglikelihood is to be maximized, the log likelihood is-

=sum(weights*log of the double integral)

you showed that this double integral doesnt converge. then where is the error.

looking forward to your kind help.

regards

Haider Mannan

[hrmannan]

dear Dr. Rick, I have sent you link for the paper. after taking log of the function and double-integrating it still seems to give unbounded function. the graph below shows this. how do I get a bounded function so I can maximise sum(weights*integrand) with respect to rho using my data when weights are survey weights? kindly help. The SAS codes & outputs are:

/* double integration for loglikelihood for estimating weighted polychoric correlation */
proc iml;
start inner(y) global(g_x, rho);
x = g_x;
return(-.5*x#x/(1-rho**2)-.5*y#y+rho#x#y/(1-rho**2)) ;
finish;

/* TEST THE INNER MODULE */
rho=0.5; /* global value of RHO */
g_x = 0; /* for testing: evaluate inner integral when x=g_x */
yinterval = {.5 .8};
call quad(z, "inner", yinterval);
print z;
free g_x; /* undefine g_x; it'll be set in the OUTER module */

/* g_x is the global variable that holds the value of x for the outer integral.
This is a constant when we evaluate the inner integral. */
start outer(x) global(g_x, rho);
g_x = x; /* save copy of x so the inner integral can use it */
yinterval = 1.5 || 3;
/* evaluate w = inner integral when x is held constant */
call quad(w, "inner", yinterval);
return (w);
finish;

/* TEST THE OUTER MODULE */
xinterval= {.6 .9};
call quad(v, "outer", xinterval);
print v;

rhoVals = do(-0.9, .6, 0.05);
IntegVals = j(1, ncol(rhoVals), .);
do i = 1 to ncol(rhoVals);
rho = rhoVals[i];
call quad(v, "outer", xinterval);
IntegVals[i] = v;
end;

title "Double Integral as a Function of rho";
call series(rhoVals, IntegVals) grid={x y};

Double Integral as a Function of rho

z

-0.0645

v

-0.846

---

[hrmannan]

Dear Dr Rick

the methodology for estimating weighted polychoric correlation taken from the doc

written by Paul Bailey, Ahmad Emad, Ting Zhang, Qingshu Xie,Source: vignettes/wCorrFormulas.Rmd.

the link for the doc is-

https://american-institutes-for-research.github.io/wCorr/articles/wCorrFormulas.html

 i couldnt paste the formulas as the format is not readable.  kindly see where i went wrong with the formulas . First the double integral of a bivariate normal density over finite limits (both are proportions so below 1) has to be solved and then the loglikelihood is to be maximized, the log likelihood is-

=sum(weights*log of the double integral)

you showed that this double integral doesnt converge. then where is the error.

looking forward to your kind help.

regards

Haider Mannan

[Rick_SAS]

Thanks for the link. You should remove the "Problem is solved" indicator if you want people to work on this problem. If you mark it as solved, people assume your issue is resolved.

[hrmannan]

dear Dr. Rick

maybe I didnt send the full link for the R doc on calculating weighted polychoric correlation. here it is-

file://ad.uws.edu.au/dfshare/HomesCMB$/30043787/My%20Documents/Phillipa's%20ABS%20data%20based%20study%20aug%202022/wCorr%20Formulas%20%E2%80%A2%20wCorr.html

[Rick_SAS]

Let me make sure I understand what you want. In the article

Bailey, Emad, Zhang, and Xie (2023)
"wCorr Formulas"
https://american-institutes-for-research.github.io/wCorr/articles/wCorrFormulas.html

you want to calculate the probability given by the following double integral inside the red box?

[hrmannan]

Dear Dr Rick
Using the article formulas I want to determine the weighted polychoric correlation (rho). The first step is to find the joint probability function which is a double integral of a bivariate density of two standard normal (latent) variables over finite boundaries as given which have to be determined comparing the observed ordinal variables against sum of survey weights in my data. The next step is to use this result of double integral to my data to maximise the log likelihood function with respect to rho when loglikelihood is sum(weight*log(double integral)). In the lit search I found Simpsons one third rule being suggested to find the double integral of a bivariate normal density. And there are other ways.thanks.
Haider

[hrmannan]

Yes the formulas you mention from the article are the ones to be used. I have just explained the steps as I understood. Thanks

[Rick_SAS]

Great. The integrand in that formula is the standard bivariate normal density function (PDF). You do not need to use the QUAD function to compute the double integral. Instead, you can use the PROBBNRM function in SAS to integrate the bivariate PDF. The PROBBNRM function assumes that the lower limits of integration are minus infinity (just like the CDF function does), but you can use the rules of integration to obtain the integral over a rectangular region.

Let (a,b) be the limits of integration for the X variables. (These are the "theta prime" limits in the paper.) Let (c,d) be the limits of integration for the Y variable. (These are the "theta" limits in the paper.) Then you can use the following function to evaluate the double integral on the rectangle [a,b]x[c,d]:

proc iml;
/* compute the probability under the bivariate normal cdf with 
   correlation parameter rho for the rectangular region [a,b]x[c,d] */
start BivarProb(a, b, c, d, rho);
   Prob = PROBBNRM(b,d,rho) - PROBBNRM(a,d,rho) - PROBBNRM(b,c,rho) + PROBBNRM(a,c,rho);
   return Prob;
finish;

/* define region of integration to be the rectangular region ([a,b]x[c,d] */
a = -1;   b = 0;
c =  0;   d = 0.5;

/* prob when rho = 0 */
rho = 0;
Prob = BivarProb(a, b, c, d, rho);
print rho Prob;
/* prob when rho = -0.5 */
rho = -0.5;
Prob = BivarProb(a, b, c, d, rho);
print rho Prob;

/* compute the prob for a range of rho values */
rhos = do(-0.9, 0.9, 0.05);
Probs = j(1, ncol(rhos), .);
do i = 1 to ncol(rhos);
   rho = rhos[i];
   Probs[i] = BivarProb(a, b, c, d, rho);
end;
title "Pr[ (x,y) in [a,b]x[c,d] ]";
title2 "(a,b,c,d) = (-1,0,0,0.5)";
call series(rhos, Probs) grid={x y};

For the full program, call the BivarProb function from inside the objective function that you are optimizing. 

If this is the information you want, then please mark this response as the correct solution so that future programmers will be able to find the solution, too.

[hrmannan]

Dear Dr Rick, thanks very much. I am trying to understand the outputs. Does the graph mean that the ML estimate of rho is aroundr -0.9 where the probability is maximum? 

another point is that the regions of integration you have assigned are (-1,0) and (0,0.5). According to the paper i sent, shouldn't these inner and outer limits of integration be non-negative fractions? please let me know.

Haider

[Rick_SAS]

My response shows how you can compute the bivariate normal probability by using the PROBBNRM function in SAS. You can use this function instead of the QUAD function to compute the double integrals in your problem. My response does not do any optimization, nor does it compute the polychoric correlation. You will have to write that part of the program yourself. 

> Does the graph mean that the ML estimate of rho is aroundr -0.9 where the probability is maximum? 

No. This is not the graph of the LL. It is a graph of the bivariate normal probability (the double integral) in the rectangle [-1, 0] x [0, 0.5] for various choices of rho.

> According to the paper i sent, shouldn't these inner and outer limits of integration be non-negative fractions?

No, I don't think so. The limits of integration are the cut points of the underlying (latent) normal variate which gives rise to the (observed) ordinal variable. Yes, they are based on the observed cumulative probabilities of the marginal counts, which are nonnegative fractions, but the actual cut points are of the form quantile("Normal", cumprob[i]), where cumprob[i] are the cumulative probabilities. So, the cut points are in (-infinity, infinity).

[hrmannan]

Dear Dr Rick
Thanks. Please explain how you got the limits (-0.5,0) and (0,1). The paper stated two ways or is it just one method. Kindly explain.
Haider

[hrmannan]

Correction I meant how did you find the limits (-1,0) and (0,0.5) for the probability integral of a bivariate normal density?
Haider