I sent Ricardo <[email protected]> some Mata code that
computes the orthogonal polynomial codings for an integer
variable and I thought that some people learning Mata might
find it a useful.
-Rich
[email protected]
. clear
. set obs 9
obs was 0, now 9
. gen int cat = mod(_n-1,3)+1
. mata:orthogonal_poly("cat")
2
3
. list
+-------------------------------+
| cat _cat_1 _cat_2 |
|-------------------------------|
1. | 1 -.70710678 .40824829 |
2. | 2 0 -.81649658 |
3. | 3 .70710678 .40824829 |
4. | 1 -.70710678 .40824829 |
5. | 2 0 -.81649658 |
|-------------------------------|
6. | 3 .70710678 .40824829 |
7. | 1 -.70710678 .40824829 |
8. | 2 0 -.81649658 |
9. | 3 .70710678 .40824829 |
+-------------------------------+
. clear
. set obs 12
obs was 0, now 12
. gen int cat = mod(_n-1,4)+1
. mata:orthogonal_poly("cat")
2
3
4
. list
+----------------------------------------+
| cat _cat_1 _cat_2 _cat_3 |
|----------------------------------------|
1. | 1 -.67082039 .5 -.2236068 |
2. | 2 -.2236068 -.5 .67082039 |
3. | 3 .2236068 -.5 -.67082039 |
4. | 4 .67082039 .5 .2236068 |
5. | 1 -.67082039 .5 -.2236068 |
|----------------------------------------|
6. | 2 -.2236068 -.5 .67082039 |
7. | 3 .2236068 -.5 -.67082039 |
8. | 4 .67082039 .5 .2236068 |
9. | 1 -.67082039 .5 -.2236068 |
10. | 2 -.2236068 -.5 .67082039 |
|----------------------------------------|
11. | 3 .2236068 -.5 -.67082039 |
12. | 4 .67082039 .5 .2236068 |
+----------------------------------------+
mata
void function orthogonal_poly(string scalar svar)
{
real scalar m, ni
string scalar names, namej
real colvector v, lev
real matrix X
pointer (real colvector) colvector ind
pragma unset v
st_view(v,.,svar)
lev = sort(uniqrows(v),1)
m = length(lev)
if (m<2) {
errprintf("too few levels %d\n", m)
exit(198)
}
if (m>100) {
errprintf("too many levels %d\n", m)
exit(198)
}
ind = J(m,1,NULL)
X = orthogonal_poly_X(m)[.,2::m]
for (i=1; i<=m; i++) {
ind[i] = &orthogonal_poly_index(v,lev[i])
}
names = ""
for (j=1; j<m; j++) {
namej = "_"+svar+"_"+strofreal(j)
stata("cap drop "+namej)
st_addvar("double",namej)
names = names + " " + namej
for (i=1; i<=m; i++) {
ni = length(*ind[i])
if (ni > 0) st_store(*ind[i],namej,J(ni,1,X[i,j]))
}
}
st_global("r(polyvars)",names)
}
/* generate orthogonal polynomials */
real matrix orthogonal_poly_X(real scalar n)
{
real scalar i
real colvector z
real matrix X, Q, R
z = (1::n):-mean((1::n),J(n,1,1))
X = J(n,n,1)
for (i=2; i<=n; i++) X[.,i] = X[.,i-1]:*z
pragma unset Q
pragma unset R
qrd(X,Q,R)
X = Q*diag(R)
X = X:/sqrt(colsum(X:^2))
X = X:*(abs(X):>1.0e-15)
return(X)
}
real colvector orthogonal_poly_index(real colvector v, real scalar lev)
{
real scalar n, k, i, j
real colvector ind, iv
n = length(v)
iv = (v:==lev)
k = sum(iv)
if (k == 0) return(J(0,1,.))
ind = J(k,1,.)
j = 0
for (i=1; i<=n; i++) {
if (iv[i]>0) ind[++j] = i
}
return(ind)
}
end
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
