* version 1.2.0 13Dec2011 MLB * version 1.0.0 14Nov2011 MLB mata void function marg_invert(real scalar l, real scalar u, real scalar e, struct margdata scalar data, real vector x, pointer scalar f, pointer scalar F, string scalar Q, | tl, tu) { real matrix res real vector o, op if (args() == 8 ) { res = marg_initknots(l,u,e,data, f, F) } else if (args() == 9) { res = marg_initknots(l,u,e,data, f, F, tl) } else if (args() == 10) { res = marg_initknots(l,u,e,data, f, F, tl, tu) } res = marg_modeknots(res, data, f, F) res = marg_knots(res,e,data, f, F) o = order(x, 1) op = invorder(o) st_store(.,Q,st_local("touse"),(marg_eval(res,x[o]))[op] ) } mata mlib add lmargdistfit marg_invert() real matrix function marg_eval(real matrix hermite, real vector x){ real scalar n, m, t, j, i, m0 real vector y real matrix temp n = length(x) m = rows(hermite) y = J(n,1,.) /* ignore values in x that are less than the lowest value in hermite */ t = hermite[1,2] for (m0=j=1;j<=n;j++) { if (x[j]>=t) break } /* begin approximations */ for (i=j;i<=n;i++) { for (j=m0;j<m;j++) { /* sic */ if (x[i]>=hermite[j,2] & x[i]<hermite[j+1,2]) break } if (j>=m) { if (x[i]==hermite[m,2]) j = m-1 } if (j<m) { temp = hermite[|j,1\j+1,6|] y[i] = marg_as_subeval(x[i],temp) } m0 = j } return(y) } mata mlib add lmargdistfit marg_eval() real matrix function marg_as(real matrix puf) { real scalar p, u, f real matrix res p=1 u=2 f=3 res = J(1,3,.) res[1,1] = (puf[2,u]-puf[1,u])/puf[1,f] res[1,2] = 3*(puf[2,p]-puf[1,p]) - (puf[2,u]-puf[1,u])*(2/puf[1,f]+1/puf[2,f]) res[1,3] = 2*(puf[1,p]-puf[2,p]) + (puf[2,u]-puf[1,u])*(1/puf[1,f]+1/puf[2,f]) return(res) } mata mlib add lmargdistfit marg_as() real scalar function marg_as_subeval(real scalar u, real matrix pufaaa) { real scalar a0, ui, a1, a2, a3, ut real matrix res a0 = 1 // = p ui = 2 a1 = 4 a2 = 5 a3 = 6 ut = (u - pufaaa[1,ui])/(pufaaa[2,ui]-pufaaa[1,ui]) res = pufaaa[1,a0] + pufaaa[1,a1]*ut + pufaaa[1,a2]*ut^2+ pufaaa[1,a3]*ut^3 return(res) } mata mlib add lmargdistfit marg_as_subeval() real matrix function marg_initknots(real scalar l, real scalar u, real scalar e, struct margdata scalar data , pointer scalar f, pointer scalar F, | real scalar tl , real scalar tu) { real scalar i, ok, diff, p1, p2, p, FF real matrix res, ins res = l \ u res = res,(*F)(data, res) // split till each interval is less than .05 i = 1 do { if (res[i+1,2]-res[i,2] > .05) { ins = res[i,1] + (res[i+1,1] - res[i,1]):/2 ins = ins, (*F)(data,ins) res = res[|1,1 \i,2|] \ ins \ res[|i+1,1 \ rows(res), 2|] } else { i = i + 1 } } while (i < rows(res) ) // expand till cover (e, 1-e) if (res[1,2] > e) { ok = 0 diff = (abs(res[2,1] - res[1,1])) do { ins = res[1,1] - diff if (args() >= 7) { if (ins < tl) { ins = tl ok = 1 } } res = (ins, (*F)(data, ins)) \ res if (res[1,2] < e) { ok = 1 } } while (ok == 0) } if (res[rows(res),2] < (1-e)) { ok = 0 diff = (abs(res[rows(res),1] - res[rows(res)-1,1])) do { ins = res[rows(res),1] + diff if (args() == 8) { if (ins>tu) { ins = tu ok = 1 } } res = res \ (ins, (*F)(data, ins)) if (res[rows(res),2] > (1-e)) { ok = 1 } } while (ok == 0) } // ensure largest and smallest value aren't too small (F<.1e, F>1-.1e) // but also not to large (F>e, F<1-e) if (res[1,2] < .01*e) { p1 = res[1,1] p2 = res[2,1] ok = 0 do { p = (p2-p1)/2 + p1 FF = (*F)(data,p) if (FF < .1*e) { p1 = p p2 = p2 } else if (FF > e) { p1 = p1 p2 = p } else { res[|1,1\1,2|] = (p , FF) ok = 1 } } while(ok == 0) } if (res[rows(res),2] > (1-.01*e)) { p1 = res[(rows(res)-1),1] p2 = res[rows(res),1] ok = 0 do { p = (p2-p1)/2 + p1 FF = (*F)(data,p) if (FF > (1-.1*e)) { p1 = p1 p2 = p } else if (FF < (1-e)) { p1 = p p2 = p2 } else { res[|rows(res),1\rows(res),2|] = (p , FF) ok = 1 } } while(ok == 0) } // add density function res = res, (*f)(data, res[.,1]) return(res) } mata mlib add lmargdistfit marg_initknots() real matrix function marg_knots(real matrix res, real scalar e, struct margdata scalar data, pointer scalar f, pointer scalar F ) { real scalar i, c, u, p real matrix temp, ins i = 1 res = res, J(rows(res),3,.) do { temp = res[|i,1\(i+1), 3|] temp = temp, ( marg_as(temp) \ J(1,3,.) ) c = 3*(temp[2,1] - temp[1,1])/(temp[2,2] - temp[1,2]) u = (temp[2,2] - temp[1,2])/2 + temp[1,2] if ((reldif((*F)(data,marg_as_subeval(u,temp)), u) < e ) & 1/temp[1,3]<= c & 1/temp[2,3]<=c) { res[|i,4 \ i, 6|] = temp[|1,4\1,6|] i = i + 1 } else { p = (temp[2,1] - temp[1,1])/2 + temp[1,1] ins = p, (*F)(data,p), (*f)(data,p), ., ., . res = res[|1,1 \ i, 6|] \ ins \ res[|i+1, 1 \ rows(res),6|] } } while ( i <= (rows(res)-1) ) return(res) } mata mlib add lmargdistfit marg_knots() end