*! 2.0.3 NJC 22 July 2004 *! 2.0.2 NJC 4 July 2004 *! 2.0.1 NJC 22 June 2004 *! 2.0.0 NJC 16 June 2004 *! 1.0.0 NJC 19 June 2002 program dpplot, sort version 8 syntax varname [if] [in] [, GENerate(str) dist(str) param(numlist) /// a(real 0.5) plot(str asis) line(str asis) diff *] if "`generate'" != "" { cap confirm new var `generate' if _rc { di as err "generate() must specify new variables" exit _rc } local nvars : word count `generate' if `nvars' != 2 { di as err "generate() must specify two variables" exit 198 } } // observations to use marksample touse qui count if `touse' if r(N) == 0 error 2000 qui replace `touse' = - `touse' sort `touse' `varlist' // plotting positions (i - a ) / (n - 2a + 1) // e.g. default a = 0.5: (i - 0.5) / n // a = 0: i / (n + 1) tempvar P ft fo qui gen `P' = (_n - `a') / (`r(N)' - 2 * `a' + 1) if `touse' // TLAs (may be) allowed for distribution options local l = max(3, length("`dist'")) // case-insensitive (allowing Gaussian, Weibull, etc.) local dist = lower("`dist'") // distribution defaults to normal (Gaussian) if "`dist'" == "" local dist normal else if "`dist'" == substr("normal", 1, `l') | /// "`dist'" == substr("gaussian", 1, `l') { local dist normal } else if "`dist'" == substr("lognormal", 1, `l') { local dist lognormal } else if "`dist'" == substr("gamma", 1, `l') { local dist gamma } else if "`dist'" == substr("exponential", 1, `l') { local dist exponential } else if "`dist'" == substr("weibull", 1, `l') { local dist weibull } // any other name falls through: will be treated as call to dp_`dist' // dp_`dist' calculates theoretical and "observed" densities qui dp_`dist' `ft' `fo' `P' `varlist' `touse' `param' _crcslbl `fo' `varlist' // graph defaults // caption set by called program local ytitle "Probability density" if "`diff'" != "" { tempvar d gen `d' = `fo' - `ft' label var `d' "difference" } // graph twoway mspline `ft' `varlist' if `touse', bands(200) `line' || /// scatter `fo' `d' `varlist' if `touse', /// yti("`ytitle'") caption("`caption'", size(medsmall)) /// yla(, ang(h)) legend(order(2 1 "`dist'") off) `options' || /// `plot' // messages about missing values will be visible if "`generate'" != "" { tokenize `generate' gen `1' = `ft' if `touse' label var `1' "`dist' density, `varlist' (direct)" gen `2' = `fo' if `touse' label var `2' "`dist' density, `varlist' (indirect)" } end program dp_laplace args ft fo P varlist touse theta sigma garbage if "`garbage'" != "" { di as err "too many parameters specified" exit 198 } if "`theta'" != "" & "`sigma'" == "" { di as err "need to specify two parameters" exit 198 } qui if "`theta'" == "" { su `varlist' if `touse', d local theta = r(p50) tempvar dev gen `dev' = abs(`varlist' - r(p50)) if `touse' su `dev', meanonly local sigma = r(mean) } gen double `ft' = /// (1 / (2 * `sigma')) * exp(-abs(`varlist' - `theta') / `sigma') if `touse' gen double `fo' = `theta' + /// `sigma' * cond(`P' <= 0.5, log(2 * `P'), -log(2 * (1 - `P'))) if `touse' replace `fo' = /// (1 / (2 * `sigma')) * exp(-abs(`fo' - `theta') / `sigma') if `touse' local theta = trim("`: di %7.0g `theta''") local sigma = trim("`: di %7.0g `sigma''") c_local caption : di "reference Laplace, theta `theta' sigma `sigma'" end program dp_alaplace args ft fo P varlist touse theta sigma kappa garbage if "`garbage'" != "" { di as err "too many parameters specified" exit 198 } if ("`theta'" == "") | ("`sigma'" == "" | "`kappa'" == "") { di as err "need to specify three parameters" exit 198 } tempname K1 K2 scalar `K1' = `kappa' / (1 + (`kappa')^2) scalar `K2' = ((`kappa')^2) / (1 + (`kappa')^2) local v `varlist' gen double `ft' = (`K1' / `sigma') * exp((`v' - `theta') * /// cond(`v' >= `theta', -`kappa' / `sigma', 1 / (`sigma' * `kappa'))) /// if `touse' gen double `fo' = `theta' + (`sigma' * `kappa' / sqrt(2)) * /// (ln(`P') - ln(`K2')) if `touse' & `P' <= `K2' replace `fo' = `theta' - (`sigma' / (`kappa' * sqrt(2))) * /// ln((1 - `P') * (1 + (`kappa')^2)) if `touse' & `P' > `K2' replace `fo' = (`K1' / `sigma') * exp((`fo' - `theta') * /// cond(`fo' >= `theta', -`kappa' / `sigma', 1 / (`sigma' * `kappa'))) /// local theta = trim("`: di %7.0g `theta''") local sigma = trim("`: di %7.0g `sigma''") local kappa = trim("`: di %7.0g `kappa''") c_local caption : di "reference asymmetric Laplace, theta `theta' sigma `sigma' kappa `kappa'" end program dp_normal args ft fo P varlist touse mean sd garbage if "`garbage'" != "" { di as err "too many parameters specified" exit 198 } if "`mean'" != "" & "`sd'" == "" { di as err "need to specify two parameters" exit 198 } if "`mean'" == "" { su `varlist' if `touse' local mean = `r(mean)' local sd = `r(sd)' } gen double `ft' = (1 / `sd') * normden((`varlist' - `mean') / `sd') if `touse' gen double `fo' = (1 / `sd') * normden(invnorm(`P')) if `touse' local mean = trim("`: di %7.0g `mean''") local sd = trim("`: di %7.0g `sd''") c_local caption : di "reference normal, mean `mean' sd `sd'" end program dp_beta args ft fo P varlist touse alpha beta garbage if "`garbage'" != "" { di as err "too many parameters specified" exit 198 } if "`alpha'" == "" { di as err "need to specify parameters" exit 198 } if "`alpha'" != "" & "`beta'" == "" { di as err "need to specify two parameters" exit 198 } tempvar x gen double `ft' = betaden(`alpha', `beta', `varlist') if `touse' gen double `x' = invibeta(`alpha',`beta',`P') if `touse' gen double `fo' = betaden(`alpha', `beta', `x') local alpha = trim("`: di %7.0g `alpha''") local beta = trim("`: di %7.0g `beta''") c_local caption : di "reference beta, alpha `alpha' beta `beta'" end program dp_gamma args ft fo P varlist touse alpha beta garbage if "`garbage'" != "" { di as err "too many parameters specified" exit 198 } if "`alpha'" == "" { di as err "need to specify parameters" exit 198 } if "`alpha'" != "" & "`beta'" == "" { di as err "need to specify two parameters" exit 198 } tempvar x gen double `ft' = gammaden(`alpha', `beta', 0, `varlist') if `touse' gen double `x' = `beta' * invgammap(`alpha',`P') if `touse' gen double `fo' = gammaden(`alpha', `beta', 0, `x') if `touse' local alpha = trim("`: di %7.0g `alpha''") local beta = trim("`: di %7.0g `beta''") c_local caption : di "reference gamma, alpha `alpha' beta `beta'" end program dp_gumbel args ft fo P varlist touse alpha mu garbage if "`garbage'" != "" { di as err "too many parameters specified" exit 198 } if "`alpha'" == "" { di as err "need to specify parameters" exit 198 } if "`alpha'" != "" & "`mu'" == "" { di as err "need to specify two parameters" exit 198 } tempvar x gen double `ft' = (1 / `alpha') * exp(-(`varlist' - `mu') / `alpha') * /// exp(-exp(-(`varlist' - `mu') / `alpha')) if `touse' gen double `x' = `mu' - `alpha' * log(-log(`P')) if `touse' gen double `fo' = (1 / `alpha') * exp(-(`x' - `mu') / `alpha') * /// exp(-exp(-(`x' - `mu') / `alpha')) local alpha = trim("`: di %7.0g `alpha''") local mu = trim("`: di %7.0g `mu''") c_local caption : di "reference Gumbel, alpha `alpha' mu `mu'" end program dp_kappa2 args ft fo P varlist touse alpha beta garbage if "`garbage'" != "" { di as err "too many parameters specified" exit 198 } if "`alpha'" == "" { di as err "need to specify parameters" exit 198 } if "`alpha'" != "" & "`beta'" == "" { di as err "need to specify two parameters" exit 198 } tempvar x gen double `ft' = (`alpha' / `beta') * /// (`alpha' + (`varlist' / `beta')^(`alpha'))^(-(`alpha' + 1) / `alpha') if `touse' gen double `x' = `beta' * `P' * /// ((`alpha')/(1 - `P'^`alpha'))^(1 / `alpha') if `touse' gen double `fo' = (`alpha' / `beta') * /// (`alpha' + (`x' / `beta')^(`alpha'))^(-(`alpha' + 1) / `alpha') local alpha = trim("`: di %7.0g `alpha''") local beta = trim("`: di %7.0g `beta''") c_local caption : di "reference kappa, alpha `alpha' beta `beta'" end program dp_dagum args ft fo P varlist touse a b p garbage if "`garbage'" != "" { di as err "too many parameters specified" exit 198 } if "`a'" == "" { di as err "need to specify parameters" exit 198 } if "`a'" != "" & ("`b'" == "" | "`p'" == "") { di as err "need to specify three parameters" exit 198 } tempvar x gen double `ft' = 1 + (`b' / `varlist')^`a' if `touse' replace `ft' = /// `a' * `p' * (`b' / `varlist')^`a' * (1 / `varlist') / `ft'^(`p' + 1) gen double `x' = `b' * ((`P')^(-1 / `p') - 1)^(-1 / `a') if `touse' gen double `fo' = 1 + (`b' / `x')^`a' replace `fo' = /// `a' * `p' * (`b' / `x')^`a' * (1 / `x') / `fo'^(`p' + 1) local a = trim("`: di %7.0g `a''") local b = trim("`: di %7.0g `b''") local p = trim("`: di %7.0g `p''") c_local caption : di "reference Dagum, a `a' b `b' p `p'" end program dp_sm args ft fo P varlist touse a b q garbage if "`garbage'" != "" { di as err "too many parameters specified" exit 198 } if "`a'" == "" { di as err "need to specify parameters" exit 198 } if "`a'" != "" & ("`b'" == "" | "`q'" == "") { di as err "need to specify three parameters" exit 198 } tempvar x gen double `ft' = 1 + (`varlist' / `b')^`a' if `touse' replace `ft' = /// `a' * `q' / `b' * `ft'^(-(`q' + 1)) * (`varlist' / `b')^(`a' - 1) gen double `x' = /// `b' * (( 1 / (1 - `P'))^(1 / `q') - 1)^(1 / `a') if `touse' gen double `fo' = 1 + (`x' / `b')^`a' replace `fo' = /// `a' * `q' / `b' * `fo'^(-(`q' + 1)) * (`x' / `b')^(`a' - 1) local a = trim("`: di %7.0g `a''") local b = trim("`: di %7.0g `b''") local q = trim("`: di %7.0g `q''") c_local caption : di "reference Singh-Maddala, a `a' b `b' q `q'" end program dp_gb2 args ft fo P varlist touse a b p q garbage if "`garbage'" != "" { di as err "too many parameters specified" exit 198 } if "`a'" == "" { di as err "need to specify parameters" exit 198 } if "`a'" != "" & ("`b'" == "" | "`p'" == "" | "`q'" == "") { di as err "need to specify four parameters" exit 198 } tempname factor tempvar v1 v2 x scalar `factor' = lngamma(`p') + lngamma(`q') - lngamma(`p' + `q') gen double `v1' = (`a' * `p' - 1) * log(`varlist') if `touse' gen double `v2' = (`a' * `p' * log(`b')) + `factor' + /// ((`p' + `q') * log(1 + (`varlist' / `b')^`a')) gen double `ft' = exp(log(`a') + `v1' - `v2') gen double `x' = invibeta(`p',`q',`P') if `touse' replace `x' = `b' * (`x' / (1 - `x'))^(1 / `a') replace `v1' = (`a' * `p' - 1) * log(`x') replace `v2' = (`a' * `p' * log(`b')) + `factor' + /// ((`p' + `q') * log(1 + (`x' / `b')^`a')) gen double `fo' = exp(log(`a') + `v1' - `v2') local a = trim("`: di %7.0g `a''") local b = trim("`: di %7.0g `b''") local p = trim("`: di %7.0g `p''") local q = trim("`: di %7.0g `q''") c_local caption : di "reference Generalized Beta, a `a' b `b' p `p' q `q'" end program dp_lognormal args ft fo P varlist touse mean sd garbage if "`garbage'" != "" { di as err "too many parameters specified" exit 198 } if "`mean'" != "" & "`sd'" == "" { di as err "need to specify two parameters" exit 198 } qui count if `varlist' <= 0 & `touse' if r(N) > 0 error 411 tempvar logvar logvar2 gen double `logvar' = log(`varlist') if `touse' if "`mean'" == "" { su `logvar' if `touse' local mean = `r(mean)' local sd = `r(sd)' } gen double `ft' = /// (1 / `varlist' * `sd' * sqrt(2 * _pi)) * /// exp(-(`logvar' - `mean')^2 / (2 * `sd'^2)) if `touse' gen double `fo' = exp(invnorm(`P') * `sd' + `mean') if `touse' gen double `logvar2' = log(`fo') replace `fo' = (1 / `fo' * `sd' * sqrt(2 * _pi)) * /// exp(-(`logvar2' - `mean')^2 / (2 * `sd'^2)) local mean = trim("`: di %7.0g `mean''") local sd = trim("`: di %7.0g `sd''") c_local caption : /// di "reference lognormal, mean of logs `mean' sd of logs `sd'" end program dp_weibull args ft fo P varlist touse b c garbage if "`garbage'" != "" { di as err "too many parameters specified" exit 198 } if "`b'" != "" & "`c'" == "" { di as err "need to specify two parameters" exit 198 } qui count if `varlist' < 0 & `touse' if r(N) > 0 { di as err "negative values encountered" exit 411 } tempname eb weibull `varlist' if `touse' nlcom (b : exp(_b[_cons])) (c : exp([ln_p][_cons])), post mat `eb' = e(b) local b = `eb'[1,1] local c = `eb'[1,2] gen double `ft' = (`c' / `b') * (`varlist' / `b')^(`c' - 1) * /// exp(-(`varlist'/ `b')^`c') if `touse' gen double `fo' = `b' * (-ln(1 - `P'))^(1 / `c') if `touse' replace `fo' = (`c' / `b') * (`fo' / `b')^(`c' - 1) * /// exp(-(`fo'/ `b')^`c') local b = trim("`: di %7.0g `b''") local c = trim("`: di %7.0g `c''") c_local caption : di "reference Weibull, b `b' c `c'" end program dp_exponential args ft fo P varlist touse mean garbage if "`garbage'" != "" { di as err "too many parameters specified" exit 198 } count if `varlist' < 0 & `touse' if r(N) > 0 { di as err "negative values encountered" exit 411 } if "`mean'" == "" { su `varlist' if `touse', meanonly local mean = `r(mean)' } gen double `ft' = (1 / `mean') * exp(-(`varlist' / `mean')) if `touse' gen double `fo' = (1 - `P') / `mean' if `touse' local mean = trim("`: di %7.0g `mean''") c_local caption : di "reference exponential, mean `mean'" end