*! 0.1 HS, Sep 7, 2016 *! 0.2 KBN, Jan 9, 2019 pr define wtdtttdiag version 14.0 syntax varname [if] [in], /// [nbins(integer 0) ytitle(passthru) /// legend(passthru) nq(integer 300) /// lpattern(passthru) fcolor(passthru) /// replace *] qui{ tempname coefs tstart tend delta converged matrix `coefs' = e(b) local disttype = r(disttype) scalar `delta' = r(delta) local reverse = r(reverse) scalar `converged' = e(converged) if r(start) != . & r(end) != . { scalar `tstart' = r(start) scalar `tend' = r(end) } else { scalar `tstart' = r(samplestart) scalar `tend' = r(sampleend) } preserve tokenize `varlist' local obstime `1' if "`reverse'" == "" { replace `obstime' = `obstime' - `tstart' } if "`reverse'" == "reverse" { replace `obstime' = `tend' - `obstime' } local xtitletxt : var label `obstime' local xformat : format `obstime' if "`if'" != "" | "`in'" != "" { keep `if' `in' } count local nobs = r(N) if `nbins' == 0 { local nbins = int(min(sqrt(`nobs'), 10*log(`nobs')/log(10))) } local interwidth = `delta' / `nbins' tempvar lowcutlev obsdens tfit fitdens gen `lowcutlev' = int(`obstime' / `interwidth') * `interwidth' collapse (count) `obstime', by(`lowcutlev') gen `obsdens' = `obstime' / `nobs' / `interwidth' * / `delta' * `interwidth' la var `obsdens' "Observed" count local np1 = r(N) local np1 = `np1' + 1 set obs `np1' if "`reverse'" == "reverse" { replace `lowcutlev' = `delta' if missing(`lowcutlev') } else { replace `lowcutlev' = `tstart' + `delta' if missing(`lowcutlev') } format `lowcutlev' `xformat' sort `lowcutlev' local np200 = max(`np1', `nq') set obs `np200' gen `tfit' = (_n - 1) / (_N - 1) * `delta' format `tfit' `xformat' if "`disttype'" == "exp" { gen `fitdens' = (invlogit(`coefs'[1, 1]) /// * exp(- `tfit' * exp(`coefs'[1, 2])) /// * exp(`coefs'[1, 2]) /// + invlogit(- `coefs'[1, 1]) / `delta') } if "`disttype'" == "lnorm" { gen `fitdens' = (invlogit(`coefs'[1, 1]) /// * normal(- (log(`tfit') - `coefs'[1, 2]) /// / exp(`coefs'[1, 3])) /// / exp(`coefs'[1, 2] + exp(2 * `coefs'[1, 3]) / 2) /// + invlogit(- `coefs'[1, 1]) / `delta') replace `fitdens' = invlogit(`coefs'[1, 1]) /// / exp(`coefs'[1, 2] + exp(2 * `coefs'[1, 3]) / 2) /// + invlogit(- `coefs'[1, 1]) / `delta' if `tfit' == 0 } if "`disttype'" == "wei" { gen `fitdens' = (invlogit(`coefs'[1, 1]) /// * exp(- ((`tfit' * exp(`coefs'[1, 2]))^exp(`coefs'[1, 3])) /// -lngamma(1 + 1/exp(`coefs'[1, 3])))* exp(`coefs'[1, 2]) /// + invlogit(- `coefs'[1, 1]) / `delta') } la var `fitdens' "Fitted density" if "`ytitle'" == "" { local ytitle = "ytitle(Density)" } if "`reverse'" == "reverse" { replace `lowcutlev' = `tend' - `lowcutlev' replace `tfit' = `tend' - `tfit' } else { replace `tfit' = `tstart' + `tfit' } if `converged' == 1 { twoway (bar `obsdens' `lowcutlev', bartype(spanning) /// bstyle(histogram) `fcolor') /// (line `fitdens' `tfit', `lpattern'), /// `ytitle' `legend' `options' xtitle(`xtitletxt') } else { twoway (bar `obsdens' `lowcutlev', bartype(spanning) /// bstyle(histogram) `fcolor'), /// `ytitle' `legend' `options' xtitle(`xtitletxt') } if "`replace'" != "" { tempfile diagdat gen obstime = `lowcutlev' la var obstime "Observed time variable (lower limit of bins)" rename `obsdens' obsdens la var obsdens "Observed density (bin height)" gen fitdens = `fitdens' la var fitdens "Fitted density" gen fittime = `tfit' la var fittime "Time variable for fitted density" keep obstime obsdens fittime fitdens order obstime obsdens fittime fitdens sa `diagdat' restore use `diagdat', clear } else { restore } } end