*! Bottai M, Orsini N 27jun16

capture program drop  mi_impute_cmd_quantile
program mi_impute_cmd_quantile, sortpreserve
version 14.1
  
tempvar p1 p0 pa pb j s t y imputed observed  

qui gen `observed' = $MI_IMPUTE_user_miss != 1
qui gen double `imputed' = .

qui count if `observed' 
local N = r(N) 

qui gen double `y' = $MI_IMPUTE_user_ivar

local df = 3
 
qui foreach v of global MI_IMPUTE_user_xvars {
	sort `v'
	gen `t' = 1 in 1
	replace `t' = `t'[_n-1] + (`v'!= `v'[_n-1]) in 2/l
	
	if (`t'[_N] > (`df'+2)) {
		mkspline `s'_`v'  = `v', cubic nk(`df')
	}
	else {
		tab `v', gen(`s'_`v')
	}
	capture drop `t'
}

sort `observed' `y', stable
by `observed': gen `p1' = _n/_N
qui by `observed' `y': replace `p1' = `p1'[_N]
qui by `observed': gen `p0' = 0 if _n==1
qui by `observed': replace `p0' = cond(`p1'==`p1'[_n-1],`p0'[_n-1],`p1'[_n-1]) if _n>1
qui mlexp (ln(ibeta(exp({a:`s'_* _cons}),exp({b:`s'_* _cons}),`p1')-ibeta(exp({a:}),exp({b:}),`p0'))) $MI_IMPUTE_user_weight if $MI_IMPUTE_user_touse & `observed'
qui predictnl `pa' = exp(xb(a))  
qui predictnl `pb' = exp(xb(b))  
gen `j' = ceil(invibeta(`pa',`pb',runiform())*`N')  
qui putmata  `y' if  `observed' , view replace 
qui putmata  `j' if  !`observed' , view replace  
mata: _000q_Q_ = `y'[`j',1]
qui getmata `imputed' = _000q_Q_  , replace force
qui replace `imputed' = `y' if `imputed' == .
qui replace $MI_IMPUTE_user_ivar  = `imputed' if $MI_IMPUTE_user_touse
ereturn clear
end