* version 1.2.0 13Dec2011 MLB * version 1.0.0 14Nov2011 MLB mata real matrix function marg_hermitegen_ac(real vector lu, real scalar e, struct margdata scalar data, pointer scalar f, pointer scalar F ) { real matrix hermite hermite = marg_initknots_ac(lu,data, f, F) hermite = marg_modeknots(hermite, data, f, F) hermite = marg_knots_ac(hermite,e,data, f, F) return(hermite) } mata mlib add lmargdistfit marg_hermitegen_ac() real matrix function marg_eval_ac(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,1] 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,1] & x[i]=m) { if (x[i]==hermite[m,1]) 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) ) // add density function res = res, (*f)(data, res[.,1]) return(res) } mata mlib add lmargdistfit marg_initknots_ac() real matrix function marg_modeknots( real matrix res, struct margdata scalar data, pointer scalar f, pointer scalar F) { real scalar k, m, j, i, dens real matrix res2, temp, ins res2 = res for (k=1; k <= 4; k++) { res = res2 j = 1 for (i=1; i <= rows(res)-1; i++) { // guess of mode based on quadratic spline using info from CDF temp = res[|i,1\(i+1), 3|] temp = temp, ( marg_as_ac(temp) \ J(1,3,.) ) if (abs(3*temp[1,6]) > 1e-6) { // check for extrema only if enough curveture m = -temp[1,5]/(3*temp[1,6]) if (m > 1e-4 & m < .05*k) { // shift left knot m = res[i,1] + (res[i+1,1] - res[i,1])*m dens = (*f)(data,m) if (temp[1,6] <= 0) { // approximation finds maximum if ( dens > res[i,3] ) { // new is improvement if new is larger than old ins = m, (*F)(data,m), dens res2[|j,1 \ j,3|] = ins } } else { // approximation finds minimum if ( dens < res[i,3] ) { // new is improvement if new is smaller than old ins = m, (*F)(data,m), dens res2[|j,1 \ j,3|] = ins } } } else if (m < (1-1e-4) & m > (1-.05*k) ) { // shift right knot m = res[i,1] + (res[i+1,1] - res[i,1])*m dens = (*f)(data,m) if (temp[1,6] <= 0) { // approximation finds a maximum if ( dens > res[i+1,3] ) { // new is improvement if new is larger than old ins = m, (*F)(data,m), dens res2[|j+1,1 \ j+1,3|] = ins } } else { // approximation finds a minimum if ( dens < res[i+1,3] ) { // new is improvement if new is smaller than old ins = m, (*F)(data,m), dens res2[|j+1,1 \ j+1,3|] = ins } } } else if ( m >= .05*k & m <= (1-.05*k)) { // add a knot m = res[i,1] + (res[i+1,1] - res[i,1])*m ins = m, (*F)(data,m), (*f)(data,m) res2 = res2[|1,1 \j,3|] \ ins \ res2[|j+1,1 \ rows(res2), 3|] j = j + 1 } } j = j+1 } } return(res2) } mata mlib add lmargdistfit marg_modeknots() real matrix function marg_knots_ac(real matrix res, real scalar e, struct margdata scalar data, pointer scalar f, pointer scalar F ) { real scalar i, u, c real matrix temp, ins i = 1 res = res, J(rows(res),3,.) do { temp = res[|i,1\(i+1), 6|] temp[|1,4 \1,6|] = marg_as_ac(temp) u = (temp[2,1] - temp[1,1])/2 + temp[1,1] c = 3*(temp[2,2] - temp[1,2])/(temp[2,1] - temp[1,1]) if ( reldif(temp[1,2], 1) < e ) { res[|i,4 \ i, 6|] = J(1,3,0) i = i + 1 } else if ( reldif(marg_as_subeval_ac(u,temp), (*F)(data,u)) < e & temp[1,3]<= c & temp[2,3]<=c ) { res[|i,4 \ i, 6|] = temp[|1,4\1,6|] i = i + 1 } else { ins = u, (*F)(data,u), (*f)(data,u), ., ., . 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_ac() end