* 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=hermite[j,2] & x[i]=m) { if (x[i]==hermite[m,2]) j = m-1 } if (j .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