*! *! mm_kern.mata *! version 1.0.1, Ben Jann, 30aug2006 version 9.2
local kernels /// epanechnikov /// epan2 /// biweight /// triweight /// cosine /// gaussian /// parzen /// rectangle /// triangle
local klist tokens("`kernels'")
local kname `"(k!="" ? k : "epan2")"'
mata:
string scalar _mm_unabkern(|string scalar k) { return(mm_strexpand(k, `klist', "epan2", 0, k+": invalid kernel")) }
pointer scalar _mm_findkern(|string scalar k) { pointer(real matrix function) scalar K K = findexternal("mm_kern_"+`kname'+"()") if (K==NULL) _error(3499,"mm_kern_"+`kname'+"() not found") return(K) }
pointer scalar _mm_findkint(|string scalar k) { pointer(real matrix function) scalar K K = findexternal("mm_kint_"+`kname'+"()") if (K==NULL) _error(3499,"mm_kint_"+`kname'+"() not found") return(K) }
pointer scalar _mm_findkdel0(|string scalar k) { pointer(real scalar function) scalar K K = findexternal("mm_kdel0_"+`kname'+"()") if (K==NULL) _error(3499,"mm_kdel0_"+`kname'+"() not found") return(K) }
real matrix mm_kern(string scalar k, real matrix x) { return((*_mm_findkern(_mm_unabkern(k)))(x)) }
real matrix mm_kint(string scalar k, real scalar r, |real matrix x) { return(mm_callf(mm_callf_setup(_mm_findkint(_mm_unabkern(k)), args()-2, x), r) > ) }
real matrix mm_kdel0(string scalar k) { return((*_mm_findkdel0(_mm_unabkern(k)))()) }
// mm_kern_*: kernel function // mm_kint_*(1,x): integral of K(z) from -infty to x // mm_kint_*(2,x): integral of K(z)^2 from -infty to x // mm_kint_*(2): kernel "roughness" // mm_kint_*(3,x): integral of z*K(z) from -infty to x // mm_kint_*(4,x): integral of z^2*K(z) from -infty to x // mm_kint_*(4): kernel variance // mm_kvar_*: variance // mm_kdel0_*: canonical bandwidth
// (scaled) Epanechnikov kernel function real matrix mm_kern_epanechnikov(real matrix x) { return((.75:-.15*x:^2)/sqrt(5):*(abs(x):<sqrt(5))) } real matrix mm_kint_epanechnikov(real scalar r, |real matrix x) { if (r==1 & args()<2) return(1) if (r==1) return((.5:+(.75*x-.05*x:^3)/sqrt(5)):*(abs(x):<sqrt(5)) + (x:>=sqrt(5))) if (r==2 & args()<2) return(3/5^1.5) if (r==2) return((.3/sqrt(5):+(.1125*x-.015*x:^3+.0009*x:^5)):*(abs(x):<sqrt(5)) + 3/5^ > 1.5*(x:>=sqrt(5))) if (r==3 & args()<2) return(0) if (r==3) return((-.1875*sqrt(5):+.375/sqrt(5)*x:^2-.0375/sqrt(5)*x:^4):*(abs(x):<sqrt( > 5))) if (r==4 & args()<2) return(1) if (r==4) return((.5:+.25/sqrt(5)*x:^3-.03/sqrt(5)*x:^5):*(abs(x):<sqrt(5)) + (x:>=sqrt > (5))) else _error(3300) } //real scalar mm_kvar_epanechnikov() return(1) real scalar mm_kdel0_epanechnikov() return((3/5^1.5)^.2)
// Epanechnikov kernel function real matrix mm_kern_epan2(real matrix x) { return((.75:-.75*x:^2):*(abs(x):<1)) } real matrix mm_kint_epan2(real scalar r, |real matrix x) { if (r==1 & args()<2) return(1) if (r==1) return((.5:+.75*x-.25*x:^3):*(abs(x):<1) + (x:>=1)) if (r==2 & args()<2) return(.6) if (r==2) return((.3:+.5625*x-.375*x:^3+.1125*x:^5):*(abs(x):<1) + .6*(x:>=1)) if (r==3 & args()<2) return(0) if (r==3) return((-.1875:+.375*x:^2-.1875*x:^4):*(abs(x):<1)) if (r==4 & args()<2) return(.2) // 1/5 if (r==4) return((.1:+.25*x:^3-.15*x:^5):*(abs(x):<1) + .2*(x:>=1)) else _error(3300) } //real scalar mm_kvar_epan2() return(.2) real scalar mm_kdel0_epan2() return(15^.2)
// Biweight kernel function real matrix mm_kern_biweight(real matrix x) { return(.9375*(1:-x:^2):^2:*(abs(x):<1)) } real matrix mm_kint_biweight(| real scalar r, real matrix x) { if (r==1 & args()<2) return(1) if (r==1) return((.5:+.9375*(x-2/3*x:^3+.2*x:^5)):*(abs(x):<1) + (x:>=1)) if (r==2 & args()<2) return(5/7) if (r==2) return((5/14:+225/256*x-75/64*x:^3+135/128*x:^5-225/448*x:^7+25/256*x:^9):*(a > bs(x):<1) + 5/7*(x:>=1)) if (r==3 & args()<2) return(0) if (r==3) return((-.15625:+.46875*x:^2-.46875*x:^4 +.15625*x:^6):*(abs(x):<1)) if (r==4 & args()<2) return(1/7) if (r==4) return((1/14:+.3125*x:^3-.375*x:^5+15/112*x:^7):*(abs(x):<1) + (x:>=1)/7) else _error(3300) } //real scalar mm_kvar_biweight() return(1/7) real scalar mm_kdel0_biweight() return(35^.2)
// Triweight kernel function real matrix mm_kern_triweight(real matrix x) { return(1.09375*(1:-x:^2):^3:*(abs(x):<1)) } real matrix mm_kint_triweight(| real scalar r, real matrix x) { if (r==1 & args()<2) return(1) if (r==1) return((.5:+1.09375*(x-x:^3+.6*x:^5-1/7*x:^7)):*(abs(x):<1) + (x:>=1)) if (r==2 & args()<2) return(350/429) if (r==2) return((175/429:+1225/1024*(x-2*x:^3+3*x:^5-20/7*x:^7+5/3*x:^9-6/11*x:^11+1/1 > 3*x:^13)):*(abs(x):<1) + 350/429*(x:>=1)) if (r==3 & args()<2) return(0) if (r==3) return((-35/256:+35/64*x:^2-105/128*x:^4+35/64*x:^6-35/256*x:^8):*(abs(x):<1) > ) if (r==4 & args()<2) return(1/9) if (r==4) return((1/18:+35/96*x:^3-21/32*x:^5+15/32*x:^7-35/288*x:^9):*(abs(x):<1) + (x > :>=1)/9) else _error(3300) } //real scalar mm_kvar_triweight() return(1/9) real scalar mm_kdel0_triweight() return((9450/143)^.2)
// Cosine trace real matrix mm_kern_cosine(real matrix x) { return((1:+cos(2*pi()*x)):*(abs(x):<.5)) } real matrix mm_kint_cosine(| real scalar r, real matrix x) { if (r==1 & args()<2) return(1) if (r==1) return((.5:+x+sin(2*pi()*x)/(2*pi())):*(abs(x):<.5) + (x:>=.5)) if (r==2 & args()<2) return(1.5) if (r==2) return((.75:+1.5*x+sin(2*pi()*x)/pi()+cos(2*pi()*x):*sin(2*pi()*x)/(4*pi())): > *(abs(x):<.5) + 1.5*(x:>=.5)) if (r==3 & args()<2) return(0) if (r==3) return((-.125+.25/pi()^2:+cos(2*pi()*x)/(4*pi()^2)+x:*sin(2*pi()*x)/(2*pi())+ > .5*x:^2):*(abs(x):<.5)) if (r==4 & args()<2) return(.5*(1/6-1/pi()^2)) if (r==4) return((1/24-.25/pi()^2:-sin(2*pi()*x)/(4*pi()^3)+x:*cos(2*pi()*x)/(2*pi()^2) > +x:^2:*sin(2*pi()*x)/(2*pi())+x:^3/3):*(abs(x):<.5) + .5*(1/6-1/pi()^2)*(x:>= > .5)) else _error(3300) } //real scalar mm_kvar_cosine() return(.5*(1/6-1/pi()^2)) real scalar mm_kdel0_cosine() return((3/(.5*(1/6-1/pi()^2)^2))^.2)
// Gaussian kernel function real matrix mm_kern_gaussian(real matrix x) { return(normalden(x)) // alternatively: exp(-0.5*((x)^2))/sqrt(2*pi()) } real matrix mm_kint_gaussian(| real scalar r, real matrix x) { if (r==1 & args()<2) return(1) if (r==1) return(normal(x)) if (r==2 & args()<2) return(.5/sqrt(pi())) // = normalden(0,sqrt(2)) if (r==2) return(.5/sqrt(pi()) * normal(sqrt(2)*x)) if (r==3 & args()<2) return(0) if (r==3) return(-1/sqrt(2*pi())*exp(-.5*x:^2)) if (r==4 & args()<2) return(1) if (r==4) return(-1/sqrt(2*pi())*x:*exp(-.5*x:^2) + normal(x)) else _error(3300) } //real scalar mm_kvar_gaussian() return(1) real scalar mm_kdel0_gaussian() return((1/(4*pi()))^.1)
// Parzen kernel function real matrix mm_kern_parzen(real matrix x) { return(((4/3:-8*x:^2+8*abs(x):^3):*(abs(x):<=.5)) + ((8*(1:-abs(x)):^3/3):*(abs(x):>1/2:&abs(x):<=1))) } real matrix mm_kint_parzen(| real scalar r, real matrix x) { if (r==1 & args()<2) return(1) if (r==1) return((2/3:+8/3*(x+1.5*x:^2+x:^3+.25*x:^4)):*(x:>=-1:&x:<-.5) + (.5:+4/3*x-8/3*x:^3-2*x:^4):*(x:>=-.5:&x:<0) + (.5:+4/3*x-8/3*x:^3+2*x:^4):*(x:>=0:&x:<=.5) + (1/3:+8/3*(x-1.5*x:^2+x:^3-.25*x:^4)):*(x:>.5:&x:<=1) + (x:>1)) if (r==2 & args()<2) return(302/315) if (r==2) return((64/63:+64/9*(x+3*x:^2+5*x:^3+5*x:^4+3*x:^5+x:^6+x:^7/7)):*(x:>=-1:&x: > <-.5) + (151/315:+16/9*(x-4*x:^3-3*x:^4+36/5*x:^5+12*x:^6+36/7*x:^7)):*(x:>=-.5:&x:< > 0) + (151/315:+16/9*(x-4*x:^3+3*x:^4+36/5*x:^5-12*x:^6+36/7*x:^7)):*(x:>=0:&x:<=. > 5) + (-2/35:+64/9*(x-3*x:^2+5*x:^3-5*x:^4+3*x:^5-x:^6+x:^7/7)):*(x:>.5:&x:<=1) + > 302/315*(x:>1)) if (r==3 & args()<2) return(0) if (r==3) return((-2/15:+4/3*x:^2+8/3*x:^3+2*x:^4+8/15*x:^5):*(x:>=-1:&x:<-.5) + (-7/60:+2/3*x:^2-2*x:^4-1.6*x:^5):*(x:>=-.5:&x:<0) + (-7/60:+2/3*x:^2-2*x:^4+1.6*x:^5):*(x:>=0:&x:<=.5) + (-2/15:+4/3*x:^2-8/3*x:^3+2*x:^4-8/15*x:^5):*(x:>.5:&x:<=1)) if (r==4 & args()<2) return(1/12) if (r==4) return((2/45:+8/9*x:^3+2*x:^4+1.6*x:^5+4/9*x:^6):*(x:>=-1:&x:<-.5) + (1/24:+4/9*x:^3-1.6*x:^5-4/3*x:^6):*(x:>=-.5:&x:<0) + (1/24:+4/9*x:^3-1.6*x:^5+4/3*x:^6):*(x:>=0:&x:<=.5) + (7/180:+8/9*x:^3-2*x:^4+1.6*x:^5-4/9*x:^6):*(x:>.5:&x:<=1) + (x:>1)/12) else _error(3300) } //real scalar mm_kvar_parzen() return(1/12) real scalar mm_kdel0_parzen() return((4832/35)^.2)
// Rectangle kernel function real matrix mm_kern_rectangle(real matrix x) { return(.5*(round(abs(x),1e-8):<1)) } // 1e-8 same as in kdensity.ado, version 2.3.6 26jun2000 (Stata 7) real matrix mm_kint_rectangle(| real scalar r, real matrix x) { if (r==1 & args()<2) return(1) if (r==1) return((.5:+.5*x):*(round(abs(x),1e-8):<1) + (round(x,1e-8):>=1)) if (r==2 & args()<2) return(.5) if (r==2) return((.25:+.25*x):*(round(abs(x),1e-8):<1) + .5*(round(x,1e-8):>=1)) if (r==3 & args()<2) return(0) if (r==3) return((-.25:+.25*x:^2):*(round(abs(x),1e-8):<1)) if (r==4 & args()<2) return(1/3) if (r==4) return((1/6:+x:^3/6):*(round(abs(x),1e-8):<1) + (round(x,1e-8):>=1)/3) else _error(3300) } //real scalar mm_kvar_rectangle() return(1/3) real scalar mm_kdel0_rectangle() return((9/2)^.2)
// Triangle kernel function real matrix mm_kern_triangle(real matrix x) { return((1:-abs(x)):*(abs(x):<1)) } real matrix mm_kint_triangle(| real scalar r, real matrix x) { if (r==1 & args()<2) return(1) if (r==1) return((.5:+x+.5*x:^2):*(x:>-1:&x:<0) + (.5:+x-.5*x:^2):*(x:>=0:&x:<1) + (x:>=1)) if (r==2 & args()<2) return(2/3) if (r==2) return((x+x:^2+(1:+x:^3)/3):*(x:>-1:&x:<0) + (x-x:^2+(1:+x:^3)/3):*(x:>=0:&x:<1) + 2/3*(x:>=1)) if (r==3 & args()<2) return(0) if (r==3) return((-1/6:+.5*x:^2+x:^3/3):*(x:>-1:&x:<0) + (-1/6:+.5*x:^2-x:^3/3):*(x:>=0:&x:<1)) if (r==4 & args()<2) return(1/6) if (r==4) return((1/12:+x:^3/3+.25*x:^4):*(x:>-1:&x:<0) + (1/12:+x:^3/3-.25*x:^4):*(x:>=0:&x:<1) + (x:>=1)/6) else _error(3300) } //real scalar mm_kvar_triangle() return(1/6) real scalar mm_kdel0_triangle() return(24^.2)
end
*! *! mm_quantile.mata *! version 1.0.8 20dec2007 Ben Jann version 9.2 mata:
real matrix mm_quantile(real matrix X, | real colvector w, real matrix P, real scalar altdef) { real rowvector result real scalar c, cX, cP, r, i
if (args()<2) w = 1 if (args()<3) P = (0, .25, .50, .75, 1)' if (args()<4) altdef = 0 if (cols(X)==1 & cols(P)!=1 & rows(P)==1) return(mm_quantile(X, w, P', altdef)') if (missing(P) | missing(X) | missing(w)) _error(3351) if (rows(w)!=1 & rows(w)!=rows(X)) _error(3200) r = rows(P) c = max(((cX=cols(X)), (cP=cols(P)))) if (cX!=1 & cX<c) _error(3200) if (cP!=1 & cP<c) _error(3200) if (rows(X)==0 | r==0 | c==0) return(J(r,c,.)) if (c==1) return(_mm_quantile(X, w, P, altdef)) result = J(r, c, .) if (cP==1) for (i=1; i<=c; i++) result[,i] = _mm_quantile(X[,i], w, P, altdef) else if (cX==1) for (i=1; i<=c; i++) result[,i] = _mm_quantile(X, w, P[,i], altdef) else for (i=1; i<=c; i++) result[,i] = _mm_quantile(X[,i], w, P[,i], altdef) return(result) }
real colvector _mm_quantile( real colvector X, real colvector w, real colvector P, real scalar altdef) { real colvector g, j, j1, p real scalar N
if (w!=1) return(_mm_quantilew(X, w, P, altdef)) N = rows(X) p = order(X,1) if (altdef) g = P*N + P else g = P*N j = floor(g) if (altdef) g = g - j else g = 0.5 :+ 0.5*((g - j):>0) j1 = j:+1 j = j :* (j:>=1) _editvalue(j, 0, 1) j = j :* (j:<=N) _editvalue(j, 0, N) j1 = j1 :* (j1:>=1) _editvalue(j1, 0, 1) j1 = j1 :* (j1:<=N) _editvalue(j1, 0, N) return((1:-g):*X[p[j]] + g:*X[p[j1]]) }
real colvector _mm_quantilew( real colvector X, real colvector w, real colvector P, real scalar altdef) { real colvector Q, pi, pj real scalar i, I, j, jj, J, rsum, W pointer scalar ww
I = rows(X) ww = (rows(w)==1 ? &J(I,1,w) : &w) if (altdef) return(_mm_quantilewalt(X, *ww, P)) W = quadsum(*ww) pi = order(X, 1) if (anyof(*ww, 0)) { pi = select(pi,(*ww)[pi]:!=0) I = rows(pi) } pj = order(P, 1) J = rows(P) Q = J(J, 1, .) j = 1 jj = pj[1] rsum = 0 for (i=1; i<=I; i++) { rsum = rsum + (*ww)[pi[i]] if (i<I) { if (rsum<P[jj]*W) continue if (X[pi[i]]==X[pi[i+1]]) continue } while (1) { if (rsum>P[jj]*W | i==I) Q[jj] = X[pi[i]] else Q[jj] = (X[pi[i]] + X[pi[i+1]])/2 j++ if (j>J) break jj = pj[j] if (i<I & rsum<P[jj]*W) break } if (j>J) break } return(Q) }
real colvector _mm_quantilewalt( real colvector X, real colvector w, real colvector P) { real colvector Q, pi, pj real scalar i, I, j, jj, J, rsum, rsum0, W, ub, g
W = quadsum(w) + 1 pi = order(X, 1) if (anyof(w, 0)) pi = select(pi, w[pi]:!=0) I = rows(pi) pj = order(P, 1) J = rows(P) Q = J(J, 1, .) rsum = w[pi[1]] for (j=1; j<=J; j++) { jj = pj[j] if (P[jj]*W <= rsum) Q[jj] = X[pi[1]] else break } for (i=2; i<=I; i++) { rsum0 = rsum rsum = rsum + w[pi[i]] if (i<I & rsum < P[jj]*W) continue while (1) { ub = rsum0+1 if (P[jj]*W>=ub | X[pi[i]]==X[pi[i-1]]) Q[jj] = X[pi[i]] else { g = (ub - P[jj]*W) / (ub - rsum0) Q[jj] = X[pi[i-1]]*g + X[pi[i]]*(1-g) } j++ if (j>J) break jj = pj[j] if (i<I & rsum < P[jj]*W) break } if (j>J) break } return(Q) } end
*! *! mm_median.mata *! version 1.0.3, Ben Jann, 22jun2006 version 9.2 mata:
real rowvector mm_median(real matrix X, | real colvector w) { if (args()==1) w = 1 return(mm_quantile(X, w, .5)) } end
*! *! mm_iqrange.mata *! version 1.0.4, Ben Jann, 22jun2006 version 9.2 mata:
real matrix mm_iqrange(real matrix X, | real colvector w, real scalar altdef) { real matrix q
if (args()<2) w = 1 if (args()<3) altdef = 0 q = mm_quantile(X, w, (.25 \ .75), altdef) return(q[2,]-q[1,]) } end
*! *! mm_ecdf.mata *! version 1.0.7 03aug2007 Ben Jann version 9.2 mata:
real matrix mm_ecdf(real matrix X, | real colvector w, real scalar mid) { real scalar W
if (args()<2) w = 1 if (args()<3) mid = 0
return(mm_ranks(X, w, 3, mid, 1)) }
end
*! *! mm_relrank.mata *! version 1.0.7 03aug2007 Ben Jann version 9.2 mata:
real matrix mm_relrank( real matrix X, real colvector w, real matrix Q, | real scalar mid) { real rowvector result real scalar c, cX, cQ, r, i
if (args()<4) mid = 0
if (cols(X)==1 & cols(Q)!=1 & rows(Q)==1) return(mm_relrank(X, w, Q', mid)' > ) if (rows(w)!=1 & rows(w)!=rows(X)) _error(3200) r = rows(Q) c = max(((cX=cols(X)), (cQ=cols(Q)))) if (cX!=1 & cX<c) _error(3200) if (cQ!=1 & cQ<c) _error(3200) if (rows(X)==0 | r==0 | c==0) return(J(r,c,.)) if (c==1) return(_mm_relrank(X, w, Q, mid)) result = J(r, c, .) if (cQ==1) for (i=1; i<=c; i++) result[,i] = _mm_relrank(X[,i], w, Q, mid) else if (cX==1) for (i=1; i<=c; i++) result[,i] = _mm_relrank(X, w, Q[,i], mid) else for (i=1; i<=c; i++) result[,i] = _mm_relrank(X[,i], w, Q[,i], mid) return(result) }
real colvector _mm_relrank( real colvector X, real colvector w, real colvector Q, real scalar mid) { real colvector P, cdf, pi, pj, sortx real scalar i, I, j, J
I = rows(X) pi = order((X,(1::I)),(1,2)) //stable sort order sortx = X[pi] cdf = __mm_ranks(sortx, (rows(w)!=1 ? w[pi] : w), 3, 0, 1)
J = rows(Q) pj = order(Q, 1) P = J(J, 1, 0) i = I for(j=J; j>=1; j--) { for (; i>=1; i--) { if (sortx[i]<=Q[pj[j]]) break } if (i<1) break P[pj[j]] = cdf[i] } if (mid==0) return(P)
cdf = _mm_ranks_mid(sortx, cdf, mid) i = 1 for (j=1; j<=J; j++) { for (; i<I; i++) { if (sortx[i]>=Q[pj[j]]) break } if (sortx[i]==Q[pj[j]]) P[pj[j]] = cdf[i] } return(P) } end
*! *! mm_ranks.mata *! version 1.0.6 03aug2007 Ben Jann version 9.2 mata:
real matrix mm_ranks(real matrix X, | real colvector w, real scalar method, real scalar mid, real scalar norm) { real rowvector result real scalar c, r, i
if (args()<2) w = 1 if (args()<3) method = 0 if (args()<4) mid = 0 if (args()<5) norm = 0 c = cols(X) r = rows(X) if (rows(w)!=1 & rows(w)!=r) _error(3200) if (c<1) return(J(r,0,.)) if (c==1) return(_mm_ranks(X, w, method, mid, norm)) result = J(r, c, .) for (i=1; i<=c; i++) { result[,i] = _mm_ranks(X[,i], w, method, mid, norm) } return(result) }
real colvector _mm_ranks(real colvector x, real colvector w, real scalar method, real scalar mid, real scalar norm) { real scalar I real colvector p
I = rows(x) if (I==0) return(J(0,1,.)) if (method==4 & rows(w)!=1) p = order((x,w),(1,2)) // -set sortseed- for re > producibility else if (method==0 | method==4) p = order(x,1) // -set sortseed- for reprod > ucibility else if (method==1 & rows(w)!=1) p = order((x,w,(1::I)),(1,2,3)) // stable > sort order else p = order((x,(1::I)),(1,2)) // stable sort order return(__mm_ranks(x[p],(rows(w)!=1 ? w[p] : w), method, mid, norm)[invorder > (p)]) }
real colvector __mm_ranks(real colvector x, real colvector w, real scalar method, real scalar mid, real scalar norm) // sorted input assum > ed { return(_mm_ranks_mid(x, ___mm_ranks(x, w, method, norm), mid)) }
real colvector ___mm_ranks(real colvector x, real colvector w, real scalar method, real scalar norm) // sorted input assum > ed { real scalar i, I, i0, i1 real colvector ranks, R pointer scalar wR
I = rows(x) if (I==0) return(J(0,1,.)) if (rows(w)!=1) { if (missing(w)) return(J(I,1,.)) ranks = mm_colrunsum(w) } else ranks = (1::I) * w if (norm) ranks = ranks / ranks[I] if (method==0 | method==4) return(ranks) // no treatment of ties if (method==3) { // ties lowest for(i=I-1; i>=1; i--) { if (x[i]==x[i+1]) ranks[i] = ranks[i+1] } if (method==3) return(ranks) } if (method==1) { // ties highest for(i=2; i<=I; i++) { if (x[i]==x[i-1]) ranks[i] = ranks[i-1] } return(ranks) } if (method!=2) _error(3498,"ties = " + strofreal(method) + " not allowed") i0 = i1 = 1 // ties mean if (rows(w)==1) wR = &1 else wR = &R for(i=2; i<=I; i++) { if (x[i0]==x[i]) { i1++ if (i<I) continue } if (i0<i1) { R = (i0 \ i1) ranks[|R|] = J(i1-i0+1,1,1) :* mean(ranks[|R|], w[|*wR|]) } i0 = i1 = i } return(ranks) }
real colvector _mm_ranks_mid(real colvector x, real colvector ranks, real scalar mid) { real scalar r, lr, i
if (rows(ranks)<1 | mid==0) return(ranks) r = ranks[1] ranks[1] = r/2 lr = r for(i=2; i<=rows(ranks); i++) { if (x[i]==x[i-1]) ranks[i] = ranks[i-1] else { r = ranks[i] ranks[i] = lr + (r-lr)/2 lr = r } } return(ranks) } end
*! *! mm_freq.mata *! version 1.0.2, Ben Jann, 16may2007 version 9.1 mata:
real colvector mm_freq(transmorphic matrix x, | real colvector w, transmorphic matrix levels) { real colvector p
if (args()<2) w = 1 if (args()<3) levels = . if (cols(x)==0) return(_mm_freq(x, w, levels)) if (rows(w)==1) return(_mm_freq(sort(x,1..cols(x)), w, levels)) p = order(x,1..cols(x)) return(_mm_freq(x[p,], w[p,], levels)) }
real colvector _mm_freq(transmorphic matrix x, | real colvector w, transmorphic matrix levels) { real scalar i, j, l real colvector result
if (args()<2) w = 1 if (args()<3) levels = . if (rows(w)!=1 & rows(w)!=rows(x)) _error(3200) if (levels==.) levels = _mm_uniqrows(x) if (rows(x)==0) return(J(0,1, .)) l = rows(levels) result = J(l,1,0) j = 1 for (i=1; i<=rows(x); i++) { for (;j<=l;j++) { if (x[i,]==levels[j,]) break } if (j>l) break result[j] = result[j] + (rows(w)!=1 ? w[i] : w) } return(result) }
real colvector mm_freq2(transmorphic matrix x, | real colvector w) { real colvector p
if (args()<2) w = 1 if (cols(x)==0) return(_mm_freq2(x, w)) p = order(x,1..cols(x)) if (rows(w)==1) return(_mm_freq2(x[p,],w)[invorder(p)]) return(_mm_freq2(x[p,],w[p,])[invorder(p)]) }
real colvector _mm_freq2(transmorphic matrix x, | real colvector w) { real scalar i, j real colvector result
if (args()<2) w = 1 if (rows(w)!=1 & rows(w)!=rows(x)) _error(3200) if (rows(x)==0) return(J(0, 1, .)) result = J(rows(x),1,.) j = 1 for (i=2; i<=rows(x); i++) { if (x[i,]!=x[i-1,]) { result[|j \ i-1|] = J(i-j, 1, (rows(w)==1 ? (i-j)*w : sum(w[|j \ i-1|]))) j = i } } result[|j \ i-1|] = J(i-j, 1, (rows(w)==1 ? (i-j)*w : sum(w[|j \ i-1|]))) return(result) }
transmorphic matrix _mm_uniqrows( // uniqrows() for sorted X transmorphic matrix x) { real scalar i, j, n, ns transmorphic matrix res
if (rows(x)==0) return(J(0,cols(x), missingof(x))) if (cols(x)==0) return(J(1,0, missingof(x)))
ns = 1 n = rows(x) for (i=2;i<=n;i++) { if (x[i-1,]!=x[i,]) ns++ } res = J(ns, cols(x), x[1,1]) res[1,] = x[1,] for (i=j=2;i<=n;i++) { if (x[i-1,]!=x[i,]) res[j++,] = x[i,] } return(res) }
end
*! *! mm_histogram.mata *! version 1.0.3, Ben Jann, 22jun2006 version 9.0 mata:
real matrix mm_histogram( real colvector x, // data | real colvector w, // weights real colvector g, // grid (interval borders) (default: 10 intervals) real scalar dir) // 0 right inclusive (default), else left inclusive { real matrix H real scalar i
if (args()<2) w = 1 if (args()<3) g = rangen(min(x), max(x), 11) if (args()<4) dir = 0
H = J(rows(g)-1, 3, 0) for (i=1; i<=rows(H); i++) { H[i,2] = g[i+1]-g[i] // width of bins H[i,1] = g[i] + H[i,2]/2 // bin midpoints } H[.,3] = mm_exactbin(x, w, g, dir) :/ (mm_nobs(x,w) * H[.,2]) return(H) } end
*! *! mm_collapse.mata *! version 1.0.1 18feb2009 Ben Jann version 9.2
local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"
mata:
real matrix mm_collapse( real matrix X, // data real colvector w, // weights (or 1) real colvector ID, // subgroup IDs | pointer(function) scalar f, // function; default: mean() `opts') // arguments to pass through to f() { real colvector p pointer scalar Xs, IDs, ws
if ((rows(X)!=rows(ID)) | (rows(w)!=1 & rows(X)!=rows(w))) _error(3200) p = order(ID, 1) if (isfleeting(X)) Xs = &(X = X[p,]) else Xs = &(X[p,]) if (rows(w)==1) ws = &w else if (isfleeting(w)) ws = &(w = w[p]) else ws = &(w[p]) if (isfleeting(ID)) IDs = &(ID = ID[p]) else IDs = &(ID[p])
if (args()==3) return(_mm_collapse(*Xs, *ws, *IDs)) else if (args()==4) return(_mm_collapse(*Xs, *ws, *IDs, f)) else if (args()==5) return(_mm_collapse(*Xs, *ws, *IDs, f, o1)) else if (args()==6) return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2)) else if (args()==7) return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3)) else if (args()==8) return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3, o4) > ) else if (args()==9) return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3, o4, > o5)) else if (args()==10) return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3, o4, > o5, o6)) else if (args()==11) return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3, o4, > o5, o6, o7)) else if (args()==12) return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3, o4, > o5, o6, o7, o8)) else if (args()==13) return(_mm_collapse(*Xs, *ws, *IDs, f, o1, o2, o3, o4, > o5, o6, o7, o8, o9)) else if (args()==14) return(_mm_collapse(*Xs, *ws, *IDs, f, `opts')) }
real matrix _mm_collapse( real matrix X, // data real colvector w, // weights (or 1) real colvector ID, // subgroup IDs | pointer(function) scalar f0, // function; default: mean() `opts') // arguments to pass through to f() { real scalar i, j, c, a, b, ww real matrix info, res, key pointer(function) scalar f transmorphic setup
if (rows(X)!=rows(ID)) _error(3200) ww = rows(w)!=1 if (ww & rows(w)!=rows(X)) _error(3200) if (rows(X)<1) return(J(0, cols(X)+1, missingof(X)))
f = args()<4 ? &mean() : f0 setup = mm_callf_setup(f, args()-4, `opts')
info = _mm_panels(ID) key = J(rows(info), 1, missingof(ID)) res = J(rows(info), c=cols(X), missingof(X)) b = 0 for (i=1; i<=rows(info); i++) { a = b + 1 b = b + info[i] key[i] = ID[a] for (j=1; j<=c; j++) { res[i, j] = mm_callf(setup, X[|a,j \ b,j|], ww ? w[|a \ b|] : w) } } return(key, res) }
end
*! *! mm_gini.mata *! version 1.0.3 26mar2008 Ben Jann version 9.2 mata:
real rowvector mm_gini(real matrix X, | real colvector w) { real rowvector result real scalar c, i
if (args()==1) w = 1 c = cols(X) if (c<1) return(J(1,0,.)) if (c==1) return(_mm_gini(X, w)) result = J(1, c, .) for (i=1; i<=c; i++) { result[1, i] = _mm_gini(X[,i], w) } return(result) }
real scalar _mm_gini(real colvector x, real colvector w) { real matrix mv
if (rows(x)<1) return(.) mv = mm_meanvariance0((x, mm_ranks(x, w, 3, 1, 1)), w) return(mv[3,1] * 2 / mv[1,1]) }
end
*! *! mm_colvar.mata *! version 1.0.1, Ben Jann, 22jun2006 version 9.0 mata:
real rowvector mm_colvar(real matrix X, |real colvector w) { real rowvector result real scalar j, k
if (args()==1) w = 1
result = J(1, k=cols(X), .) for (j=1; j<=k; j++) { result[j] = variance(X[,j], w) } return(result) }
real matrix mm_meancolvar(real matrix X, |real colvector w) { real matrix result real scalar j, k
if (args()==1) w = 1
result = J(2, k=cols(X), .) for (j=1; j<=k; j++) { result[.,j] = meanvariance(X[,j], w) } return(result) }
end
*! *! mm_variance0.mata *! version 1.0.1, Ben Jann, 22jun2006 version 9.0 mata:
real matrix mm_variance0(real matrix X, |real colvector w) { real rowvector CP real rowvector means real scalar n
if (args()==1) w = 1
CP = quadcross(w,0, X,1) n = cols(CP) means = CP[|1\n-1|] :/ CP[n] if (missing(CP) | CP[n]==0) return(J(cols(X), cols(X), .)) return(crossdev(X,0,means, w, X,0,means) :/ CP[n]) }
real matrix mm_meanvariance0(real matrix X, |real colvector w) { real rowvector CP real rowvector means real scalar n
if (args()==1) w = 1
CP = quadcross(w,0, X,1) n = cols(CP) means = CP[|1\n-1|] :/ CP[n] if (missing(means)) return(means \ J(cols(X),cols(X),.)) return(means \ crossdev(X,0,means, w, X,0,means) :/ CP[n]) }
end
*! *! mm_mse.mata *! version 1.0.0, Ben Jann, 22jun2006 version 9.0 mata:
real matrix mm_mse(real matrix X, real colvector w, real rowvector mu) { real rowvector CP
CP = cross(w,0, X,1) if (missing(CP) | CP[cols(CP)]==0) return(J(cols(X), cols(X), .)) return(mm_sse(X, w, mu) :/ CP[cols(CP)]) }
real matrix mm_sse(real matrix X, real colvector w, real rowvector mu) { real rowvector CP
CP = cross(w,0, X,1) if (missing(CP) | CP[cols(CP)]==0) return(J(cols(X), cols(X), .)) return(crossdev(X,0,mu, w, X,0,mu)) }
real rowvector mm_colmse(real matrix X, real colvector w, real rowvector mu) { real rowvector result real scalar j, k
result = J(1, k=cols(X), .) for (j=1; j<=k; j++) { result[j] = mm_mse(X[,j], w, mu[j]) } return(result) }
real rowvector mm_colsse(real matrix X, real colvector w, real rowvector mu) { real rowvector result real scalar j, k
result = J(1, k=cols(X), .) for (j=1; j<=k; j++) { result[j] = mm_sse(X[,j], w, mu[j]) } return(result) }
end
*! *! mm_sample.mata *! version 1.0.1, Ben Jann, 14apr2006 version 9.1 mata:
real colvector mm_sample( real colvector n, real matrix strata, | real colvector cluster, real colvector w, real scalar wor, real scalar count, real scalar fast) { real scalar i, b, e, n0, n1, N, cb, ce real colvector s, si, ci, wi, nn, R pointer scalar f
if (args()<3) cluster=. if (args()<4) w=1 if (args()<5) wor=0 if (args()<6) count=0 if (args()<7) fast=0 if (cols(strata)<1) _error(3200) if (fast==0) { if (rows(strata)>1 | (cluster==. & rows(w)==1)) { if (missing(strata)) _error(3351, "'strata' has missing values") } } if (rows(n)<1) _error(3200)
//choose sampling function and check w if (rows(w)==1) { if (wor) f = &mm_srswor() else f = &mm_srswr() } else { if (cluster!=. & rows(w)!=rows(cluster)) _error(3200, "rows('w') unequal number of clusters") if (fast==0) { if (missing(w)) _error(3351, "'w' has missing values") if (colsum(w:<0)) _error(3498, "'w' has negative values") if (colsum(w)<=0) _error(3498, "sum('w') is zero") } if (wor) f = &mm_upswor() else f = &mm_upswr() }
//simple sample, unstratified if (rows(strata)==1) { N = strata[1,1] if (cluster==.) { if (N>=.) N = rows(w) else if (rows(w)!=1 & N!=rows(w)) _error(3200, "rows('w') unequal population size") return((*f)(n, (rows(w)==1 ? N : w), count)) }
//cluster sample, unstratified if (rows(cluster)<1) _error(3200) if (N>=.) N = colsum(cluster) else if (fast==0) { if (N!=colsum(cluster)) _error(3200, "sum of cluster sizes unequal population size") } s = (*f)(n, (rows(w)==1 ? rows(cluster) : w), count) return(_mm_expandclusters(s, cluster, count, N)) }
//stratified sample if (rows(strata)<1) _error(3200) if (cluster!=. & cols(strata)<2) _error(3200) //-sample sizes if (rows(n)==1) { if (n>=.) { if (cluster==.) nn = strata[.,1] else nn = strata[.,2] } else nn = J(rows(strata),1,n) } else { if (rows(n)!=rows(strata)) _error(3200) nn = n for (i=1;i<=rows(nn);i++) { if (nn[i]>=.) { if (cluster==.) nn[i] = strata[i,1] else nn[i] = strata[i,2] } } } //-prepare sample vector if (count) s = J(colsum(strata[.,1]),1,.) else s = J(colsum(nn),1,.) //-draw samples within strata if (count==0) n0 = 1 e = 0 ce = 0 for (i=1;i<=rows(strata);i++) { b = e + 1 e = e + strata[i,1] //--simple sample if (cluster==.) { si = (*f)(nn[i], (rows(w)==1 ? strata[i,1] : w[|b \ e|]), count) } //--cluster sample else { cb = ce + 1 ce = ce + strata[i,2] ci = cluster[|cb \ ce|] wi = (rows(w)==1 ? rows(ci) : w[|cb \ ce|]) if (fast==0) { if (strata[i,1]!=colsum(ci)) _error(3200, "sum of cluster sizes unequal size of stratum") } si = (*f)(nn[i], wi, count) if (count) si = _mm_expandclusters(si, ci, 1, strata[i,1]) } //--add subsample to sample vector if (count) R = (b \ e) else { n1 = n0 + rows(si) - 1 R = (n0 \ n1) if (cluster==.) si = si :+ (b-1) else si = si :+ (cb-1) n0 = n1 + 1 } if (R[1]<=R[2]) s[|R|] = si } if (count==0&cluster!=.) s = _mm_expandclusters(s, cluster) return(s) }
real colvector _mm_expandclusters(real colvector s0, real colvector cluster, | real scalar count, real scalar N) { real scalar i, e, b, n0, n1 real colvector s, eclust
if (args()<3) count = 0 if (count) { s = J(N,1,.) e = 0 for (i=1;i<=rows(cluster);i++) { b = e + 1 e = e + cluster[i] s[|b \ e|] = J(e-b+1, 1, s0[i]) } return(s) } s = J(colsum(cluster[s0,]),1,.) eclust = mm_colrunsum(cluster) n0 = 1 for (i=1;i<=rows(s0);i++) { e = eclust[s0[i]] b = e - cluster[s0[i]] + 1 n1 = n0 + e-b s[|n0 \ n1|] = (b::e) n0 = n1+1 } return(s) }
end
*! *! mm_srswr.mata *! version 1.0.0, Ben Jann, 27mar2006 version 9.1 mata:
// simple random sample, with replacement real colvector mm_srswr(real scalar n, real scalar N, | real scalar count) { real colvector res, u real scalar i, nn
// check args if (args()<3) count = 0 if (N>=.) _error(3351) if (N<1) _error(3300) if (n>=.) nn = N else nn = n
// no sampling if (N==1) { if (count) return(nn) return(J(nn,1,1)) }
// draw sample u = ceil(uniform(nn,1)*N) if (count==0) return(u) res = J(N,1,0) for (i=1;i<=nn;i++) { res[u[i]] = res[u[i]] + 1 } return(res) }
end
*! *! mm_srswor.mata *! version 1.0.0, Ben Jann, 29mar2006 version 9.1 mata:
// simple random sample, without replacement real colvector mm_srswor(real scalar n, real scalar N, | real scalar count) { real colvector res, u real scalar i, nn
// check args if (args()<3) count = 0 if (N>=.) _error(3351) if (N<1) _error(3300) if (n<0) _error(3300) if (n>=.) nn = N else nn = n if (N<nn) _error(3300, "n may not be larger than N")
// no sampling if (N==1) { if (count) return(nn) return(J(nn,1,1)) } if (nn==0) { if (count) return(J(N,1,0)) return(J(0,1,.)) }
// draw sample u = mm_unorder2(N)[|1 \ nn|] //=> stable results even for very large N if (count==0) return(u) res = J(N,1,0) for (i=1;i<=nn;i++) { res[u[i]] = 1 } return(res) } end
*! *! mm_upswr.mata *! version 1.0.0, Ben Jann, 01apr2006 version 9.1 mata:
// unequal probability sampling, with replacement real colvector mm_upswr(real scalar n, real colvector w, | real scalar count) { real colvector ub, res, u real scalar i, j, nn
// check args if (args()<3) count = 0 if (rows(w)<1) _error(3498, "no cases") if (n>=.) nn = rows(w) else nn = n
// no sampling if (rows(w)==1) { if (count) return(nn) return(J(nn,1,1)) }
// draw sample ub = mm_colrunsum(w) ub = ub/ub[rows(ub)] u = sort(uniform(nn,1),1) j=1 if (count) { res = J(rows(w),1,0) for (i=1;i<=nn;i++) { while (u[i]>ub[j]) j++ res[j] = res[j]+1 } return(res) } res = J(nn,1,0) for (i=1;i<=nn;i++) { while (u[i]>ub[j]) j++ res[i] = j } return(mm_jumble2(res)) //=> stable results even for very large n }
end
*! *! mm_upswor.mata *! version 1.0.1, Ben Jann, 12jun2007 version 9.1 mata:
// unequal probability sampling, without replacement real colvector mm_upswor(real scalar n, real colvector w, | real scalar count, real scalar nowarn) { real colvector p, ub, res real scalar i, j, nn, u, z
// check args if (args()<4) nowarn = 0 if (args()<3) count = 0 if (rows(w)<1) _error(3498, "no cases") if (n<0) _error(3300) if (n>=.) nn = rows(w) else nn = n if (nowarn==0) { if (rows(w)<nn) _error(3300, "n may not be larger than rows(w)") }
// no sampling if (rows(w)==1) { if (count) return(nn) return(J(nn,1,1)) }
// draw sample p = mm_unorder2(rows(w)) //=> stable results even for very large N ub = w[p] / colsum(w) * nn if (nowarn==0) { z = colsum(ub:>1) if (z==1) _error(3300, "1 case has w_i*n/sum(w)>1") if (z>0) _error(3300, strofreal(z) + " cases have w_i*n/sum(w)>1") } ub = mm_colrunsum(ub) u = uniform(1,1) j=1 if (count) { res = J(rows(w),1,0) for (i=1;i<=nn;i++) { while (u>ub[j]) j++ res[j] = res[j]+1 u = u + 1 } return(res[invorder(p)]) } res = J(nn,1,0) for (i=1;i<=nn;i++) { while (u>ub[j]) j++ res[i] = p[j] u = u + 1 } return(res) }
end
*! *! mm_jk.mata *! version 1.0.0, Ben Jann, 04jul2006 version 9.2
local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"
mata:
struct mm_jkstats { real colvector reps, ndrop, fpc real rowvector stat real matrix rstat }
// note: reps is total n. of replications // ndrop is n. of unsuccessful replications // rows(rstat) = reps-ndrop is n. of successful replications
struct mm_jkstats scalar mm_jk( pointer(real rowvector function) scalar f, real matrix X, | real colvector w, real scalar nodots, // 0 real colvector strata, // . real colvector cluster, // . real colvector subpop, // . real colvector fpc, // . real matrix stat, // . `opts') // opts to pass to f { real scalar i, ii, i0, i1, h, h0, h1, j real colvector C, wjk, wjki real matrix S transmorphic fsetup struct mm_jkstats scalar jk pointer scalar ws
if (args()<3) w = 1 if (args()<4) nodots = 0 if (args()<7) subpop = . if (args()<8) fpc = . if (args()<9) stat = .
if (subpop!=. & rows(subpop)>0) ws = &(w :* (subpop:!=0)) else ws = &w if (rows(*ws)==1) ws = &J(rows(X), 1, *ws)
fsetup = mm_callf_setup(f, args()-9, `opts')
if (stat==.) jk.stat = mm_callf(fsetup, X, *ws)[1,] else jk.stat = stat[1,]
mm_panels(strata, S=J(1,2,rows(X)), cluster, C=.) if (anyof(S, 1)) _error(460, "singleton cluster detected")
if (nodots==0) { printf("{txt}\nJackknife replications ({res}%g{txt})\n", colsum(S[,2])) display("{txt}{hline 4}{c +}{hline 3} 1 " + "{hline 3}{c +}{hline 3} 2 " + "{hline 3}{c +}{hline 3} 3 " + "{hline 3}{c +}{hline 3} 4 " + "{hline 3}{c +}{hline 3} 5 ") } jk.rstat = J(colsum(S[,2]), cols(jk.stat), .) jk.ndrop = J(rows(S), 1, 0) jk.fpc = J(rows(S), 1, 0) h1 = i1 = 0 i = ii = 1 for (h=1; h<=rows(S); h++) { // for each stratum h0 = h1 + 1 h1 = h1 + S[h,1] if (fpc!=. & rows(fpc)>0) jk.fpc[h] = fpc[h1] wjk = *ws wjk[| h0 \ h1 |] = wjk[| h0 \ h1 |] * (S[h,2] / (S[h,2] - 1)) for (j=1; j<=S[h,2]; j++) { // for each PSU i0 = i1 + 1 i1 = i1 + (C!=. ? C[i] : 1) wjki = wjk wjki[| i0 \ i1 |] = J(i1-i0+1,1,0) jk.rstat[ii,] = mm_callf(fsetup, X, wjki)[1,] if (missing(jk.rstat[ii,])) { jk.ndrop[h] = jk.ndrop[h] + 1 if (nodots==0) printf("{err}x{txt}") } else { ii++ if (nodots==0) printf(".") } if (nodots==0) { if (mod(i,50)==0) printf(" %5.0f\n",i) displayflush() } i++ } } if (nodots==0 & mod(i-1,50)) display("") if (ii<i) { if (ii>1) jk.rstat = jk.rstat[|1,1 \ ii-1,.|] else jk.rstat = J(0, cols(jk.rstat), .) } jk.reps = S[,2] return(jk) }
real matrix mm_jk_report( struct mm_jkstats scalar jk, // jackknife estimates | string vector what, // b,pseudo,mean,bias,v,se,ci real scalar level, // real scalar mse, // !=0 => use mse formula real vector fpc0) // sampling fraction (for each stratum) { real scalar theta, pseudo, mean, bias, se, v, ci, l, i, alpha, tval real rowvector jkse, jkv, jkmean real colvector fpc, nh, mh real matrix res, jkpseudo
if (args()<2) what = "se" if (args()<3) level = st_numscalar("c(level)") if (args()<4) mse = 0
theta = pseudo = mean = bias = v = se = ci = l = 0 for (i=1; i<=length(what); i++) { if (what[i]=="theta" | what[i]=="b") theta = (l = l + 1) else if (what[i]=="pseudo") pseudo = (l = l + rows(jk.rstat)) else if (what[i]=="mean" ) mean = (l = l + 1) else if (what[i]=="bias" ) bias = (l = l + 1) else if (what[i]=="v" ) v = (l = l + cols(jk.rstat)) else if (what[i]=="se" ) se = (l = l + 1) else if (what[i]=="ci" ) ci = (l = l + 2) else _error(3498, "'" + what[i] + "' not allowed") }
res = J(l, cols(jk.rstat), .) nh = jk.reps - jk.ndrop fpc = fpc0 if (cols(fpc)!=1) fpc = fpc' if (rows(fpc)==0 | fpc==.) fpc = jk.fpc if (rows(fpc)==0) fpc = J(rows(jk.reps),1,0) // observed parameter vector if (theta) res[theta,] = jk.stat // pseudo values if (pseudo | (mse==0 & (mean | bias))) { jkpseudo = jk.rstat + mm_expand(nh, nh, 1, 1):*(jk.stat :- jk.rstat) if (pseudo) res[|pseudo-rows(jk.rstat)+1,1 \ pseudo,.|] = jkpseudo } // jackknife mean if (mean | bias) { if (mse) jkmean = mean(jk.rstat, 1) else jkmean = mean(jkpseudo, 1) if (mean) res[mean,] = jkmean } // bias if (bias) res[bias,] = jkmean - jk.stat // variance if (v) { mh = mm_expand((1:-fpc) :* (nh :- 1) :/ nh, nh, 1, 1) if (mse) jkv = mm_sse(jk.rstat, mh, jk.stat) else jkv = mm_sse(jk.rstat, mh, mean(jk.rstat,1)) res[|v-cols(res)+1,1 \ v,.|] = jkv } // standard errors if (se | ci) { if (v) jkse = sqrt(diagonal(jkv))' else { mh = mm_expand((1:-fpc) :* (nh :- 1) :/ nh, nh, 1, 1) if (mse) jkse = sqrt(mm_colsse(jk.rstat, mh, jk.stat)) else jkse = sqrt(mm_colsse(jk.rstat, mh, mean(jk.rstat,1))) } if (se) res[se,] = jkse } // confidence interval if (ci) { alpha = (100-level)/200 tval = invttail(rows(jk.rstat)-rows(jk.reps), alpha) res[ci-1,] = jk.stat - tval*jkse res[ci,] = jk.stat + tval*jkse } return(res) }
end
*! *! mm_bs.mata *! version 1.0.0, Ben Jann, 10jul2006 version 9.2
local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"
mata:
struct mm_bsstats { real scalar reps, ndrop real rowvector stat, se real matrix rstat, rse }
// note: reps is total n. of replications // ndrop is n. of unsuccessful replications // rows(rstat) = reps-ndrop is n. of successful replications
struct mm_bsstats scalar mm_bs( pointer(real rowvector function) scalar f, real matrix X, | real colvector w, real scalar reps, // 50 real scalar d, // 0 real scalar nodots, // 0 real colvector strata, // . real colvector cluster, // . real matrix stat, // . `opts') // opts to pass to f { real scalar i, ii, hasse real colvector C, N, s real matrix S, res transmorphic setup struct mm_bsstats scalar bs
if (args()<3) w = 1 if (args()<4) reps = 50 if (args()<5) nodots = 0 if (args()<6) d = 0 if (args()<9) stat = .
setup = mm_callf_setup(f, args()-9, `opts')
if (stat==.) { res = mm_callf(setup, X, w) bs.stat = res[1,] if (hasse=(rows(res)>1)) bs.se = res[2,] } else { bs.stat = stat[1,] if (hasse=(rows(stat)>1)) bs.se = stat[2,] }
mm_panels(strata, S=J(1,2,rows(X)), cluster, C=.) if (anyof(S, 1)) _error(460, "singleton cluster detected") if (d==0|d>=.) N = . else N = S[,2] :- d
if (nodots==0) { printf("{txt}\nBootstrap replications ({res}%g{txt})\n", reps) display("{txt}{hline 4}{c +}{hline 3} 1 " + "{hline 3}{c +}{hline 3} 2 " + "{hline 3}{c +}{hline 3} 3 " + "{hline 3}{c +}{hline 3} 4 " + "{hline 3}{c +}{hline 3} 5 ") } bs.rstat = J(reps, cols(bs.stat), .) if (hasse) bs.rse = J(reps, cols(bs.se), .) for (i=ii=1; i<=reps; i++) { s = mm_sample(N, S, C) res = mm_callf(setup, X[s,], rows(w)==1 ? w : w[s,]) if (missing(res)) { if (nodots==0) printf("{err}x{txt}") } else { bs.rstat[ii,] = res[1,] if (hasse) bs.rse[ii,] = res[2,] ii++ if (nodots==0) printf(".") } if (nodots==0) { if (mod(i,50)==0) printf(" %5.0f\n",i) displayflush() } } if (nodots==0 & mod(i-1,50)) display("") if (ii<i) { if (ii>1) { bs.rstat = bs.rstat[|1,1 \ ii-1,.|] if (hasse) bs.rse = bs.rse[|1,1 \ ii-1,.|] } else { bs.rstat = J(0, cols(bs.rstat), .) if (hasse) bs.rse = J(0, cols(bs.rse), .) } } bs.reps = reps bs.ndrop = reps - rows(bs.rstat) return(bs) }
struct mm_bsstats scalar mm_bs2( pointer(real rowvector function) scalar f, real matrix X, | real colvector w, real scalar reps, // 50 real scalar d, // 0 real scalar nodots, // 0 real colvector strata, // . real colvector cluster, // . real matrix stat, // . `opts') // opts to pass to f { real scalar i, ii, hasse real colvector C, N, ws real matrix S, res transmorphic setup struct mm_bsstats scalar bs
if (args()<3) w = 1 if (args()<4) reps = 50 if (args()<5) nodots = 0 if (args()<6) d = 0 if (args()<9) stat = .
setup = mm_callf_setup(f, args()-9, `opts')
if (stat==.) { res = mm_callf(setup, X, w) bs.stat = res[1,] if (hasse=(rows(res)>1)) bs.se = res[2,] } else { bs.stat = stat[1,] if (hasse=(rows(stat)>1)) bs.se = stat[2,] }
mm_panels(strata, S=J(1,2,rows(X)), cluster, C=.) if (anyof(S, 1)) _error(460, "singleton cluster detected") if (d==0|d>=.) N = . else N = S[,2] :- d
if (nodots==0) { printf("{txt}\nBootstrap replications ({res}%g{txt})\n", reps) display("{txt}{hline 4}{c +}{hline 3} 1 " + "{hline 3}{c +}{hline 3} 2 " + "{hline 3}{c +}{hline 3} 3 " + "{hline 3}{c +}{hline 3} 4 " + "{hline 3}{c +}{hline 3} 5 ") } bs.rstat = J(reps, cols(bs.stat), .) if (hasse) bs.rse = J(reps, cols(bs.se), .) for (i=ii=1; i<=reps; i++) { ws = mm_sample(N, S, C, 1, 0, 1):*w res = mm_callf(setup, X, ws) if (missing(res)) { if (nodots==0) printf("{err}x{txt}") } else { bs.rstat[ii,] = res[1,] if (hasse) bs.rse[ii,] = res[2,] ii++ if (nodots==0) printf(".") } if (nodots==0) { if (mod(i,50)==0) printf(" %5.0f\n",i) displayflush() } } if (nodots==0 & mod(i-1,50)) display("") if (ii<i) { if (ii>1) { bs.rstat = bs.rstat[|1,1 \ ii-1,.|] if (hasse) bs.rse = bs.rse[|1,1 \ ii-1,.|] } else { bs.rstat = J(0, cols(bs.rstat), .) if (hasse) bs.rse = J(0, cols(bs.rse), .) } } bs.reps = reps bs.ndrop = reps - rows(bs.rstat) return(bs) }
real matrix mm_bs_report( struct mm_bsstats scalar bs, // boostrap estimates | string vector what, // b,mean,bias,v,se,n,basic,p,bc,bca,t real scalar level, // real scalar mse, // !=0 => use mse formula struct mm_jkstats jk) // jackknife estimates { real scalar theta, mean, bias, v, se, n, basic, p, bc, bca, t, l, i, alpha, tval real rowvector bsv, bsse, bsmean, z0, a real matrix res
if (args()<2) what = "se" if (args()<3) level = st_numscalar("c(level)") if (args()<4) mse = 0
theta = mean = bias = v = se = n = basic = p = bc = bca = t = l = 0 for (i=1; i<=length(what); i++) { if (what[i]=="theta" | what[i]=="b") theta = (l = l + 1) else if (what[i]=="mean" ) mean = (l = l + 1) else if (what[i]=="bias" ) bias = (l = l + 1) else if (what[i]=="v" ) v = (l = l + cols(bs.rstat)) else if (what[i]=="se" ) se = (l = l + 1) else if (what[i]=="n" | what[i]=="ci") n = (l = l + 2) else if (what[i]=="basic") basic = (l = l + 2) else if (what[i]=="p" ) p = (l = l + 2) else if (what[i]=="bc" ) bc = (l = l + 2) else if (what[i]=="bca" ) bca = (l = l + 2) else if (what[i]=="t" ) t = (l = l + 2) else _error(3498, "'" + what[i] + "' not allowed") }
res = J(l, cols(bs.rstat), .) // observed parameter vector if (theta) res[theta,] = bs.stat // mean of bs replicates if (mean | bias) { bsmean = mean(bs.rstat, 1) if (mean) res[mean,] = bsmean } // bias if (bias) res[bias,] = bsmean - bs.stat // variance if (v) { if (mse) bsv = mm_mse(bs.rstat, 1, bs.stat) else bsv = variance(bs.rstat, 1) res[|v-cols(res)+1,1 \ v,.|] = bsv } // standard errors if (se | n) { if (v) bsse = sqrt(diagonal(bsv))' else { if (mse) bsse = sqrt(mm_colmse(bs.rstat, 1, bs.stat)) else bsse = sqrt(mm_colvar(bs.rstat, 1)) } if (se) res[se,] = bsse } // normal ci alpha = (100-level)/200 if (n) { tval = invttail(bs.reps-bs.ndrop-1, alpha) res[n-1,] = bs.stat - tval*bsse res[n,] = bs.stat + tval*bsse } // basic ci if (basic) res[|basic-1,1 \ basic,.|] = 2*bs.stat :- mm_quantile(bs.rstat, 1, ((1-alpha) \ alpha)) // percentile ci if (p) res[|p-1,1 \ p,.|] = mm_quantile(bs.rstat, 1, (alpha \ (1-alpha))) // bias-corrected ci if (bc | bca) z0 = editmissing(invnormal( colsum(bs.rstat:<=bs.stat)/rows(bs.rstat) ), 0) if (bc) { res[|bc-1,1 \ bc,.|] = mm_quantile(bs.rstat, 1, (normal(2*z0 :+ invnormal(alpha)) \ normal(2*z0 :- invnormal(alpha)))) } // bias-corrected and accelerated ci if (bca) { a = colsum((mean(jk.rstat, 1):-jk.rstat):^3) :/ ( 6 * colsum((mean(jk.rstat, 1):-jk.rstat):^2):^1.5 ) res[|bca-1,1 \ bca,.|] = mm_quantile(bs.rstat, 1, (normal(z0 :+ (z0 :+ invnormal(alpha)) :/ (1:-a:*(z0 :+ invnormal(alpha)))) \ normal(z0 :+ (z0 :- invnormal(alpha)) :/ (1:-a:*(z0 :- invnormal(alpha)))))) } // percentile-t ci if (t) { res[|t-1,1 \ t,.|] = bs.stat :- bs.se :* mm_quantile( editmissing((bs.rstat:-bs.stat):/bs.rse, 0), 1, ((1-alpha) \ alpha) ) } return(res) }
end
*! *! mm_rbinomial.mata *! version 1.0.1, Ben Jann, 14apr2006 version 9.1 mata:
real matrix mm_rbinomial(real matrix n, real matrix p) { real matrix x real scalar r, R, c, C transmorphic scalar rn, cn, rp, cp
R = max((rows(n),rows(p))) C = max((cols(n),cols(p))) rn = (rows(n)==1 ? &1 : (rows(n)<R ? _error(3200) : &r)) cn = (cols(n)==1 ? &1 : (cols(n)<C ? _error(3200) : &c)) rp = (rows(p)==1 ? &1 : (rows(p)<R ? _error(3200) : &r)) cp = (cols(p)==1 ? &1 : (cols(p)<C ? _error(3200) : &c)) x = J(R,C,.) for (r=1;r<=R;r++) { for (c=1;c<=C;c++) { x[r,c] = _mm_rbinomial(n[*rn,*cn], p[*rp,*cp]) } } return(x) }
real scalar _mm_rbinomial(real scalar n, real scalar p) { if (p<0|p>1) return(.) if (n<=0|(n-trunc(n))!=0) return(.) //rejection technique if (n<50|p>.03) { return(colsum(uniform(n,1):<p)) } //geometric distribution technique real scalar sum, x sum = 0 for (x=1;x<=n;x++) { sum = sum + ceil((ln(uniform(1,1))/ln(1-p))-1) if (sum > n-x) break } return(x-1) }
end
*! *! mm_cebinomial.mata *! version 1.0.1, Ben Jann, 14apr2006 version 9.1 mata:
real matrix mm_cebinomial(real matrix n, real matrix k, real matrix p) { real matrix x real scalar r, R, c, C transmorphic scalar rn, cn, rk, ck, rp, cp
R = max((rows(n),rows(k),rows(p))) C = max((cols(n),cols(k),cols(p))) rn = (rows(n)==1 ? &1 : (rows(n)<R ? _error(3200) : &r)) cn = (cols(n)==1 ? &1 : (cols(n)<C ? _error(3200) : &c)) rk = (rows(k)==1 ? &1 : (rows(k)<R ? _error(3200) : &r)) ck = (cols(k)==1 ? &1 : (cols(k)<C ? _error(3200) : &c)) rp = (rows(p)==1 ? &1 : (rows(p)<R ? _error(3200) : &r)) cp = (cols(p)==1 ? &1 : (cols(p)<C ? _error(3200) : &c)) x = J(R,C,.) for (r=1;r<=R;r++) { for (c=1;c<=C;c++) { x[r,c] = _mm_cebinomial(n[*rn,*cn], k[*rk,*ck], p[*rp,*cp]) } } return(x) }
real scalar _mm_cebinomial(real scalar n, real scalar k, real scalar p) { real scalar e, e0, i
if (p<0|p>1) return(.) //success prob. out of range if (n<=0|(n-trunc(n))!=0) return(.) //n out of range or non-integer if (k<0|k>n|(k-trunc(k))!=0) return(.) //k out of range or non-integer if (k==0) return(n*p) e = 0 e0 = 0 for (i=1;i<=n-k;i++) { e = e + Binomial(n,k+i,p) if (e==e0) break e0 = e } return(e / Binomial(n,k,p) + k) } end
*! *! mm_benford.mata *! version 1.0.0, Ben Jann, 19jun2007 version 9.2 mata:
real vector mm_benford( real vector digit0, | real scalar pos0, real scalar base0) { real scalar pos, base real vector digit
// check arguments and set defaults if (args()<3 | base0>=.) base = 10 else { base = trunc(base0) if (base <0 | base >10) _error(3300, "base must be in [2,10]") } if (args()<2 | pos0>=.) pos = 1 else { pos = trunc(pos0) if (pos<1) _error(3300, "pos must be 1 or larger") } digit = trunc(digit0) if (any(digit:<0) | any(digit:>=base)) _error(3300, "digit must be in [0,base-1] (base is "+strofreal(base)+")")
// case 1: pos==1 (simple formula) if (pos==1) return(ln(1 :+ 1:/digit)/ln(base))
// case 2: pos > 1 if (cols(digit)==1) return( rowsum(ln(1 :+ 1:/( J(rows(digit),1,base)*(base^(pos-2)..base^(pos-1)-1) :+ digit))) / ln(base > )) else return( colsum(ln(1 :+ 1:/( (base^(pos-2)::base^(pos-1)-1)*J(1,cols(digit),base) :+ digit))) / ln(base > )) }
end
*! *! mm_cauchy.mata version 9.2 // CFBaum 3aug2007 // minor changes by Ben Jann, 09aug2007 // from Wikipedia // http://en.wikipedia.org/wiki/Cauchy_distribution mata: // density function of the Cauchy-Lorentz distribution real matrix mm_cauchyden(real matrix x, real scalar x0, real scalar gamma) { real matrix cauchyden
if (gamma <= 0) { return(x*.) } cauchyden = 1 :/ ( pi() * gamma :* ( 1 :+ ( (x :- x0) :/ gamma ) :^2) ) return(cauchyden) }
// cumulative distribution function of the Cauchy-Lorentz distribution real matrix mm_cauchy(real matrix x, real scalar x0, real scalar gamma) { real matrix cauchy
if (gamma <= 0) { return(x*.) } cauchy = 1 / pi() :* atan( (x :- x0) :/ gamma ) :+ 0.5 return(cauchy) }
// tail probability of the Cauchy-Lorentz distribution // note: cauchy(x,0,1) is the standard Cauchy distribution // with cauchytail(x,0,1) = ttail(1,x) real matrix mm_cauchytail(real matrix x, real scalar x0, real scalar gamma) { return(1 :- mm_cauchy(x, x0, gamma)) }
// inverse cumulative distribution function of the Cauchy-Lorentz distribution real matrix mm_invcauchy(real matrix p, real scalar x0, real scalar gamma) { real matrix invcauchy
if ( gamma <= 0 | min(p) < 0 | max(p) > 1 ) { return(p*.) } invcauchy = x0 :+ gamma :* tan( pi() :* ( p :- 0.5 ) ) :+ ln(p:>=0:&p:<=1) return(invcauchy) } end
*! *! mm_subset.mata *! version 1.0.1, Ben Jann, 02may2007 version 9.2 mata:
/*-SUBSETS---------------------------------------------------------------*/ // Generate all combinations (subsets) one-by-one (lexicographic) // based on Algorithm 5.8 from Reingold et al. (1977) // Usage: // // info = mm_subsetsetup(n,k) // while ((c = mm_subset(info)) != J(0,1,.)) { // ... c ... // } // // The total number of combinations can be computed using // Mata's -comb()- function, i.e. // // Ctot = comb(n, k)
struct mm_subsetinfo { real scalar n, k, i, j, counter, algorithm real colvector x, y }
struct mm_subsetinfo scalar mm_subsetsetup( real scalar n, | real scalar k) { struct mm_subsetinfo scalar s // setup s.n = trunc(n) s.k = trunc(k) if (s.n>=.) _error(3351) if (s.k>=.) s.k = s.n else if (s.k>s.n) _error(3300,"k may not be larger than n") if (s.n<1|s.k<1) s.i = 0 // => return empty set else { s.x = -1 \ (1::s.k) s.i = 1 } s.i = s.i + 1 // offset i by 1 s.k = s.k + 1 // offset k by 1 s.counter = 0 return(s) }
real colvector mm_subset(struct mm_subsetinfo scalar s) { real scalar j real colvector subset
if (s.i==1) return(J(0,1,.)) /*done*/ // get current subset (to be returned at end) subset = s.x[|2 \ s.k|] // generate next subset s.i = s.k while (s.x[s.i] == s.n-s.k+s.i) { s.i = s.i - 1 } s.x[s.i] = s.x[s.i] + 1 for (j=s.i+1;j<=s.k;j++) s.x[j] = s.x[j-1] + 1 /* the following would improve speed but does not work for some reason if (s.i<s.k) s.x[|s.i+1 \ s.k|] = s.x[|s.i \ s.k-1|] :+ 1 */ s.counter = s.counter + 1 return(subset) }
// wrapper for mm_subsetinfo()/mm_subset() real matrix mm_subsets(real scalar n, | real scalar k) { real scalar i real colvector set real matrix res struct mm_subsetinfo scalar s
s = mm_subsetsetup(n, (k<. ? k : n)) if ((set = mm_subset(s)) == J(0,1,.)) return(set) res = J(rows(set), comb(n, rows(set)), .) res[,i=1] = set while ((set = mm_subset(s)) != J(0,1,.)) res[,++i] = set return(res) }
/*-COMPOSITIONS----------------------------------------------------------*/ // Generate all k-way compositions of n, one-by-one // Default algorithm: direct (less elegant but faster) (anti-lexicographic) // Alternate algorithm: indirect via subsets (lexicographic) // // Usage: // // info = mm_compositionsetup(n,k) // while ((c = mm_composition(info)) != J(0,1,.)) { // ... c ... // } // // or: // // info = mm_compositionsetup(n,k,1) // while ((c = mm_composition(info)) != J(0,1,.)) { // ... c ... // } // // mm_ncompositions() returns the total number of compositions
real scalar mm_ncompositions( real scalar n, | real scalar k0) { real scalar k
k = trunc(k0<. ? k0 : n) return(comb(trunc(n)+k-1,k-1)) }
struct mm_subsetinfo scalar mm_compositionsetup( real scalar n, | real scalar k, real scalar alt) // !=0 => alternate algorithm { struct mm_subsetinfo scalar s
if (args()<3) alt = 0 if (alt) { s = _mm_composition2setup(n,k) s.algorithm = 2 return(s) } s = _mm_compositionsetup(n,k) s.algorithm = 1 return(s) }
real colvector mm_composition(struct mm_subsetinfo scalar s) { if (s.algorithm==1) return(_mm_composition(s)) if (s.algorithm==2) return(_mm_composition2(s)) _error(3498, "unknown algorithm") }
// default algorithm
struct mm_subsetinfo scalar _mm_compositionsetup( real scalar n, | real scalar k) { struct mm_subsetinfo scalar s
s.n = trunc(n) s.k = (k<. ? trunc(k) : s.n) if (s.n>=.) _error(3351) if (s.k<1) { s.i = 0 s.x = J(0,1,.) } else if (s.n<1) { s.i = 0 s.x = J(s.k,1,0) } else if (s.k<2) { s.i = 0 s.x = J(s.k,1,s.n) } else { s.x = s.n \ J(s.k-1,1,0) s.i = 1 } s.j = 2 s.counter = 0 return(s) }
real colvector _mm_composition(struct mm_subsetinfo scalar s) { real colvector subset
subset = s.x if (s.i<1) { if (subset!=J(0,1,.)) s.counter = s.counter + 1 s.x = J(0,1,.) return(subset) /*done*/ } while (s.x[s.i]==0 & s.i>1) { // move back s.x[s.i] = s.x[s.j] s.x[s.j] = 0 s.j = s.i-- } if (s.x[s.i]>0) { s.x[s.i] = s.x[s.i]-1 s.x[s.j] = s.x[s.j]+1 if (s.j<s.k) s.i = s.j++ } else { s.i = 0 s.x = J(0,1,.) } s.counter = s.counter + 1 return(subset) }
// alternate algorithm
struct mm_subsetinfo scalar _mm_composition2setup( real scalar n, | real scalar k0) { real scalar k
k = trunc(k0<. ? k0 : n) return(mm_subsetsetup(trunc(n)+k-1,k-1)) }
real colvector _mm_composition2(struct mm_subsetinfo scalar s) { real colvector subset
if (s.i==1) { if (s.k==1) { // length of composition is 1 s.k = 0 s.counter = s.counter + 1 return(s.n) } return(J(0,1,.)) /*done*/ } subset = mm_subset(s) return( (subset \ s.n+1) - (0 \ subset) :- 1 ) }
// wrapper for mm_compositioninfo()/mm_composition() real matrix mm_compositions(real scalar n, | real scalar k, real scalar alt) { real scalar i real colvector set real matrix res struct mm_subsetinfo scalar s
if (args()<3) alt = 0 s = mm_compositionsetup(n, (k<. ? k : n), alt) if ((set = mm_composition(s)) == J(0,1,.)) return(set) res = J(rows(set), mm_ncompositions(n, rows(set)), .) res[,i=1] = set while ((set = mm_composition(s)) != J(0,1,.)) res[,++i] = set return(res) }
/*-RANDOM SUBSETS--------------------------------------------------------*/ // Return random combination (subsets) // Usage: // // for (i=1;i<=reps;i++) { // c = mm_rsubset(n, k) // ... c ... // } // // mm_rsubset() returns jumbled sets; apply sort() if you need // ordered sets, e.g. // // c = sort(mm_rsubset(n, k), 1)
real colvector mm_rsubset( real scalar n, | real scalar k) { if (k<. & k>n) _error(3300, "k may not be larger than n") if (n<1|k<1) return(J(0,1,.)) return(mm_unorder2(n)[|1 \ (k<. ? k : n)|]) }
/*-RANDOM COMPOSITIONS---------------------------------------------------*/ // Return random combination (subsets) and random composition // Usage: // // for (i=1;i<=reps;i++) { // c = mm_rcomposition(n, k) // ... c ... // }
real colvector mm_rcomposition( real scalar n0, | real scalar k0) { real scalar n, k real colvector subset
n = trunc(n0) k = (k0<. ? trunc(k0) : n) if (n>=.) _error(3351) if (k<1) return(J(0,1,.)) if (n<1) return(J(k,1,0)) if (k<2) return(J(k,1,n)) subset = sort(mm_rsubset(n+k-1,k-1), 1) return( (subset \ n+k) - (0 \ subset) :- 1 ) }
/*-PARTITIONS------------------------------------------------------------*/ // Generate all k-way partitions of n, one-by-one // Default algorithm: based on ZS1 algorithm in Zoghbi&Stojmenovic (1998) // (anti-lexicographic) // Alternate algorithm: based on Algorithm 5.12 in Reingold et al. (1977) // (dictionary order) // Usage: // // info = mm_partitionsetup(n,k) // while ((c = mm_partition(info)) != J(0,1,.)) { // ... c ... // } // // or: // // info = mm_partitionsetup(n,k,1) // while ((c = mm_partition(info)) != J(0,1,.)) { // ... c ... // } // // mm_npartitions() returns the total number of compositions, based // on algorithms from // http://home.att.net/~numericana/answer/numbers.htm#partitions
// number of partitions of n with k or fewer addends (slow) real scalar mm_npartitions( real scalar n, | real scalar k) { real scalar u, i, j, m real colvector p, a
m = min((n,k)) if (m==n) return(_mm_npartitions(n)) p = J(n+1,1,1) for (u=2;u<=m;u++) { a = p p = J(n+1,1,0) for (i=0;i<=n;i=i+u) { for (j=i+1;j<=n+1;j++) { p[j] = p[j] + a[j-i] } } } return(p[n+1]) } // total number of partitions of n (fast) real scalar _mm_npartitions( real scalar n) { real scalar i, j, s, k real colvector p
p = J(n+1,1,1) for (i=1;i<=n;i++) { j = 1 ; k = 1 ; s = 0 while (j>0) { j = i - (3*k*k + k)/2 if (j>=0) s = s - (-1)^k * p[j+1] j = i - (3*k*k - k)/2 if (j>=0) s = s - (-1)^k * p[j+1] k = k + 1 } p[i+1] = s } return(p[n+1]) }
struct mm_subsetinfo scalar mm_partitionsetup( real scalar n, | real scalar k, real scalar pad, // zero-padded output real scalar alt) // !=0 => alternate algorithm { struct mm_subsetinfo scalar s
if (args()<3) pad = 0 if (args()<4) alt = 0 if (alt) { s = _mm_partition2setup(n,k) s.algorithm = 3 + (pad!=0) return(s) } s = _mm_partitionsetup(n,k) s.algorithm = 1 + (pad!=0) return(s) }
real colvector mm_partition(struct mm_subsetinfo scalar s) { if (s.algorithm==1) return(_mm_partition(s,0)) if (s.algorithm==2) return(_mm_partition(s,1)) if (s.algorithm==3) return(_mm_partition2(s,0)) if (s.algorithm==4) return(_mm_partition2(s,1)) _error(3498, "unknown algorithm") }
// default algorithm
struct mm_subsetinfo scalar _mm_partitionsetup( real scalar n, | real scalar k) { struct mm_subsetinfo scalar s
s.n = trunc(n) s.k = trunc(k) if (s.n>=.) _error(3351) if (s.k>=.) s.k = s.n // else if (s.k>s.n) _error(3300,"k may not be larger than n") if (s.n<1 | s.k<1 ) s.i = -1 else { s.x = s.n \ J(s.k-1, 1, 1) s.i = 1 } s.j = 1 s.counter = 0 return(s) }
real colvector _mm_partition( // original ZS1 Algorithm does not struct mm_subsetinfo scalar s, // support k-way partitions; changes real scalar pad) // are indicated { real scalar r, t, l real colvector subset
if (s.i==-1) return(J(0,1,.)) /*done*/ subset = s.x[|1 \ s.j|] if (pad) subset = subset \ J(s.k-rows(subset),1,0) l = (s.k<s.n ? s.k : s.n) // added if (s.j>=l) { // added; original stopping rule is (s.x[1]==1) if (s.x[1]<=s.x[l]+1) { // changed s.i = -1 s.counter = s.counter + 1 return(subset) } while (s.x[s.i]<=s.x[l]+1) { // step back loop added s.j = s.j + s.x[s.i] - 1 s.i = s.i - 1 } } if (s.x[s.i]==2) { s.j = s.j + 1 s.x[s.i] = 1 s.x[s.j] = 1 // added s.i = s.i - 1 } else { r = s.x[s.i] - 1 t = s.j - s.i + 1 s.x[s.i] = r while (t>=r) { s.i = s.i + 1 s.x[s.i] = r t = t - r } if (t==0) { s.j = s.i } else { s.j = s.i + 1 if (t==1) s.x[s.i+1] = 1 // added if (t>1) { s.i = s.i + 1 s.x[s.i] = t } } } s.counter = s.counter + 1 return(subset) }
// alternate algorithm
struct mm_subsetinfo scalar _mm_partition2setup( real scalar n, | real scalar k) { real scalar offset struct mm_subsetinfo scalar s
s.n = trunc(n) s.k = trunc(k) if (s.n>=.) _error(3351) if (s.k>=.) s.k = s.n offset = 2 if (s.n<1 | s.k<1 ) s.i = 0 + offset else { s.i = 1 + offset s.x = s.y = J(s.n+offset, 1, .) s.x[-1+offset] = 0 s.y[-1+offset] = 0 s.x[0+offset] = 0 s.y[0+offset] = s.n + 1 s.x[1+offset] = s.n s.y[1+offset] = 1 } s.counter = 0 return(s) }
real colvector _mm_partition2( struct mm_subsetinfo scalar s, real scalar pad) { real scalar skip, sum real colvector subset
skip = 1 while (skip) { if (s.i<=2) return(J(0,1,.)) /*done*/ if ((skip = (sum(s.x[|3 \ s.i|])>s.k))==0) subset = _mm_partition2_expand(s.x[|3 \ s.i|], s.y[|3 \ s.i|], (pad ? s.k : .)) sum = s.x[s.i]*s.y[s.i] if (s.x[s.i]==1) { s.i = s.i - 1 sum = sum + s.x[s.i]*s.y[s.i] } if (s.y[s.i-1]==s.y[s.i]+1) { s.i = s.i - 1 s.x[s.i] = s.x[s.i] + 1 } else { s.y[s.i] = s.y[s.i] + 1 s.x[s.i] = 1 } if (sum>s.y[s.i]) { s.y[s.i+1] = 1 s.x[s.i+1] = sum - s.y[s.i] s.i = s.i + 1 } } s.counter = s.counter + 1 return(subset) }
real colvector _mm_partition2_expand( real vector m, real vector p, real scalar k) { real scalar i, j, l real colvector res
if (k<.) res = J(k,1,0) else res = J(sum(m),1,0) l = 0 for (i=1;i<=rows(p);i++) { for (j=1;j<=m[i];j++) { res[++l] = p[i] } } return(res) }
// wrapper for mm_partitioninfo()/mm_partition() real matrix mm_partitions(real scalar n, | real scalar k, real scalar alt) { real scalar i real colvector set real matrix res struct mm_subsetinfo scalar s
if (args()<3) alt = 0 s = mm_partitionsetup(n, (k<. ? k : n), 1, alt) if ((set = mm_partition(s)) == J(0,1,.)) return(set) res = J(rows(set), mm_npartitions(n, rows(set)), .) res[,i=1] = set while ((set = mm_partition(s)) != J(0,1,.)) res[,++i] = set return(res) }
end
*! *! mm_mgof.mata *! version 1.0.7 Ben Jann 12jan2008 * mm_mgof: goodness-of-fit tests for multinomial data
version 9.2 mata:
real matrix mm_mgof( real colvector f, // observed count | real colvector h0, // expected distribution string scalar m0, // method: "approx" (default), "mc", "rc", or "ee" string vector s0, // test-stats: "x2", "lr", "cr", "mlnp", "ksmirnov" real scalar lambda, // lambda for Cressie-Read statistic real scalar arg1, // nfit for "approx"; reps for "mc" real scalar dots, // display dots with mc/ee real scalar arg2) // will be replaced; used by mgof_ee to return counter { real scalar i, n, hsum string scalar m, s string vector stats real vector p pointer(real colvector) scalar h pointer(real scalar function) vector s1
// parsing and defaults if (args()<7) dots = 0 m = mm_strexpand(m0, ("approx", "mc", "rc", "ee"), "approx", 0) stats = ("x2", "lr", "cr") if (m!="approx") stats = (stats, "mlnp", "ksmirnov") if (length(s0)<1) s = "x2" else { s = J(length(s0),1,"") for (i=1;i<=length(s0);i++) { s[i] = mm_strexpand(s0[i], stats, "x2", 1) } } if (m!="approx" & lambda<0) _error(3498, "lambda<0 not allowed with "+m)
// check data / prepare null distribution if (missing(f)) _error(3351, "f has missing values") if (any(f:<0)) _error(3498, "f has negative values") if (any(f:<0)) _error(3498, "f has negative values") if (m=="ee" | m=="rc" | (m=="mc" & anyof(s, "mlnp"))) { if (f!=trunc(f)) _error(3498, "non-integer f not allowed in this contex > t") } if (anyof(stats,"cr") & lambda<0 & anyof(f,0)) _error(3498, "some f are 0: lambda<0 not allowed") if (args()<2) h0 = 1 if (rows(f)!=rows(h0) & rows(h0)!=1) _error(3200) if (missing(h0)) _error(3351, "h has missing values") if (any(h0:<=0)) _error(3498, "some h are negative or 0") n = colsum(f) hsum = colsum(h0) if (n!=hsum) h = &(h0*n/hsum) else h = &h0
// return (0,1) if n. of obs is 0 or n. of categories < 2 if ( n==0 | rows(f)<2 ) return(J(rows(s),1,0),J(rows(s),1,1))
// gather statistics functions s1 = J(rows(s),1,NULL) for (i=1;i<=length(s);i++) { if (s[i]=="x2") s1[i] = &_mm_mgof_x2() else if (s[i]=="lr") s1[i] = &_mm_mgof_lr() else if (s[i]=="cr") s1[i] = &_mm_mgof_cr() else if (s[i]=="mlnp") s1[i] = &_mm_mgof_mlnp() else if (s[i]=="ksmirnov") s1[i] = &_mm_mgof_ksmirnov() else _error(3498,s[i]+ " invalid") }
// approximate chi2-test if (m=="approx") return(_mm_mgof_approx(s1, f, *h, lambda, arg1))
// Monte Carlo test if (m=="mc") return(_mm_mgof_mc(s1, f, *h, lambda, arg1, dots))
// exhaustive enumeration / random composition test if (s=="mlnp") { if (m=="rc") return(_mm_mgof_rc(NULL, f, *h)) return(_mm_mgof_ee(NULL, f, *h)) } p = mm_which(s:=="mlnp") if (rows(p)>0) { p = p[1] p = p \ (p>1 ? (1::p-1) : J(0,1,.)) \ (p<rows(s) ? (p+1::rows(s)) : J(0 > ,1,.)) s1 = s1[p[| 2 \ .|]] p = invorder(p) } else p = (2::length(s)+1) if (m=="rc") return(_mm_mgof_rc(s1, f, *h, lambda, arg1)[p,]) return(_mm_mgof_ee(s1, f, *h, lambda, anyof(s,"ksmirnov"), dots, arg2)[p,]) }
// large sample chi2-approximation test for multinomial distributions real matrix _mm_mgof_approx( pointer(real scalar function) vector s, // statistics to be tested real colvector f, // observed count real colvector h0, // expected count | real scalar lambda, // lambda for Cressie-Read statistic real scalar nfit) // n. of fitted parameters (default 0) { real scalar n, i, k, nstat real colvector stat pointer(real colvector) scalar h
// compute hypothetical distribution if missing n = colsum(f) k = rows(f) if (rows(h0)<=1) h = &(n :* J(k,1,1/k)) else h = &h0 // compute statistics and p-values nstat = length(s) stat = J(nstat,1,.) for (i=1; i<=nstat; i++) stat[i] = (*s[i])(f,*h,lambda) return(stat, chi2tail(rows(f) - (nfit<. ? nfit : 0) - 1,stat)) }
// Monte Carlo exact test for discrete distributions // (based on sampling from the hypothetical distribution) real matrix _mm_mgof_mc( pointer(real scalar function) vector s, // statistics to be tested real colvector f, // observed count real colvector h0, // expected count | real scalar lambda, // lambda for Cressie-Read statistic real scalar reps0, // replications (default=10000) real scalar dots) // display dots { real scalar n, j, i, k, nstat, reps, uni real scalar doti, ndots real colvector fj, H, stat, p, hp pointer(real colvector) scalar h
// set defaults and compute hypothetical distribution if (args()<5 | reps0>=.) reps = 10000 else reps = reps0 if (args()<6) dots = 0 n = round(colsum(f)) // so that, e.g., 10 are sampled if n=9.9999... k = rows(f) if (uni = (rows(h0)==1)) { uni = 1 h = &(n :* J(k,1,1/k)) H = rangen(1/k,1,k) hp = J(k,1,1/k) } else { uni = mm_isconstant(h0) H = mm_colrunsum(h0) H = H / H[rows(H)] hp = (h0) / colsum(h0) h = &h0 } // compute statistics and p-values nstat = length(s) stat = p = J(nstat,1,0) for (i=1; i<=nstat; i++) stat[i] = (*s[i])(f,*h, lambda, n, hp, H) if (reps<=0) return(stat, J(nstat,1,.)) if (dots) { printf("\n{txt}Percent completed ({res}%g{txt} replications)\n",reps) display("{txt}0 {hline 5} 20 {hline 6} 40 {hline 6} 60 {hline 6} 80 {hl > ine 5} 100") ndots = (reps>=50 ? 1 : (reps>=25 ? 2 : (reps>=10 ? 5 : (reps>=5 ? 10 : (reps>=2 ? 25 : 50))))) doti = 1 } for (j=1; j<=reps; j++) { if (uni) fj = _mm_mgof_mc_usmpl(n, k) // = mm_srswr(n,rows(f),1) else fj = _mm_mgof_mc_smpl(n, H) // = mm_upswr(n, h, 1) (slower) for (i=1; i<=nstat; i++) p[i] = p[i] + (stat[i] <= (*s[i])(fj,*h, lambda, n, hp, H)) if (dots) { if (j*50 >= doti*ndots*reps) { printf(ndots*".") displayflush() doti++ } } } if (dots) display("") return(stat, p/reps) } // Monte Carlo subroutine: sample from uniform multinomial real colvector _mm_mgof_mc_usmpl( real scalar n, // sample size real scalar k) // number of categories { real scalar i real colvector u, res
u = ceil(uniform(n,1)*k) res = J(k,1,0) if (n>1.5^(k-2)) { // faster for small k relative to n for (i=1; i<=k; i++) res[i] = colsum(u:==i) } else { for (i=1;i<=n;i++) res[u[i]] = res[u[i]] + 1 } return(res) } // Monte Carlo subroutine: sample from non-uniform multinomial real colvector _mm_mgof_mc_smpl( real scalar n, // sample size real colvector F) // distribution function (cumulative) { real scalar i, j real colvector u, res
if (n>2.3^(rows(F)-3)) { // faster for small k relative to n u = uniform(n,1) res = J(rows(F),1,0) res[1] = sum(u:<=F[1]) for (j=2;j<=rows(F);j++) { res[j] = sum(u:>F[j-1] :& u:<=F[j]) } return(res) } u = sort(uniform(n,1),1) j=1 res = J(rows(F),1,0) for (i=1;i<=n;i++) { while (u[i]>F[j]) j++ res[j] = res[j]+1 } return(res) }
// Random composition exact test for discrete distributions // -> undocumented; do not use this // the method is highly unreliable if the number of potential // compositions is large real matrix _mm_mgof_rc( pointer(real scalar function) vector s0, // statistics to be tested real colvector f, // observed count real colvector h0, // expected count | real scalar lambda, // lambda for Cressie-Read statistic real scalar reps0) // replications (default=10000) { real scalar n, i, j, k, nstat, reps, lnscale real colvector fj, H, stat, statj, p, jj, se, hp pointer(real colvector) scalar h pointer(real scalar function) vector s
// set defaults and compute hypothetical distribution if (args()<5 | reps0>=.) reps = 10000 else reps = reps0 if (s0==NULL) s = &_mm_mgof_mlnp() else { if (rows(s0)!=1) s = &_mm_mgof_mlnp() \ s0 else s = &_mm_mgof_mlnp() , s0 } n = colsum(f) k = rows(f) if (rows(h0)==1) { h = &(n :* J(k,1,1/k)) H = rangen(1/k,1,k) hp = J(k,1,1/k) } else { h = &h0 H = mm_colrunsum(*h) H = H / H[rows(H)] hp = (*h) / colsum(*h) } // compute statistics and p-values nstat = length(s) stat = statj = p = jj = se = J(nstat,1,0) for (i=1; i<=nstat; i++) stat[i] = (*s[i])(f,*h, lambda, n, hp, H) if (reps<=0) return(stat, J(nstat,2,.)) for (j=1; j<=reps; j++) { fj = mm_rcomposition(n, k) for (i=1; i<=nstat; i++) { statj[i] = (*s[i])(fj,*h, lambda, n, hp, H) if (stat[i] <= statj[i]) { jj[i] = jj[i] + 1 p[i] = p[i] + (exp(-statj[1])-p[i])/jj[i] // mean updating se[i] = se[i] + (exp(-statj[1])^2-se[i])/jj[i] // mean updating } } } p = p:*jj; se = se:*jj // lnscale = log. of inverse sampling probability lnscale = lnfactorial(n+k-1) - lnfactorial(k-1) - lnfactorial(n) - ln(reps) se = exp( lnscale :+ ln(sqrt(reps*(1/(reps-1) * se - (p/reps):^2))) ) // se = mm_ncompositions(n, k) * sqrt((1/(reps-1) * se - (p/reps):^2)/reps) p = exp( lnscale :+ ln(p) ) // p = mm_ncompositions(n, k) / reps * p return(stat, p, se) }
// exhaustive enumeration exact test (performs the exact // multinomial g.o.f. test plus, optionally, exact tests // for specified additional statistics) // WARNING: use this test for very small samples only since // computation time is linear to the number of compositions // =comb(n+k-1,k-1), where n is the sample size and k is // the number of categories. // IMPORTANT EXCEPTION: when the null distribution is uniform, // compution time is determined by the number of partitions // which is magnitudes smaller than the number of compositions real matrix _mm_mgof_ee( pointer(real scalar function) vector s0, // statistics to be tested real colvector f, // observed count real colvector h0, // expected count | real scalar lambda, // lambda for Cressie-Read statistic real scalar force, // do not use partitions (for ksmirnov) real scalar dots, // display dots real scalar counter) // will be replaced { real scalar n, i, k, nstat, uni, w real scalar doti, ndots, reps real colvector fj, j, H, stat, statj, p, hp pointer(real colvector) scalar h pointer(real scalar function) vector s struct mm_subsetinfo scalar info
// set defaults and compute hypothetical distribution if (args()<6) dots = 0 if (s0==NULL) s = &_mm_mgof_mlnp() else { if (rows(s0)!=1) s = &_mm_mgof_mlnp() \ s0 else s = &_mm_mgof_mlnp() , s0 } n = colsum(f) k = rows(f) uni = (force==1 ? 0 : 1) if (rows(h0)==1) { h = &(n :* J(k,1,1/k)) H = rangen(1/k,1,k) hp = J(k,1,1/k) } else { h = &h0 if (mm_isconstant(*h)==0) uni = 0 // reset to 0 if h not constant H = mm_colrunsum(*h) H = H / H[rows(H)] hp = (*h) / colsum(*h) } // compute statistics and p-values nstat = length(s) stat = statj = p = j = J(nstat,1,0) for (i=1; i<=nstat; i++) stat[i] = (*s[i])(f,*h, lambda, n, hp, H) if (uni) { // uniform null distribution: work with partitions if (dots) { reps = mm_npartitions(n,k) printf("\n{txt}Percent completed ({res}%g{txt} partitions)\n",reps) display("{txt}0 {hline 5} 20 {hline 6} 40 {hline 6} 60 {hline 6} 80 > {hline 5} 100") ndots = (reps>=50 ? 1 : (reps>=25 ? 2 : (reps>=10 ? 5 : (reps>=5 ? 10 : (reps>=2 ? 25 : 50))))) doti = 1 } info = _mm_partitionsetup(n, k) while ((fj = _mm_partition(info,1)) != J(0,1,.)) { w = exp(lnfactorial(info.k) - quadcolsum(lnfactorial(_mm_panels(fj) > ))) for (i=1; i<=nstat; i++) { statj[i] = (*s[i])(fj,*h, lambda, n, hp, H) if (stat[i] <= statj[i]) { j[i] = j[i] + w // statj[1] = -ln(p) => p = exp(-statj[1]) p[i] = p[i] + (exp(-statj[1])-p[i])*w/j[i] // mean updatin > g } } if (dots) { if (info.counter*50 >= doti*ndots*reps) { printf(ndots*".") displayflush() doti++ } } } if (dots) display("") counter = info.counter return(stat, p:*j) } if (dots) { reps = mm_ncompositions(n,k) printf("\n{txt}Percent completed ({res}%g{txt} compositions)\n",reps) display("{txt}0 {hline 5} 20 {hline 6} 40 {hline 6} 60 {hline 6} 80 {hl > ine 5} 100") ndots = (reps>=50 ? 1 : (reps>=25 ? 2 : (reps>=10 ? 5 : (reps>=5 ? 10 : (reps>=2 ? 25 : 50))))) doti = 1 } info = _mm_compositionsetup(n, k) while ((fj = _mm_composition(info)) != J(0,1,.)) { for (i=1; i<=nstat; i++) { statj[i] = (*s[i])(fj,*h, lambda, n, hp, H) if (stat[i] <= statj[i]) { j[i] = j[i] + 1 // statj[1] = -ln(p) => p = exp(-statj[1]) p[i] = p[i] + (exp(-statj[1])-p[i])/j[i] // mean updating } } if (dots) { if (info.counter*50 >= doti*ndots*reps) { printf(ndots*".") displayflush() doti++ } } } if (dots) display("") counter = info.counter return(stat, p:*j) }
// Pearson's chi2 statistic subroutine real scalar _mm_mgof_x2( real colvector f, // observed count real colvector h, // expected count | transmorphic matrix opt1, // not used transmorphic matrix opt2, // not used transmorphic matrix opt3, // not used transmorphic matrix opt4) // not used { return(quadcolsum((f-h):^2:/h)) }
// log likelihood ratio statistic subroutine real scalar _mm_mgof_lr( real colvector f, // observed count real colvector h, // expected count | transmorphic matrix opt1, // not used transmorphic matrix opt2, // not used transmorphic matrix opt3, // not used transmorphic matrix opt4) // not used { return(2*quadcolsum(f:*(ln(f)-ln(h)))) }
// Cressie-Read statistic real scalar _mm_mgof_cr( real colvector f, // observed count real colvector h, // expected count | real scalar lambda, // lambda parameter; default is 2/3 transmorphic matrix opt2, // not used transmorphic matrix opt3, // not used transmorphic matrix opt4) // not used { real scalar L L = (lambda<. ? lambda : 2/3) if (abs(L)<1e-6) return( 2 * quadcolsum(f:*(ln(f)-ln(h))) ) else if (abs(L+1)<1e-6) return( 2 * quadcolsum(h:*(ln(h)-ln(f))) ) else if (L>0) return( 2/(L*(L+1)) * quadcolsum(f:*((f:/h):^L :- 1)) ) return( 2/(L*(L+1)) * quadcolsum(f:*((h:/f):^-L :- 1)) ) }
// -ln(p) statistic subroutine (multinomial test) real scalar _mm_mgof_mlnp( real colvector f, // observed count real colvector h, // expected count | transmorphic matrix opt1, // not used real scalar n, // n. of obs [ =colsum(f) ] real colvector p, // hypothical propotions [ = h/colsum(h) ] transmorphic matrix opt4) // not used { if (args()>=5) return(-(lnfactorial(n) - quadcolsum(lnfactorial(f)) + quadcolsum(f:*ln(p)))) return(-(lnfactorial(quadcolsum(f)) - quadcolsum(lnfactorial(f)) + quadcolsum(f:*ln(h/quadcolsum(h))))) }
// kolmogorov-smirnov statistic subroutine real scalar _mm_mgof_ksmirnov( real colvector f, // observed count real colvector h, // expected count | transmorphic matrix opt1, // not used real scalar n, // n. of obs [ =colsum(f) ] transmorphic matrix opt3, // not used real colvector H) // hypothetical cumulative [ = runsum(h)/sum(h) > ] { if (args()>=6) return(max(abs(H-mm_colrunsum(f)/n))) return(max(abs(mm_colrunsum(h)/colsum(h)-mm_colrunsum(f)/colsum(f)))) }
end
*! *! mm_root.mata *! version 1.0.0, Ben Jann, 07jul2006 *! translation of zeroin.c from http://www.netlib.org/c/brent.shar
* Quoted from http://www.netlib.org/c/brent.shar: * X * Algorithm * X * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical * X * computations. M., Mir, 1980, p.180 of the Russian edition * X * * X * The function makes use of the bissection procedure combined with * X * the linear or quadric inverse interpolation. * X * At every step program operates on three abscissae - a, b, and c. * X * b - the last and the best approximation to the root * X * a - the last but one approximation * X * c - the last but one or even earlier approximation than a that * X * 1) |f(b)| <= |f(c)| * X * 2) f(b) and f(c) have opposite signs, i.e. b and c confine * X * the root * X * At every step Zeroin selects one of the two new approximations, the * X * former being obtained by the bissection procedure and the latter * X * resulting in the interpolation (if a,b, and c are all different * X * the quadric interpolation is utilized, otherwise the linear one). * X * If the latter (i.e. obtained by the interpolation) point is * X * reasonable (i.e. lies within the current interval [b,c] not being * X * too close to the boundaries) it is accepted. The bissection result * X * is used in the other case. Therefore, the range of uncertainty is * X * ensured to be reduced at least by the factor 1.6
* Main changes compared to the original: * 1. mm_root() returns an error code and stored the solution in x. * 2. The maximum # of iteration may be set. * 3. The tolerance argument is optional. * 4. Additional arguments may be passed on to f. * 5. Special action is taken if f returns missing, convergence is not * reached, or f(ax) and f(bx) do not have opposite signs.
version 9.2
local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"
mata:
real scalar mm_root( transmorphic x, // bj: will be replaced by solution pointer(real scalar function) scalar f, // Address of the function whose zero will be sought for real scalar ax, // Root will be sought for within a range [ax,bx] real scalar bx, // | real scalar tol, // Acceptable tolerance for the root value (default 0) real scalar maxit, // bj: maximum # of iterations (default: 1000) `opts') // bj: additional args to pass on to f { transmorphic fs // setup for f real scalar a, b, c // Abscissae, descr. see above real scalar fa, fb, fc // f(a), f(b), f(c) real scalar prev_step // Distance from the last but one real scalar tol_act // Actual tolerance real scalar p // Interpolation step is calcu- real scalar q // lated in the form p/q; divi- // sion operations is delayed // until the last moment real scalar new_step // Step at this iteration real scalar t1, cb, t2 real scalar itr
if (args()<5) tol = 0 // bj: set tolerance if (args()<6) maxit = 1000 // bj: maximum # of iterations
fs = mm_callf_setup(f, args()-6, `opts') // bj: prepare function call
x = . // bj: initialize output
a = ax; b = bx; fa = mm_callf(fs, a); fb = mm_callf(fs, b) c = a; fc = fa
if ( fa==. ) return(0) // bj: abort if fa missing
if ( (fa > 0 & fb > 0) | // bj: f(ax) and f(bx) do not (fa < 0 & fb < 0) ) { // have opposite signs if ( abs(fa) < abs(fb) ) { x = a; return(2) // bj: fa closer to zero than fb } x = b; return(3) // bj: fb closer to zero than fa }
for (itr=1; itr<=maxit; itr++) { if ( fb==. ) return(0) // bj: abort if fb missing
prev_step = b-a
if( abs(fc) < abs(fb) ) { // Swap data for b to be the a = b; b = c; c = a; // best approximation fa = fb; fb = fc; fc = fa }
tol_act = 2*epsilon(b) + tol/2 new_step = (c-b)/2
if( abs(new_step) <= tol_act | fb == 0 ) { x = b // Acceptable approx. is found return(0) }
// Decide if the interpolation can be tried if( abs(prev_step) >= tol_act // If prev_step was large enough & abs(fa) > abs(fb) ) { // and was in true direction, // Interpolation may be tried cb = c-b if( a==c ) { // If we have only two distinct t1 = fb/fa // points linear interpolation p = cb*t1 // can only be applied q = 1.0 - t1 } else { // Quadric inverse interpolation q = fa/fc; t1 = fb/fc; t2 = fb/fa p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) ) q = (q-1.0) * (t1-1.0) * (t2-1.0) } if( p>0 ) // p was calculated with the op- q = -q // posite sign; make p positive else // and assign possible minus to p = -p // q if( p < (0.75*cb*q-abs(tol_act*q)/2) // If b+p/q falls in [b,c] & p < abs(prev_step*q/2) ) // and isn't too large new_step = p/q // it is accepted // If p/q is too large then the // bisection procedure can // reduce [b,c] range to more // extent }
if( abs(new_step) < tol_act ) { // Adjust the step to be not less if( new_step > 0 ) // than tolerance new_step = tol_act else new_step = -tol_act }
a = b; fa = fb // Save the previous approx. b = b + new_step; fb = mm_callf(fs, b) // Do step to a new approxim. if( (fb > 0 & fc > 0) | (fb < 0 & fc < 0) ) { // Adjust c for it to have a sign opposite to that of b c = a; fc = fa } } x = b return(1) // bj: convergence not reached }
end
*! *! mm_nrroot.mata *! version 1.0.0, Ben Jann, 07jul2006
version 9.2
local opts "o1, o2, o3, o4, o5, o6, o7, o8, o9, o10"
mata:
real scalar mm_nrroot( real scalar x, // initial guess; will be replaced by solution pointer(function) scalar f, // the function; must return f() and f'() | real scalar tol, // acceptable tolerance (default 0) real scalar maxit, // maximum # of iterations (default: 1000) `opts') // additional args to pass on to f { transmorphic fs // setup for f real scalar itr real scalar fx // fx = ( f(x), f'(x) ) real scalar dx // dx = fx[1] / fx[2] real scalar tol_act // actual tolerance
if (args()<3) tol = 0 // default tolerance if (args()<4) maxit = 1000 // default maximum # of iterations
fs = mm_callf_setup(f, args()-4, `opts') // prepare function call for (itr=1; itr<=maxit; itr++) { tol_act = 2e+4*epsilon(x) + tol/2 fx = mm_callf(fs, x) dx = fx[1] / fx[2] x = x - dx if ( abs(dx) <= tol_act | missing(x) ) return(0) } return(1) // convergence not reached }
end
*! *! mm_finvert.mata *! version 1.0.0, Ben Jann, 07jul2006
version 9.2
mata:
real vector mm_finvert( y, f, // nr: brent: o1, // fd lo | o2, // x0 up tol, maxit) { if (args()<5) tol = 0 if (args()<6) maxit = 1000 if (ispointer(o1)) { if (args()<4) o2 = y // set initial guess to y return(_mm_finvert_nr(y, f, o1, o2, tol, maxit)) } if (args()==3) _error(3001, "expected 4 to 6 arguments but received 3") return(_mm_finvert_brent(y, f, o1, o2, tol, maxit)) }
real vector _mm_finvert_nr( real vector y, pointer(real scalar function) scalar f, pointer(real scalar function) scalar df, real vector x, real scalar tol, real scalar maxit) { real scalar i, I, resi, rc real vector res pointer scalar ix
I = length(y) ix = (length(x)==1 ? &1 : (length(x)<I ? _error(3200) : &i))
res = J(rows(y), cols(y), .) for (i=1; i<=I; i++) { rc = mm_nrroot(resi=x[*ix], &_mm_finvert_nr_fdf(), tol, maxit, y[i], f, df) if (rc) _error(3360, "failure to converge for element " + strofreal(i)) res[i] = resi } return(res) }
real rowvector _mm_finvert_nr_fdf( real scalar x, real scalar y, pointer(real scalar function) scalar f, pointer(real scalar function) scalar df) { return( (*f)(x) - y, (*df)(x) ) }
real vector _mm_finvert_brent( real vector y, pointer(real scalar function) scalar f, real vector lo, real vector up, real scalar tol, real scalar maxit) { real scalar i, I, resi, rc real vector res pointer scalar il, iu
I = length(y) il = (length(lo)==1 ? &1 : (length(lo)<I ? _error(3200) : &i)) iu = (length(up)==1 ? &1 : (length(up)<I ? _error(3200) : &i))
res = J(rows(y), cols(y), .) for (i=1; i<=I; i++) { rc = mm_root(resi=., &_mm_finvert_brent_f(), lo[*il], up[*iu], tol, maxit, y[i], f) if (rc==1) _error(3360, "failure to converge for element " + strofreal(i)) if (rc) _error(3498, "invalid search range for element " + strofreal(i)) res[i] = resi } return(res)
}
real scalar _mm_finvert_brent_f( real scalar x, real scalar y, pointer(real scalar function) scalar f) { return( (*f)(x) - y ) }
end
*! *! mm_plot.mata *! version 1.0.0, Ben Jann, 03aug2006 version 9.2 mata:
void mm_plot( real matrix X, | string scalar type, string scalar opts) { real scalar rc
rc = _mm_plot(X, type, opts) if (rc) _error(rc) }
real scalar _mm_plot( real matrix X, | string scalar type, string scalar opts) { real scalar n, N, k, j, rc, updata string scalar vars, plottype real vector vid
updata = st_updata() plottype = type if (plottype=="") plottype = "scatter" n = rows(X) k = cols(X) N = st_nobs() if (N<n & k>0) st_addobs(n-N) vid = J(1, k, .) for (j=1; j<=k; j++) { st_store((1,n), vid[j]=st_addvar("double", st_tempname()), X[,j]) st_varlabel(vid[j], "Column " + strofreal(j)) vars = vars + " " + st_varname(vid[j]) } st_updata(updata) rc = _stata("twoway " + plottype + vars + ", " + opts) updata = st_updata() // twoway might change data-have-changed flag if (N<n & k>0) st_dropobsin((N+1,n)) st_dropvar(vid) st_updata(updata) return(rc) }
end
*! *! mm_panels.mata *! version 1.0.3, Ben Jann, 03may2007 version 9.2 mata:
void function mm_panels( transmorphic vector X1, // upper level ID variable (e.g. strata) transmorphic matrix info1, // will be replaced (unless X1 and X2 missing) | transmorphic vector X2, // lower level ID variable (e.g. clusters) transmorphic matrix info2) // will be replaced (unless X2 missing) { real scalar i, j, nc, b, e, b2, e2
if (length(X1)>0 & X1!=.) info1 = _mm_panels(X1) if (args()<3) return if (length(X2)==0 | X2==.) { if (length(X1)>0 & X1!=.) info1 = info1, info1 return } if (length(X1)==0 | X1==.) { info2 = _mm_panels(X2) info1 = colsum(info2), rows(info2) return } if (length(X1)!=length(X2)) _error(3200) info1 = info1, J(rows(info1),1,.) e = 0 for (i=1; i<=rows(info1); i++) { nc = 1 b = e + 2 e = e + info1[i,1] for (j=b; j<=e; j++) { if (X2[j]!=X2[j-1]) nc++ } info1[i,2] = nc } if (args()<4) return info2 = J(colsum(info1[.,2]), 1, .) e = e2 = 0 for (i=1; i<=rows(info1); i++) { b = e + 1 e = e + info1[i,1] b2 = e2 + 1 e2 = e2 + info1[i,2] info2[|b2 \ e2|] = _mm_panels(X2[|b \ e|], info1[i,2]) } }
real colvector _mm_panels(transmorphic vector X, | real scalar np) { real scalar i, j, r real colvector res
r = length(X) if (r<1) return(J(0,1,.)) if (args()<2) np = r res = J(np, 1, 1) j = 1 for (i=2; i<=r; i++) { if (X[i]==X[i-1]) res[j] = res[j] + 1 else j++ } if (j==r) return(res) return(res[|1 \ j|]) }
// /*Old (slow) version of _mm_panels()*/ //real colvector _mm_panels(transmorphic vector X, | real scalar np) //{ // real scalar i, j, n // real colvector res // // if (args()<2) np = mm_npanels(X) // if (length(X)==0) return(J(0,1,.)) // res = J(np, 1, .) // n = j = 1 // for (i=2; i<=length(X); i++) { // if (X[i]!=X[i-1]) { // res[j++] = n // n = 1 // } // else n++ // } // res[j] = n // return(res) //}
real scalar mm_npanels(vector X) { real scalar i, np
if (length(X)==0) return(0) np = 1 for (i=2; i<=length(X); i++) { if (X[i]!=X[i-1]) np++ } return(np) }
end
*! *! mm_nunique.mata *! version 1.0.2 16may2007 version 9.2 mata:
real scalar mm_nunique(transmorphic vector x) { if (cols(x)==1) return(mm_npanels(sort(x,1))) else { if (iscomplex(x)) return(mm_npanels(sort(transposeonly(x),1))) else return(mm_npanels(sort(x',1))) } }
real scalar mm_nuniqrows(transmorphic matrix x) { real scalar i, n, ns real colvector p
if (rows(x)==0) return(0) if (cols(x)==0) return(1)
p = order(x, 1..cols(x)) ns = 1 n = rows(x) for (i=2;i<=n;i++) { if (x[p[i-1],]!=x[p[i],]) ns++ } return(ns) }
end
*! *! mm_isconstant.mata *! version 1.0.0, Ben Jann, 01mar2007 version 9.0 mata:
real scalar mm_isconstant(X) { if (length(X)<2) return(1) return(all(X:==X[1,1])) }
end
*! *! mm_nobs.mata *! version 1.0.0, Ben Jann, 07jun2006 version 9.2 mata:
real scalar mm_nobs(transmorphic matrix x, real colvector w) { return(rows(w)==1 ? rows(x)*w : (rows(x)==rows(w) ? colsum(w) : _error(3200))) } end
*! *! mm_colrunsum.mata *! version 1.0.8 11jan2008 Ben Jann version 9.0 mata: numeric matrix mm_colrunsum(numeric matrix A) { numeric matrix B
if (stataversion()>=1000) return(_mm_colrunsum_goto10(A))
if (isfleeting(A)) { _mm_colrunsum(A) return(A) } _mm_colrunsum(B=A) return(B) }
void _mm_colrunsum(numeric matrix Z) { real scalar i
_editmissing(Z, 0) for (i=2; i<=rows(Z); i++) Z[i,] = Z[i-1,] + Z[i,]
// // running sum using mean-update formula; see Gould 2006, SJ 6(4) // _editmissing(Z, 0) // if (rows(Z)<2) return // for (i=2; i<=rows(Z); i++) Z[i,] = Z[i-1,] + (Z[i,]-Z[i-1,])/i // Z = Z :* (1::rows(Z))
}
numeric matrix _mm_colrunsum_goto10(numeric matrix Z) { return(_mm_colrunsum10(Z)) }
end
*! *! mm_linbin.mata *! version 1.0.8 03aug2007 Ben Jann version 9.0 mata:
real colvector mm_linbin( real colvector x, // data points real colvector w, // weights real colvector g) // grid points { real colvector c, p, ng, nx real scalar i, j, d, ww
ng = rows(g) nx = rows(x) ww = rows(w)!=1 if ((ww & rows(w)!=nx) | ng<1) _error(3200) if (nx<1) return(J(ng, 1, 0)) if (ng<2) return(ww ? colsum(w) : w*nx) p = order((x,(1::nx)),(1,2)) //stable sort order c = J(ng, 1, 0) for (i=1; i<=nx; i++) { // data below grid if (x[p[i]]>=g[1]) break c[1] = c[1] + (ww ? w[p[i]] : w) } j = 2 for (; i<=nx; i++) { if (g[ng]<x[p[i]]) break while (g[j]<x[p[i]]) j++ d = (ww ? w[p[i]] : w) / (g[j] - g[j-1]) c[j-1] = c[j-1] + (g[j] - x[p[i]]) * d c[j] = c[j] + (x[p[i]] - g[j-1]) * d } if (i<=nx) { // data above grid c[ng] = c[ng] + (ww ? colsum(w[p[|i \ .|]]) : w*(nx-i+1)) } return(c) }
end
*! *! mm_fastlinbin.mata *! version 1.0.2 03aug2007 Ben Jann version 9.1 mata:
real colvector mm_fastlinbin( real colvector x, // data points real colvector w, // weights real colvector g) // grid points { real colvector c real scalar g1, i, j, d, z, ng, nx pointer scalar wi
ng = rows(g) nx = rows(x) if (ng<1) _error(3200) if (nx<1) return(J(ng, 1, 0)) if (ng<2) return(mm_nobs(x,w)) wi = (rows(w)==1 ? &1 : (rows(w)!=nx ? _error(3200) : &i)) g1 = g[1] d = (g[ng] - g1)/(ng-1) c = J(ng, 1, 0) for (i=1; i<=rows(x); i++) { z = (x[i] - g1) / d if (z<=0) { c[1] = c[1] + w[*wi] continue } if (z>=(ng-1)) { c[ng] = c[ng] + w[*wi] continue } j = trunc(z) + 1 c[j] = c[j] + w[*wi] * (j-z) j++ c[j] = c[j] + w[*wi] * (2-j+z) } return(c) } end
*! *! mm_exactbin.mata *! version 1.0.5 03aug2007 Ben Jann version 9.0 mata:
real colvector mm_exactbin( real colvector x, // data real colvector w, // weights real colvector g, // grid | real scalar dir, // 0 left inclusive (default), else right inclusive real scalar include) // include data outside grid in first/last bin { real colvector c, p, nx, ng real scalar i, j, ww
if (args()<4) dir = 0 if (args()<5) include = 0 nx = rows(x) ng = rows(g) ww = rows(w)!=1 if ((ww & rows(w)!=nx) | rows(g)<2) _error(3200) if (nx<1) return(J(ng-1, 1, 0)) p = order((x,(1::nx)),(1,2)) //stable sort order if (include==0) { if (x[p[1]]<g[1] | x[p[nx]]>g[ng]) _error("data out of grid range") } if (ng<3) return(ww ? colsum(w) : w*nx) c = J(ng-1, 1, 0) if (dir==0) { j = rows(c) for (i=nx; i>=1; i--) { while (g[j]>x[p[i]]) { if (j<=1) break j-- } c[j] = c[j] + (ww ? w[p[i]] : w) } } else { j = 2 for (i=1; i<=nx; i++) { while (g[j]<x[p[i]]) { if (j>=ng) break j++ } c[j-1] = c[j-1] + (ww ? w[p[i]] : w) } } return(c) } end
*! *! mm_makegrid.mata *! version 1.0.4, Ben Jann, 26jun2006 version 9.0 mata:
real colvector mm_makegrid( real colvector x, // data points | real scalar M, // number of bins (default: 401) real scalar e, // extend range by e on each side real scalar min, // min real scalar max) // max { real scalar a, b
a = (min<. ? min : min(x) - (e<. ? e : 0)) b = (max<. ? max : max(x) + (e<. ? e : 0)) return( rangen(a, b, (M<. ? M : 512)) ) }
end
*! *! mm_cut.mata *! version 1.0.0, Ben Jann, 16mar2006 version 9.1 mata:
real colvector mm_cut(real colvector x, real vector at, | real scalar sorted) { real scalar i, j real colvector result, p
result = J(rows(x),1,.) j = length(at) if (args()<3) sorted = 0 if (sorted) { for (i=rows(x); i>0; i--) { if (x[i]>=.) continue for (; j>0; j--) { if (at[j]<=x[i]) break } if (j>0) result[i,1] = at[j] } return(result) } p = order(x,1) for (i=rows(x); i>0; i--) { if (x[p[i]]>=.) continue for (; j>0; j--) { if (at[j]<=x[p[i]]) break } if (j>0) result[p[i],1] = at[j] } return(result) } end
*! *! mm_posof.mata *! version 1.0.0, Ben Jann, 24mar2006 version 9.0 mata:
real scalar mm_posof(transmorphic vector haystack, transmorphic scalar needle) { real scalar i
for (i=1; i<=length(haystack); i++) { if (haystack[i]==needle) return(i) } return(0) } end
*! *! mm_which.mata *! version 1.0.2, Ben Jann, 17apr2007 version 9.2 mata:
real matrix mm_which(real vector I) { if (cols(I)!=1) return(select(1..cols(I), I)) else return(select(1::rows(I), I)) }
end
*! *! mm_locate.mata *! version 1.0.0, Ben Jann, 07jul2006
version 9.2
mata:
void mm_locate( // translation of -locate- from Press et al. (1992), Numerical // Recipes in C, p. 117, http://www.numerical-recipes.com // Description: "Given an array xx[1..n], and given a value x, // returns a value j such that x is between xx[j] and xx[j+1]. // xx must be monotonic, either increasing or decreasing. j=0 // or j=n is returned to indicate that x is out of range." numeric vector xx, numeric scalar x, real scalar j) { real scalar ju, jm, jl real scalar ascnd real scalar n
n = length(xx) jl = 0 ju = n + 1 ascnd = (xx[n] >= xx[1]) while (ju-jl > 1) { jm = trunc((ju+jl)/2) if (x >= xx[jm] == ascnd) jl = jm else ju = jm } if (x == xx[1]) j = 1 else if(x == xx[n]) j = n - 1 else j = jl }
void mm_hunt( // translation of -hunt- from Press et al. (1992), Numerical // Recipes in C, p. 118-9, http://www.numerical-recipes.com // Description: "Given an array xx[1..n], and given a value x, // returns a value jlo such that x is between xx[jlo] and // xx[jlo+1]. xx[1..n] must be monotonic, either increasing or // decreasing. jlo=0 or jlo=n is returned to indicate that x is // out of range. jlo on input is taken as the initial guess for // jlo on output. numeric vector xx, numeric scalar x, real scalar jlo) { real scalar jm, jhi, inc real scalar ascnd real scalar n
n = length(xx) ascnd=(xx[n] >= xx[1]) if (jlo <= 0 | jlo > n) { jlo = 0 jhi = n + 1 } else { inc = 1 if (x >= xx[jlo] == ascnd) { if (jlo == n) return jhi = jlo + 1 while (x >= xx[jhi] == ascnd) { jlo = jhi inc = inc + inc jhi = jlo + inc if (jhi > n) { jhi = n + 1 break } } } else { if (jlo == 1) { jlo = 0 return } jhi = jlo-- while (x < xx[jlo] == ascnd) { jhi = jlo inc = inc + inc if (inc >= jhi) { jlo = 0 break } else jlo = jhi - inc } } } while (jhi-jlo != 1) { jm = trunc((jhi+jlo)/2) if (x >= xx[jm] == ascnd) jlo = jm else jhi = jm } if (x == xx[n]) jlo = n - 1 if (x == xx[1]) jlo = 1 }
end
*! *! mm_ipolate.mata *! version 1.0.6, Ben Jann, 10jul2006 version 9.2 mata:
real colvector mm_ipolate(real colvector x, real colvector y, real colvector xnew, | real scalar outer) { real scalar i, j0b, j0e, j1b, j1e, r, xlo, xup, xi, y0, y1 real colvector p, pnew, ynew
r = rows(x) if (rows(y)!=r) _error(3200) if (r<1) return(J(rows(xnew), 1, .)) if (args()<4) outer = 0
p = order(x, 1) pnew = order(xnew, 1) ynew = J(rows(xnew),1,.) xlo = x[p[1]]; xup = x[p[rows(p)]] j0b = j1e = j0e = 1 for (i=1; i<=rows(xnew); i++) { xi = xnew[pnew[i]] if (outer==0) { if (xi<xlo) continue if (xi>xup) return(ynew) } while (j0e<r) { if (x[p[j0e+1]]>xi) break j0e++ if (x[p[j0e]]>x[p[j0b]]) j0b = j0e } if (j0e>=j1e) { j1b = j0e while (j1b<=r) { if (x[p[j1b]]>=xi) break j1b++ } if (j1b>r) j1b = r j1e = j1b while (j1e<r) { if (x[p[j1e+1]]>x[p[j1b]]) break j1e++ } } y0 = (j0b==j0e ? y[p[j0b]] : mean(y[p[|j0b \ j0e|]],1)) y1 = (j1b==j1e ? y[p[j1b]] : mean(y[p[|j1b \ j1e|]],1)) if (outer) { if (xi<xlo) { ynew[pnew[i]] = y1 continue } if (xi>xup) { ynew[pnew[|i \ rows(pnew)|]] = J(rows(pnew)-i+1,1,y0) return(ynew) } } ynew[pnew[i]] = ( j0e==j1e ? y0 : y0 + (y1-y0) * (xi-x[p[j0e]])/(x[p[j1b]]-x[p[j0e]]) ) } return(ynew) } end
*! *! mm_polint.mata *! version 1.0.0, Ben Jann, 12jul2006 version 9.2 mata:
real vector mm_polint( real vector xa, real vector ya, real vector x, | real scalar degree) // degree { real scalar n, m, i, j, a, b real vector y
if ( args()<4 ) degree = 1 n = length(xa) if (length(ya)!=n) _error(3200) if ( degree<1 ) _error(3498, "degree must be 1 or larger") if ( degree>=n ) _error(3498, "degree must be smaller than length of data") if ( degree!=trunc(degree) ) _error(3498, "degree must be integer")
m = degree + 1 y = J(rows(x), cols(x), .) j = 0 for (i=1; i<=length(x); i++) { mm_hunt(xa, x[i], j) a = min((max((trunc(j-(m-1)/2+.5), 1)), n+1-m)) b = a+m-1 y[i] = _mm_polint(xa[|a \ b|], ya[|a \ b|], x[i]) } return(y) }
real scalar _mm_polint( // translation of -polint- from Press et al. (1992), Numerical // Recipes in C, p. 109-110, http://www.numerical-recipes.com real vector xa, real vector ya, real scalar x) { real scalar i, m, n, ns real scalar y, den, dif, dift, ho, hp, w real vector c, d
n = length(xa) if ( length(ya) != n ) _error(3200)
ns = 1 dif = abs(x-xa[1]) c = J(1, n, .) d = J(1, n, .) for (i=1; i<=n; i++) { if ( (dift=abs(x-xa[i])) < dif) { ns = i dif = dift } c[i] = ya[i] d[i] = ya[i] } y = ya[ns--] for (m=1; m<n; m++) { for (i=1; i<=n-m; i++) { ho = xa[i] - x hp = xa[i+m] - x w = c[i+1] - d[i] if ( (den=ho-hp) == 0.0) _error(3498, "invalid input data") // This error can occur // only if two input xa’s are (to within roundoff) identical. den = w / den d[i] = hp * den c[i] = ho * den } if ( 2*ns < (n-m) ) y = y + c[ns+1] else y = y + d[ns--] } return(y) }
end
*! *! mm_cond.mata *! version 1.0.0 29feb2008 Ben Jann version 9.2 mata:
transmorphic matrix mm_cond(real matrix x, transmorphic matrix y, transmorphic > matrix z) { transmorphic matrix res real scalar r, R, c, C transmorphic scalar rx, cx, ry, cy, rz, cz
if (eltype(y) != eltype(z)) _error(3250) R = max((rows(x),rows(y),rows(z))) C = max((cols(x),cols(y),cols(z))) rx = (rows(x)==1 ? &1 : (rows(x)<R ? _error(3200) : &r)) cx = (cols(x)==1 ? &1 : (cols(x)<C ? _error(3200) : &c)) ry = (rows(y)==1 ? &1 : (rows(y)<R ? _error(3200) : &r)) cy = (cols(y)==1 ? &1 : (cols(y)<C ? _error(3200) : &c)) rz = (rows(z)==1 ? &1 : (rows(z)<R ? _error(3200) : &r)) cz = (cols(z)==1 ? &1 : (cols(z)<C ? _error(3200) : &c)) res = J(R,C, missingof(y)) for (r=1;r<=R;r++) { for (c=1;c<=C;c++) { res[r,c] = (x[*rx,*cx] ? y[*ry,*cy] : z[*rz,*cz]) } } return(res) }
end
*! *! mm_expand.mata *! version 1.0.0, Ben Jann, 07jun2006 version 9.2 mata:
transmorphic matrix mm_expand(transmorphic matrix X, real vector wr, | real vector wc, real scalar sort) { transmorphic matrix Xnew
if (args()<3) wc = 1 if (args()<4) sort = 0 if (isfleeting(X)) { _mm_expand(X, wr, wc, sort) return(X) } else { _mm_expand(Xnew=X, wr, wc, sort) return(Xnew) } }
void _mm_expand(transmorphic matrix X, real vector wr, real vector wc, real scalar sort) { real scalar i, b, e, n, add real colvector f pointer scalar fi
// expand rows if (length(wr)==1 & sort==0) _mm_repeat(X, wr, 1) else { f = trunc(wr:-1) :* (wr:>0) _editmissing(f,0) n = rows(X) if (length(f)==1) { fi = &1 add = n*f } else { if (n!=length(f)) _error(3200) fi = &i add = sum(f) } if (add>0) { X = X \ J(add,cols(X),missingof(X)) if (cols(X)>0) { if (sort) f = f :+ 1 b = rows(X)+1 for (i=n; i>=1; i--) { if (f[*fi]<1) continue e = b - 1 b = b - f[*fi] X[|b,1 \ e,.|] = X[J(f[*fi],1,i), .] } } } } // expand columns if (length(wc)==1 & sort==0) { _mm_repeat(X, 1, wc) return } f = trunc(wc:-1) :* (wc:>0) _editmissing(f,0) n = cols(X) if (length(f)==1) { fi = &1 add = n*f } else { if (n!=length(f)) _error(3200) fi = &i add = sum(f) } if (add>0) { X = X , J(rows(X),add,missingof(X)) if (rows(X)>0) { if (sort) f = f :+ 1 b = cols(X)+1 for (i=n; i>=1; i--) { if (f[*fi]<1) continue e = b - 1 b = b - f[*fi] X[|1,b \ .,e|] = X[., J(1,f[*fi],i)] } } } }
transmorphic matrix mm_repeat(transmorphic matrix X, real scalar wr, | real scalar wc) { transmorphic matrix Xnew
if (args()<3) wc = 1 if (isfleeting(X)) { _mm_repeat(X, wr, wc) return(X) } else { _mm_repeat(Xnew=X, wr, wc) return(Xnew) } }
void _mm_repeat(transmorphic matrix X, real scalar wr, real scalar wc) { real scalar i, rr, rc, r, c
r = rows(X) c = cols(X) rr = max((trunc(wr-1),0)) if (rr>0) { X = X \ J(rr*r, c, missingof(X)) if (r>0 & c>0) for (i=1; i<=rr; i++) X[|i*r+1,1 \ i*r+r,.|] = X[|1,1 \ r,.|] } rc = max((trunc(wc-1),0)) if (rc>0) { X = X , J(rows(X), rc*c, missingof(X)) if (r>0 & c>0) for (i=1; i<=rc; i++) X[|1,i*c+1 \ .,i*c+c|] = X[|1,1 \ .,c|] } }
end
*! *! mm_unorder2.mata *! version 1.0.0, Ben Jann, 29mar2006 version 9.0 mata:
real colvector mm_unorder2(real scalar n) return(order(uniform(n,2),(1,2)))
end
*! *! mm_jumble2.mata *! version 1.0.0, Ben Jann, 29mar2006 version 9.0 mata:
transmorphic matrix mm_jumble2(transmorphic matrix x) { return(x[mm_unorder2(rows(x)), .]) }
end
*! *! mm__jumble2.mata *! version 1.0.0, Ben Jann, 29mar2006 version 9.0 mata:
function mm__jumble2(transmorphic matrix x) { _collate(x,mm_unorder2(rows(x))) }
end
*! *! mm_pieces.mata *! version 1.0.1, Ben Jann, 30apr2007 version 9.2 mata:
string rowvector mm_pieces( string scalar s, real scalar w, | real scalar nobreak) { real scalar i, j, k, n, l, b, nobr string scalar c string rowvector res
nobr = ( args()>2 ? nobreak : 0 ) l = strlen(s) if (l<2 | w>=l) return(strtrim(s)) res = J(1, _mm_npieces(s, w, nobr), "") j = k = n = 0 b = 1 for (i=1; i<=l; i++) { c = substr(s, i, 1) if (j<1) { // skip to first nonblank character if (c==" ") { b++ continue } } j++ if (i==l) res[++n] = strtrim(substr(s, b, .)) else { if (c==" ") k = i if (j>=w) { if (k<1) { if (nobr) continue k = i } else { if (substr(s, i+1, 1)==" ") k = i } res[++n] = strtrim(substr(s, b, k-b+1)) j = i - k b = k + 1 k = 0 } } } return(res) }
real matrix mm_npieces( string matrix S, real scalar w, | real scalar nobreak) { real scalar i, j real matrix res
res = J(rows(S),cols(S),1) if (w>=.) return(res) for (i=1; i<=rows(S); i++) { for (j=1; j<=cols(S); j++) { res[i,j] = _mm_npieces(S[i,j], w, (args()>2 ? nobreak : 0)) } } return(res) }
real scalar _mm_npieces( string scalar s, real scalar w, | real scalar nobreak) { real scalar i, j, k, n, l, nobr string scalar c
nobr = ( args()>2 ? nobreak : 0 ) l = strlen(s) if (l<2 | w>=l) return(1) j = k = n = 0 for (i=1; i<=l; i++) { c = substr(s, i, 1) if (j<1) { // skip to first nonblank character if (c==" ") continue } j++ if (i==l) n++ else { if (c==" ") k = i if (j>=w) { if (k<1) { if (nobr) continue k = i } else { if (substr(s, i+1, 1)==" ") k = i } n++ j = i - k k = 0 } } } if (n==0) n++ return(n) } end
*! *! mm_invtokens.mata *! version 1.0.2 Ben Jann 05jan2008 version 9.0 mata:
string scalar mm_invtokens(string vector In, | real scalar noclean) { string scalar Out string scalar tok real scalar i
if (args()<2) noclean = 0 Out = "" for (i=1; i<=length(In); i++) { tok = In[i] if (noclean) tok = "`" + `"""' + tok + `"""' + "'" else if (strpos(tok, `"""')) tok = "`" + `"""' + tok + `"""' + "'" else if (strpos(tok, " ")) tok = `"""' + tok + `"""' else if (tok=="") tok = `"""' + `"""' if (i>1) tok = " " + tok Out = Out + tok } return(Out) } end
*! *! mm_realofstr.mata *! version 1.0.2, Ben Jann, 24mar2006 version 9.0 mata:
real matrix mm_realofstr(string matrix S) { real matrix R string scalar tmp real scalar i, j
R = J(rows(S), cols(S), .) tmp = st_tempname() for (i=1; i<=rows(S); i++) { for (j=1; j<=cols(S); j++) { stata("scalar " + tmp + "=real(`" + `"""' + S[i,j] + `"""' + "')") R[i,j] = st_numscalar(tmp) } } stata("capture scalar drop " + tmp) return(R) } end
*! *! mm_strexpand.mata *! version 1.0.4, Ben Jann, 30apr2007 version 9.0 mata:
string scalar mm_strexpand(string scalar s, string vector slist, | string scalar def, real scalar unique, string scalar errtxt) { real scalar err string scalar res
if (args()<5) errtxt = `"""' + s + `"" invalid"' if (args()<4) unique = 0 err = _mm_strexpand(res, s, slist, def, unique) if (err) _error(err, errtxt) return(res) }
real scalar _mm_strexpand(res, string scalar s, string vector slist, | string scalar def, real scalar unique) { real scalar i, l, match
if (s=="") { res = def return(0) } if (args()<5) unique = 0 l = strlen(s) if (unique) { match = 0 for (i=1; i<=length(slist); i++) { if (s==substr(slist[i], 1, l)) { if (match) return(3498) match = i } } if (match) { res = slist[match] return(0) } } else { for (i=1; i<=length(slist); i++) { if (s==substr(slist[i], 1, l)) { res = slist[i] return(0) } } } return(3499) } end
*! *! mm_matlist.mata *! version 1.0.4 31aug2007 Ben Jann version 9.2 mata:
void mm_matlist( real matrix X, | string matrix fmt, real scalar b, string vector rstripe, string vector cstripe, string scalar rtitle, transmorphic matrix res ) { real scalar i, j, rw, lw, j0, j1, l, r, save real rowvector wd string rowvector sfmt pointer(real scalar) scalar fi, fj pointer(string vector) scalar rstr, cstr
if (args()<2) fmt = "%9.0g" if (args()<3) b = 1 save = (args()==7) //wd = strlen(sprintf(fmt,-1/3))
if ((rows(fmt)!=1 & rows(fmt)!=rows(X)) | (cols(fmt)!=1 & cols(fmt)!=cols(X))) _error(3200) if (!anyof((0,1,2,3),b)) _error(3300)
fi = (rows(fmt)!=1 ? &i : &1) fj = (cols(fmt)!=1 ? &j : &1)
wd = J(1,cols(X),0) for (i=1;i<=rows(X);i++) { for (j=1;j<=cols(X);j++) { wd[j] = max((wd[j],strlen(sprintf(fmt[*fi,*fj],X[i,j])))) } }
if (length(X)==0) return
if (length(rstripe)<1) rstr = &strofreal(1::rows(X)) else rstr = &rstripe if (length(cstripe)<1) cstr = &strofreal(1..cols(X)) else if (rows(cstripe)!=1) cstr = &(cstripe') else cstr = &cstripe
if ((*rstr!="" & length(*rstr)!=rows(X)) | (*cstr!="" & length(*cstr)!=cols(X))) _error(3200)
if (*rstr!="") rw = max((max(strlen(*rstr)),strlen(rtitle))) else rw = -1 if (*cstr!="") wd = colmax((strlen(*cstr)\wd)) :+ 2 else wd = wd :+ 2
sfmt = "%" :+ strofreal(wd) :+ "s"
if (save) { _mm_matlist(X, fmt, b, *rstr, *cstr, rtitle, wd, rw, sfmt, 1, cols(X), 1, 1, res) return }
lw = st_numscalar("c(linesize)") - ((2+rw+2*(b>0))+2*(b==1)) if ((sum(wd)+cols(X))<=lw) { _mm_matlist(X, fmt, b, *rstr, *cstr, rtitle, wd, rw, sfmt, 1, cols(X), 1, 1) return }
j0 = 1 l = 1 r = 0 while (j0<=cols(X)) { j1 = j0 while (sum(wd[|j0\j1|]:+1)<=lw) { if (++j1>cols(X)) { r = 1 break } } j1 = max((j1-1,j0)) _mm_matlist(X, fmt, b, *rstr, *cstr, rtitle, wd, rw, sfmt, j0, j1, l, r) l = 0 j0 = j1+1 } }
string colvector mm_smatlist( real matrix X, | string matrix fmt, real scalar b, string vector rstripe, string vector cstripe, string scalar rtitle ) { string colvector res
if (args()<2) fmt = "%9.0g" if (args()<3) b = 1 if (args()<4) rstripe = J(0,1,"") if (args()<5) cstripe = J(1,0,"") if (args()<5) rtitle = "" mm_matlist(X, fmt, b, rstripe, cstripe, rtitle, res) return(res) }
void _mm_matlist( real matrix X, string matrix fmt, real scalar b, string vector rstripe, string vector cstripe, string scalar rtitle, real rowvector wd, real scalar rw, string rowvector sfmt, real scalar j0, real scalar j1, real scalar l, real scalar r, | transmorphic matrix res ) { real scalar i, j, di, ri string scalar rfmt, line pointer(real scalar) scalar fi, fj
di = args()<14 fi = (rows(fmt)!=1 ? &i : &1) fj = (cols(fmt)!=1 ? &j : &1) rfmt = "%"+strofreal(rw)+"s"
if (di==0) { ri = 0 res = J(rows(X)+(1+(b==3))*(cstripe!="")+(b>0)+(b==1) +(b==3)+(l==0&(b==0|b==2)),1,"") }
if (l==0 & (b==0 | b==2)) { if (di) printf("\n") else res[++ri] = "" } if (cstripe!="") { if (b==3) { line = sprintf("{txt}" + "{hline " + strofreal(2+rw+1) + "}{c TT}" + "{hline " + strofreal(sum(wd[|j0\j1|]:+1)-1) + "}") if (di) display(line) else res[++ri] = line } line = sprintf("{txt} ") if (rstripe!="") line = line + sprintf(rfmt+" ", rtitle) if (b>0) line = line + sprintf((b==1 ? " " : "{c |}")) for (j=j0;j<=j1;j++) { line = line + sprintf(sfmt[j], cstripe[j]) if (j<j1) line = line + sprintf(" ") } if (di) display(line) else res[++ri] = line } if (b) { line = sprintf("{txt}" + (b==1 ? (2+rw+1)*" " + (l ? "{c TLC}" : " ") : "{hline " + strofreal(2+rw+1) + "}{c +}") + "{hline " + strofreal(sum(wd[|j0\j1|]:+1)+1-2*(b!=1)) + "}" + (b==1 ? (r ? "{c TRC}" : "") : "")) if (di) display(line) else res[++ri] = line } for (i=1;i<=rows(X);i++) { line = sprintf("{txt} ") if (rstripe!="") line = line + sprintf(rfmt+" ", rstripe[i]) if (b>0) line = line + sprintf((l | b!=1 ? "{c |}" : " ")) line = line + sprintf("{res}") for (j=j0;j<=j1;j++) { line = line + sprintf(sfmt[j],sprintf(fmt[*fi,*fj], X[i,j])) if (j<j1) line = line + sprintf(" ") } if (b==1 & r) line = line + sprintf(" {txt}{c |}") if (di) display(line) else res[++ri] = line } if (b==1|b==3) { line = sprintf("{txt}" + (b==1 ? (2+rw+1)*" " + (l ? "{c BLC}" : " ") : "{hline " + strofreal(2+rw+1) + "}{c BT}") + "{hline " + strofreal(sum(wd[|j0\j1|]:+1)+1-2*(b!=1)) + "}" + (b==1 ? (r ? "{c BRC}" : "") : "")) if (di) display(line) else res[++ri] = line } }
end
*! *! mm_insheet.mata *! version 2.0.1, Ben Jann, 18mar2006 *! based on cat(), version 2.0.0 23jan2006 version 9.0 mata:
string matrix mm_insheet(string scalar filename, | string scalar del, real scalar uline1, real scalar line2) { real scalar fh, fpos, line1, i, r, c, j, delpos, ji string matrix EOF string matrix res string scalar line
/* ------------------------------------------------------------ */ /* setup */ if (del=="") del = char(9) fh = fopen(filename, "r") EOF = J(0, 0, "") line1 = floor(uline1) if (line1<1 | line1>=.) line1 = 1
/* ------------------------------------------------------------ */ /* nothing to read case */ if (line1>line2) { fclose(fh) return(J(0, 0, "")) }
/* ------------------------------------------------------------ */ /* skip opening lines */ for (i=1; i<line1; i++) { if (fget(fh)==EOF) { fclose(fh) return(J(0, 0, "")) } } fpos = ftell(fh)
/* ------------------------------------------------------------ */ /* count lines and columns */ j = 0 for (i=line1; i<=line2; i++) { if ((line=fget(fh))==EOF) break ji = 1 while (delpos = strpos(line, del)) { ji++ line = substr(line, delpos+1, .) } if (ji>j) j = ji } res = J(r = i-line1, c = j, "")
/* ------------------------------------------------------------ */ /* read lines */ fseek(fh, fpos, -1) for (i=1; i<=r; i++) { if ((line=fget(fh))==EOF) { /* unexpected EOF -- file must have changed */ fclose(fh) if (--i) return(res[|1,1 \ i,c|]) return(J(0,c,"")) } for (j=1; j<=c; j++) { delpos = strpos(line, del) if (delpos==0) { res[i,j] = line break } res[i,j] = substr(line, 1, delpos-1) line = substr(line, delpos+1, .) } } fclose(fh) return(res) }
end
*! *! mm_infile.mata *! version 2.0.2, Ben Jann, 18mar2006 *! based on cat(), version 2.0.0 23jan2006 version 9.0 mata:
string matrix mm_infile(string scalar filename, | real scalar uline1, real scalar line2) { real scalar fh, fpos, line1, i, r, c, j string matrix EOF string matrix res string scalar line
/* ------------------------------------------------------------ */ /* setup */ fh = fopen(filename, "r") EOF = J(0, 0, "") line1 = floor(uline1) if (line1<1 | line1>=.) line1 = 1
/* ------------------------------------------------------------ */ /* nothing to read case */ if (line1>line2) { fclose(fh) return(J(0, 0, "")) }
/* ------------------------------------------------------------ */ /* skip opening lines */ for (i=1; i<line1; i++) { if (fget(fh)==EOF) { fclose(fh) return(J(0, 0, "")) } } fpos = ftell(fh)
/* ------------------------------------------------------------ */ /* count lines and columns */ j = 0 for (i=line1; i<=line2; i++) { if ((line=fget(fh))==EOF) break line = tokens(line) if (cols(line)>j) j = cols(line) } res = J(r = i-line1, c = j, "")
/* ------------------------------------------------------------ */ /* read lines */ fseek(fh, fpos, -1) for (i=1; i<=r; i++) { if ((line=fget(fh))==EOF) { /* unexpected EOF -- file must have changed */ fclose(fh) if (--i) return(res[|1,1 \ i,c|]) return(J(0,c,"")) } line = tokens(line) res[|i,1 \ i,cols(line)|] = line } fclose(fh) return(res) }
end
*! *! mm_outsheet.mata *! version 1.0.4, Ben Jann, 14apr2006 version 9.0 mata:
void mm_outsheet(string scalar fn, string matrix s, | string scalar mode0, string scalar del) { string scalar line, mode, m real scalar i, j, fh
if (args()<4) del = char(9) mode = mm_strexpand(mode0, ("append", "replace")) m = "w" if (mode=="replace") unlink(fn) else if (mode=="append") m = "a" fh = fopen(fn, m) for (i=1; i<=rows(s); i++) { line = J(1,1,"") for (j=1; j<=cols(s); j++) { line = line + s[i,j] if (j<cols(s)) line = line + del } fput(fh, line) } fclose(fh) } end
*! *! mm_callf.mata *! version 1.0.1, Ben Jann, 16jun2006 version 9.2
local o1 "o1" local po1 "*p.o1" forv i=2/10 { local o`i' "`o`=`i'-1'', o`i'" local po`i' "`po`=`i'-1'', *p.o`i'" }
mata:
struct mm_callf_o10 { pointer(function) scalar f real scalar n pointer scalar `o10' }
struct mm_callf_o10 scalar mm_callf_setup( pointer(function) scalar f, real scalar n, | `o10') { struct mm_callf_o10 scalar p
p.f = f p.n = n p.o1 = &o1 ; p.o2 = &o2 ; p.o3 = &o3 ; p.o4 = &o4 ; p.o5 = &o5 p.o6 = &o6 ; p.o7 = &o7 ; p.o8 = &o8 ; p.o9 = &o9 ; p.o10 = &o10 return(p) }
transmorphic mm_callf(struct mm_callf_o10 p, | `o5') { if (args()==1) return(_mm_callf0(p)) if (args()==2) return(_mm_callf1(p, `o1')) if (args()==3) return(_mm_callf2(p, `o2')) if (args()==4) return(_mm_callf3(p, `o3')) if (args()==5) return(_mm_callf4(p, `o4')) return(_mm_callf5(p, `o5')) }
transmorphic _mm_callf0(struct mm_callf_o10 p) { if (p.n<=0) return( (*p.f)() ) if (p.n==1) return( (*p.f)(`po1') ) if (p.n==2) return( (*p.f)(`po2') ) if (p.n==3) return( (*p.f)(`po3') ) if (p.n==4) return( (*p.f)(`po4') ) if (p.n==5) return( (*p.f)(`po5') ) if (p.n==6) return( (*p.f)(`po6') ) if (p.n==7) return( (*p.f)(`po7') ) if (p.n==8) return( (*p.f)(`po8') ) if (p.n==9) return( (*p.f)(`po9') ) return( (*p.f)(`po10') ) }
transmorphic _mm_callf1(struct mm_callf_o10 p, `o1') { if (p.n<=0) return( (*p.f)(`o1') ) if (p.n==1) return( (*p.f)(`o1', `po1') ) if (p.n==2) return( (*p.f)(`o1', `po2') ) if (p.n==3) return( (*p.f)(`o1', `po3') ) if (p.n==4) return( (*p.f)(`o1', `po4') ) if (p.n==5) return( (*p.f)(`o1', `po5') ) if (p.n==6) return( (*p.f)(`o1', `po6') ) if (p.n==7) return( (*p.f)(`o1', `po7') ) if (p.n==8) return( (*p.f)(`o1', `po8') ) if (p.n==9) return( (*p.f)(`o1', `po9') ) return( (*p.f)(`o1', `po10') ) }
transmorphic _mm_callf2(struct mm_callf_o10 p, `o2') { if (p.n<=0) return( (*p.f)(`o2') ) if (p.n==1) return( (*p.f)(`o2', `po1') ) if (p.n==2) return( (*p.f)(`o2', `po2') ) if (p.n==3) return( (*p.f)(`o2', `po3') ) if (p.n==4) return( (*p.f)(`o2', `po4') ) if (p.n==5) return( (*p.f)(`o2', `po5') ) if (p.n==6) return( (*p.f)(`o2', `po6') ) if (p.n==7) return( (*p.f)(`o2', `po7') ) if (p.n==8) return( (*p.f)(`o2', `po8') ) if (p.n==9) return( (*p.f)(`o2', `po9') ) return( (*p.f)(`o2', `po10') ) }
transmorphic _mm_callf3(struct mm_callf_o10 p, `o3') { if (p.n<=0) return( (*p.f)(`o3') ) if (p.n==1) return( (*p.f)(`o3', `po1') ) if (p.n==2) return( (*p.f)(`o3', `po2') ) if (p.n==3) return( (*p.f)(`o3', `po3') ) if (p.n==4) return( (*p.f)(`o3', `po4') ) if (p.n==5) return( (*p.f)(`o3', `po5') ) if (p.n==6) return( (*p.f)(`o3', `po6') ) if (p.n==7) return( (*p.f)(`o3', `po7') ) if (p.n==8) return( (*p.f)(`o3', `po8') ) if (p.n==9) return( (*p.f)(`o3', `po9') ) return( (*p.f)(`o3', `po10') ) }
transmorphic _mm_callf4(struct mm_callf_o10 p, `o4') { if (p.n<=0) return( (*p.f)(`o4') ) if (p.n==1) return( (*p.f)(`o4', `po1') ) if (p.n==2) return( (*p.f)(`o4', `po2') ) if (p.n==3) return( (*p.f)(`o4', `po3') ) if (p.n==4) return( (*p.f)(`o4', `po4') ) if (p.n==5) return( (*p.f)(`o4', `po5') ) if (p.n==6) return( (*p.f)(`o4', `po6') ) if (p.n==7) return( (*p.f)(`o4', `po7') ) if (p.n==8) return( (*p.f)(`o4', `po8') ) if (p.n==9) return( (*p.f)(`o4', `po9') ) return( (*p.f)(`o4', `po10') ) }
transmorphic _mm_callf5(struct mm_callf_o10 p, `o5') { if (p.n<=0) return( (*p.f)(`o5') ) if (p.n==1) return( (*p.f)(`o5', `po1') ) if (p.n==2) return( (*p.f)(`o5', `po2') ) if (p.n==3) return( (*p.f)(`o5', `po3') ) if (p.n==4) return( (*p.f)(`o5', `po4') ) if (p.n==5) return( (*p.f)(`o5', `po5') ) if (p.n==6) return( (*p.f)(`o5', `po6') ) if (p.n==7) return( (*p.f)(`o5', `po7') ) if (p.n==8) return( (*p.f)(`o5', `po8') ) if (p.n==9) return( (*p.f)(`o5', `po9') ) return( (*p.f)(`o5', `po10') ) }
end