*! 1.2.0 NJC 2sep2015
* 1.1.0 NJC 27aug2015
* 1.0.0 NJC 20jul2015
* ipolate 1.3.3 21sep2004
program mipolate, byable(onecall) sort
version 12
syntax varlist(min=2 max=2 numeric) [if] [in], ///
GENerate(string) ///
[ BY(varlist) ///
Epolate ///
Linear ///
Cubic ///
Spline ///
Idw ///
Idw2(numlist min=1 max=1 >=0) ///
Pchip ///
Forward ///
Backward ///
Nearest ///
Groupwise ///
ties(str) ]
// syntax checks
if _by() {
if "`by'" != "" {
di as err /*
*/ "option by() may not be combined with by prefix"
exit 190
}
local by "`_byvars'"
}
if "`ties'" != "" & "`nearest'" == "" {
di as err "ties() applies only with nearest"
exit 198
}
// ties():
// indulge upper case:
// any abbreviations of after, before, minimum, maximum
if "`ties'" != "" {
local ties = lower("`ties'")
local nchar = length("`ties'")
local OK = 0
if "`ties'" == substr("after", 1, `nchar') {
local ties "after"
local OK 1
}
else if "`ties'" == substr("before", 1, `nchar') {
local ties "before"
local OK 1
}
else if `nchar' == 1 & "`ties'" == "m" {
di as err "m ambiguous for ties() option"
exit 198
}
else if "`ties'" == substr("minimum", 1, `nchar') {
local ties "min"
local OK 1
}
else if "`ties'" == substr("maximum", 1, `nchar') {
local ties "max"
local OK 1
}
if !`OK' {
di as err "invalid ties() option: see help"
exit 198
}
}
if "`idw2'" != "" local idw "idw"
local nopts : word count ///
`cubic' `spline' `pchip' `idw' `nearest' `linear' ///
`forward' `backward' `groupwise'
if `nopts' > 1 {
di as err "must specify one interpolation method only"
}
else if `nopts' == 0 local linear "linear"
confirm new var `generate'
tokenize `varlist'
args usery x
local flag = 0
quietly {
tempvar y
marksample touse, novarlist
replace `touse' = 0 if `x' >= .
count if `touse'
if r(N) == 0 error 2000
bysort `touse' `by' `x': ///
gen double `y' = sum(`usery') / sum(`usery' < .) if `touse'
by `touse' `by' `x': replace `y' = `y'[_N]
if "`linear'" != "" {
tempvar negy negx xhi yhi xlo ylo m b z
gen double `negx' = -`x'
gen double `negy' = -`y'
sort `touse' `by' `negx' `negy'
gen double `xhi' = `x' if `touse'
gen double `yhi' = `y' if `touse'
by `touse' `by': replace `yhi' = `yhi'[_n-1] if `y' >= . & `touse'
by `touse' `by': replace `xhi' = `xhi'[_n-1] if `y' >= . & `touse'
sort `touse' `by' `x' `y'
gen double `xlo' = `x' if `touse'
gen double `ylo' = `y' if `touse'
by `touse' `by': replace `ylo' = `ylo'[_n-1] if `y' >= . & `touse'
by `touse' `by': replace `xlo' = `xlo'[_n-1] if `y' >= . & `touse'
gen double `m' = (`yhi'-`ylo')/(`xhi'-`xlo')
drop `yhi' `xhi'
gen double `b' = `ylo' - `m'*`xlo'
drop `ylo' `xlo' `negx' `negy'
gen double `z' = `y' if `touse'
replace `z' = `m'*`x' + `b' if `touse' & `z' >= .
}
else if "`cubic'" != "" {
tempvar negx negy ok ok2 x1 x2 x3 x4 y1 y2 y3 y4 m m1 m2 m3 m4 z
* following values
gen double `negx' = -`x'
gen double `negy' = -`y'
bysort `touse' `by' (`negx' `negy') : gen `ok' = _n * (`y' < .)
by `touse' `by' : replace `ok' = `ok'[_n-1] if !`ok'
by `touse' `by' : gen `ok2' = `ok'[_n-1] * (`ok' > `ok'[_n-1])
by `touse' `by' : replace `ok2' = `ok2'[_n-1] if !`ok2'
by `touse' `by' : gen double `x4' = `x'[`ok2']
by `touse' `by' : gen double `y4' = `y'[`ok2']
by `touse' `by' : gen double `x3' = `x'[`ok']
by `touse' `by' : gen double `y3' = `y'[`ok']
* preceding values
bysort `touse' `by' (`x' `y'): replace `ok' = _n * (`y' < .)
by `touse' `by' : replace `ok' = `ok'[_n-1] if !`ok'
by `touse' `by' : replace `ok2' = `ok'[_n-1] * (`ok' > `ok'[_n-1])
by `touse' `by' : replace `ok2' = `ok2'[_n-1] if !`ok2'
by `touse' `by' : gen double `x1' = `x'[`ok2']
by `touse' `by' : gen double `y1' = `y'[`ok2']
by `touse' `by' : gen double `x2' = `x'[`ok']
by `touse' `by' : gen double `y2' = `y'[`ok']
gen double `m1' = (`x' - `x2') * (`x' - `x3') * (`x' - `x4') /*
*/ / ((`x1' - `x2') * (`x1' - `x3') * (`x1' - `x4'))
gen double `m2' = (`x' - `x1') * (`x' - `x3') * (`x' - `x4') /*
*/ / ((`x2' - `x1') * (`x2' - `x3') * (`x2' - `x4'))
gen double `m3' = (`x' - `x1') * (`x' - `x2') * (`x' - `x4') /*
*/ / ((`x3' - `x1') * (`x3' - `x2') * (`x3' - `x4'))
gen double `m4' = (`x' - `x1') * (`x' - `x2') * (`x' - `x3') /*
*/ / ((`x4' - `x1') * (`x4' - `x2') * (`x4' - `x3'))
gen double `m' = /*
*/ `m1' * `y1' + `m2' * `y2' + `m3' * `y3' + `m4' * `y4'
gen double `z' = `y' if `touse'
replace `z' = `m' if `touse' & `z' == .
}
else if "`nearest'" != "" {
tempvar negx prevy prevx nexty nextx z
* values before
gen double `prevy' = `y' if `touse' & `y' < .
gen double `prevx' = `x' if `touse' & `y' < .
gen double `nexty' = `prevy'
gen double `nextx' = `prevx'
bysort `touse' `by' (`x'): replace `prevy' = `prevy'[_n - 1] if `prevy' == .
by `touse' `by': replace `prevx' = `prevx'[_n - 1] if `prevx' == .
* values after
gen double `negx' = -`x'
bysort `touse' `by' (`negx') : replace `nexty' = `nexty'[_n - 1] if `nexty' == .
by `touse' `by' : replace `nextx' = `nextx'[_n - 1] if `nextx' == .
* interpolation
gen double `z' = `y' if `touse'
replace `z' = `nexty' if (`nextx' - `x') < (`x' - `prevx') & `z' == . & `touse'
replace `z' = `prevy' if (`x' - `prevx') < (`nextx' - `x') & `z' == . & `touse'
if "`ties'" != "" {
if "`ties'" == "after" {
replace `z' = `nexty' if (`nextx' - `x') == (`x' - `prevx') & `z' == . & `touse'
}
else if "`ties'" == "before" {
replace `z' = `prevy' if (`nextx' - `x') == (`x' - `prevx') & `z' == . & `touse'
}
else if "`ties'" == "min" {
replace `z' = min(`nexty', `prevy') if (`nextx' - `x') == (`x' - `prevx') & `z' == . & `touse'
}
else if "`ties'" == "max" {
replace `z' = max(`nexty', `prevy') if (`nextx' - `x') == (`x' - `prevx') & `z' == . & `touse'
}
}
else replace `z' = (`nexty' + `prevy')/2 if (`nextx' - `x') == (`x' - `prevx') & `z' == . & `touse'
local epolate
}
else if "`spline'" != "" {
tempvar group z guse zmiss negy negx
// only use one of any repeated x
by `touse' `by' `x': gen byte `guse' = `touse' & (_n == 1)
sort `guse' `by' `x'
if "`by'" != "" {
bysort `guse' `by': gen byte `group' = _n == 1 & `guse'
by `guse': replace `group' = sum(`group')
}
else gen byte `group' = 1
su `group', meanonly
local ng = r(max)
gen double `z' = `y' if `touse'
gen byte `zmiss' = 0
forval g = 1/`ng' {
replace `guse' = `group' == `g'
replace `zmiss' = `guse' & missing(`z')
count if `zmiss'
if r(N) > 1 {
mata : csipolate("`y'", "`x'", "`guse'", "`z'", "`zmiss'")
}
}
// copy interpolated y to repeated x
bysort `touse' `by' `x' (`z') : ///
replace `z' = `z'[1] if `touse' & `z' == .
}
else if "`pchip'" != "" {
tempvar group z guse zmiss
// only use one of any repeated x
by `touse' `by' `x': gen byte `guse' = `touse' & (_n == 1)
sort `guse' `by' `x'
if "`by'" != "" {
by `guse' `by': gen byte `group' = _n == 1 & `guse'
by `guse': replace `group' = sum(`group')
}
else gen byte `group' = 1
su `group', meanonly
local ng = r(max)
gen double `z' = `y' if `touse'
gen byte `zmiss' = 0
forval g = 1/`ng' {
replace `guse' = `group' == `g'
replace `zmiss' = `guse' & missing(`z')
count if `guse' & !missing(`z')
if r(N) > 2 {
mata : ///
pchipolate("`y'", "`x'", "`guse'", "`z'", "`zmiss'")
}
else local flag = 1
}
// copy interpolated y to repeated x
bysort `touse' `by' `x' (`z') : ///
replace `z' = `z'[1] if `touse' & `z' == .
local epolate
}
else if "`idw'`idw2'" != "" {
tempvar group z guse zmiss
// only use one of any repeated x
by `touse' `by' `x': gen byte `guse' = `touse' & (_n == 1)
sort `guse' `by' `x'
if "`by'" != "" {
bysort `guse' `by': gen byte `group' = _n == 1 & `guse'
by `guse': replace `group' = sum(`group')
}
else gen byte `group' = 1
su `group', meanonly
local ng = r(max)
gen double `z' = `y' if `touse'
gen byte `zmiss' = 0
if "`idw2'" == "" local idw2 = 2
forval g = 1/`ng' {
replace `guse' = `group' == `g'
replace `zmiss' = `guse' & missing(`z')
count if `zmiss'
mata : idwipolate("`y'", "`x'", "`guse'", "`z'", "`zmiss'", `idw2')
}
// copy interpolated y to repeated x
bysort `touse' `by' `x' (`z') : ///
replace `z' = `z'[1] if `touse' & `z' == .
local epolate
}
else if "`forward'`backward'" != "" {
tempvar z
gen double `z' = `y'
if "`forward'" != "" {
bysort `touse' `by' (`x') : ///
replace `z' = `z'[_n-1] if `touse' & `z' == .
}
else {
tempvar negx
gen double `negx' = -`x'
bysort `touse' `by' (`negx'): ///
replace `z' = `z'[_n-1] if `touse' & `z' == .
}
local epolate
}
else if "`groupwise'" != "" {
tempvar z OK
bysort `touse' `by' (`y') : ///
gen byte `OK' = (`y' == `y'[1]) | missing(`y')
count if `touse' & !`OK'
if r(N) {
if "`by'" != "" local text "within groups"
di as err ///
"different non-missing values of `usery' `text'"
exit 498
}
gen double `z' = `y'
by `touse' `by': ///
replace `z' = `z'[1] if missing(`z') & `touse'
local epolate
}
if `"`epolate'"' != "" {
sort `touse' `by' `x' `z' /* already sorted */
tempvar m b M B ismiss
by `touse' `by': gen double /*
*/ `m' = (`z'[_n+1] - `z')/(`x'[_n+1] - `x')
gen double `b' = `z' - `m'*`x'
gen double `M' = `m'
gen double `B' = `b'
by `touse' `by': replace `m' = `m'[_n-1] if `m'[_n-1]< .
by `touse' `by': replace `b' = `b'[_n-1] if `b'[_n-1]< .
gen byte `ismiss' = `z' >= .
by `touse' `by': replace `ismiss' = 0 /*
*/ if _n>1 & `ismiss'[_n-1]==0
by `touse' `by': replace `z' = `m'[_N]*`x' + `b'[_N] /*
*/ if `touse' & `ismiss'
drop `ismiss' `m' `b'
gen double `negx' = -`x'
gen double `negy' = -`z'
sort `touse' `by' `negx' `negy'
by `touse' `by': replace `M' = `M'[_n-1] if `M'[_n-1] < .
by `touse' `by': replace `B' = `B'[_n-1] if `B'[_n-1] < .
gen byte `ismiss' = `z' >= .
by `touse' `by': replace `ismiss' = 0 /*
*/ if _n>1 & `ismiss'[_n-1] == 0
by `touse' `by': replace `z' = `M'[_N]*`x' + `B'[_N] /*
*/ if `touse' & `ismiss'
}
rename `z' `generate'
compress `generate'
count if `generate' >= .
}
if `flag' {
di as txt "note: at least 3 values needed in any interpolation"
}
if r(N) > 0 {
if r(N) != 1 local pl "s"
di as txt "(" r(N) `" missing value`pl' generated)"'
}
end
mata:
void idwipolate(string scalar yvarname,
string scalar xvarname,
string scalar tousename,
string scalar zvarname,
string scalar zmissname,
real scalar power
) {
real matrix xyvar
real colvector x, y, where, newz, weight
real scalar i
st_view(xyvar, ., (xvarname, yvarname), tousename)
x = select(xyvar[,1], (xyvar[,2] :< .))
y = select(xyvar[,2], (xyvar[,2] :< .))
newz = where = select(xyvar[,1], (xyvar[,2] :== .))
for(i = 1; i <= rows(where); i++) {
weight = abs(x :- where[i]):^(-power)
newz[i] = sum(y :* weight) / sum(weight)
}
st_store(., zvarname, zmissname, newz)
}
void csipolate(string scalar yvarname,
string scalar xvarname,
string scalar tousename,
string scalar zvarname,
string scalar zmissname
) {
real matrix xyvar
real colvector x, y, where
st_view(xyvar, ., (xvarname, yvarname), tousename)
x = select(xyvar[,1], (xyvar[,2] :< .))
where = select(xyvar[,1], (xyvar[,2] :== .))
y = select(xyvar[,2], (xyvar[,2] :< .))
st_store(., zvarname, zmissname, spline3eval(spline3(x, y), where))
}
void pchipolate(string scalar yvarname,
string scalar xvarname,
string scalar tousename,
string scalar zvarname,
string scalar zmissname
) {
real matrix xyvar
real colvector x, y, where
st_view(xyvar, ., (xvarname, yvarname), tousename)
x = select(xyvar[,1], (xyvar[,2] :< .))
where = select(xyvar[,1], (xyvar[,2] :== .))
y = select(xyvar[,2], (xyvar[,2] :< .))
st_store(., zvarname, zmissname, pchip(x, y, where))
}
real colvector pchip(real colvector x, real colvector y, real colvector u)
{
real scalar n, nu, k, j
real colvector h, delta, d, c, b, which, s
n = length(x)
h = x[2::n] - x[1::n-1]
delta = (y[2::n] - y[1::n-1]) :/ h
d = pchipslopes(h, delta)
c = (3*delta - 2*d[1::n-1] - d[2::n]) :/ h
b = (d[1::n-1] - 2*delta + d[2::n]) :/ (h:^2)
nu = length(u)
k = J(nu, 1, 1)
for (j = 2; j <= n-1; j++) {
which = select((1::nu), x[j] :<= u)
k[which] = J(length(which), 1, j)
}
s = u - x[k]
return(y[k] + s :* (d[k] + s :* (c[k] + s :* b[k])))
}
real colvector function pchipslopes(real colvector h, real colvector delta) {
real scalar n
real colvector d, k, w1, w2
n = length(h) + 1
d = J(n, 1, 0)
k = 1 :+ select((1::n-2), sign(delta[1::n-2]) :* sign(delta[2::n-1]) :> 0)
w1 = 2*h[k] + h[k:-1]
w2 = h[k] + 2*h[k:-1]
d[k] = (w1 + w2) :/ (w1 :/ delta[k:-1] + w2 :/ delta[k])
d[1] = pchipend(h[1], h[2], delta[1], delta[2])
d[n] = pchipend(h[n-1], h[n-2], delta[n-1], delta[n-2])
return(d)
}
real scalar function pchipend(h1, h2, del1, del2) {
real scalar d
d = ((2*h1 + h2)*del1 - h1*del2) / (h1 + h2)
if (sign(d) != sign(del1)) d = 0
else {
if (sign(del1) != sign(del2) & (abs(d) > abs(3*del1))) {
d = 3*del1
}
}
return(d)
}
end