*! version 2.0.1 September 23, 2014 David Roodman *! version 1.0.1 November 2, 2007 Scott Merryman *! Based on on -betafit- by Cox, Jenkins, and Buis * rewritten by David Roodman 9/8/14 cap program drop gevfit_lf2 program gevfit_lf2 version 11.0 args todo b lng lng_sig lng_xi lng_mu H tempvar x sig xi mu mleval `sig'=`b', eq(1) $gev_scale_scalar if `sig'<=0 exit // works well enough if sig is a variable mleval `xi' =`b', eq(2) $gev_shape_scalar mleval `mu' =`b', eq(3) $gev_loc_scalar qui { gen double `x' = ($S_MLy - `mu') / `sig' if $ML_samp tempvar xi_x one_plus_xi_x ln_one_plus_xi_x f inv_xi absxi gumbel gumbel2 if "$gev_shape_scalar"!="" { scalar `inv_xi' = 1/`xi' scalar `gumbel' = abs(`xi')<1e-10 scalar `gumbel2' = abs(`xi')<1e-5 } else { gen `absxi' = abs(`xi') if $ML_samp gen byte `gumbel' = `absxi'<1e-10 gen byte `gumbel2' = `absxi'<1e-5 gen double `inv_xi' = 1/`xi' if $ML_samp } gen double `xi_x' = `xi' * `x' if $ML_samp gen double `one_plus_xi_x' = 1 + `xi_x' if $ML_samp gen double `ln_one_plus_xi_x' = ln(`one_plus_xi_x') if $ML_samp gen double `f' = cond(`gumbel', exp(-`x'), (`one_plus_xi_x')^-`inv_xi') if $ML_samp replace `lng' = -ln(`sig') - cond(`gumbel', `x'+`f', (`inv_xi'+1)*`ln_one_plus_xi_x' + `f') if $ML_samp if `todo' { tempvar lng_x neg_inv_sig inv_sig2 f_x f_xi x2 one_plus_xi gen double `x2' = `x'*`x' if $ML_samp if "$gev_shape_scalar"!="" scalar `one_plus_xi' = 1+`xi' else gen double `one_plus_xi' = 1+`xi' if $ML_samp if "$gev_scale_scalar"!="" scalar `neg_inv_sig' = -1/`sig' else gen double `neg_inv_sig' = -1/`sig' if $ML_samp gen double `f_x' = -1/`one_plus_xi_x' if $ML_samp // actually f_x/f gen double `f_xi' = `inv_xi'*(-`x'/`one_plus_xi_x'+`inv_xi'*`ln_one_plus_xi_x') if $ML_samp // actually f_xi/f gen double `lng_x' = cond(`gumbel', `f'-1 , `f_x' *(`one_plus_xi'-`f') ) if $ML_samp replace `lng_xi' = cond(`gumbel', `x2'*.5*(1-`f')-`x', `f_xi'*(`one_plus_xi'-`f')+ln(`f')) if $ML_samp replace `lng_mu' = `neg_inv_sig' * `lng_x' if $ML_samp replace `lng_sig' = `neg_inv_sig' + `x' * `lng_mu' if $ML_samp if `todo'==2 { tempvar f_xx f_xxi f_xixi lng_xx lng_xxi lng_xixi lng_musig lng_muxi f_xi2 one_plus_xi_x2 tempname d2mumu d2musig d2muxi d2sigsig d2sigxi d2xixi gen double `one_plus_xi_x2' = `one_plus_xi_x'*`one_plus_xi_x' if $ML_samp gen double `f_xi2' = `f_xi'*`f_xi' if $ML_samp gen double `f_xx' = `one_plus_xi'/`one_plus_xi_x2' if $ML_samp // actually f_xx/f gen double `f_xxi' = -`f_x' * ((`inv_xi'+1)*`x'/`one_plus_xi_x' - `ln_one_plus_xi_x'/(`xi'*`xi')) if $ML_samp // actually f_xxi/f gen double `f_xixi' = ((`inv_xi'*`x2'/`one_plus_xi_x2'-2/`xi'*`f_xi')+`f_xi2') if $ML_samp // actually f_xixi/f gen double `lng_xx' = cond(`gumbel2', -`f' , `one_plus_xi'*(`f_xx' -`f_x' *`f_x' ) -`f_xx' *`f') if $ML_samp gen double `lng_xxi' = cond(`gumbel2', (`x2'*.5-`x')*`f'+`x'-1 , `one_plus_xi'*(`f_xxi' -`f_x' *`f_xi') +`f_x' -`f_xxi' *`f') if $ML_samp gen double `lng_xixi' = cond(`gumbel2', `x2'*(`f'*((2/3)*`x'-`x2'*.25)-(2/3)*`x'+1), `one_plus_xi'*(`f_xixi'-`f_xi2' )+2*`f_xi'-`f_xixi'*`f') if $ML_samp if "$gev_scale_scalar"!="" { gen double `lng_musig' = `lng_x' + `x'*`lng_xx' if $ML_samp local lng_muxi `lng_xxi' // lng_muxi needs to be multiplied by -1/sig, but do that after summing for efficiency mlmatsum `lng' `d2mumu' = `lng_xx', eq(3) mlmatsum `lng' `d2musig' = `lng_musig', eq(1,3) mlmatsum `lng' `d2muxi' = `lng_muxi', eq(2,3) mlmatsum `lng' `d2sigsig' = (1/`sig') * `x' * `lng_musig' - `lng_sig', eq(1) mlmatsum `lng' `d2sigxi' = `x' * `lng_muxi', eq(1,2) mat `d2muxi' = -`d2muxi'/`sig' mat `d2sigxi' = -`d2sigxi'/`sig' mat `d2musig' = `d2musig'/`sig'^2 mat `d2sigsig' = `d2sigsig'/`sig' mat `d2mumu' = `d2mumu'/`sig'^2 } else { gen double `inv_sig2' = `neg_inv_sig' * `neg_inv_sig' if $ML_samp gen double `lng_musig' = `inv_sig2' * (`lng_x' + `x'*`lng_xx') if $ML_samp gen double `lng_muxi' = `neg_inv_sig' * `lng_xxi' if $ML_samp mlmatsum `lng' `d2mumu' = `inv_sig2' * `lng_xx', eq(3) mlmatsum `lng' `d2musig' = `lng_musig', eq(1,3) mlmatsum `lng' `d2muxi' = `lng_muxi', eq(2,3) mlmatsum `lng' `d2sigsig' = `x' * `lng_musig' + `neg_inv_sig' * `lng_sig', eq(1) mlmatsum `lng' `d2sigxi' = `x' * `lng_muxi', eq(1,2) } mlmatsum `lng' `d2xixi' = `lng_xixi', eq(2) mat `H' = `d2sigsig',`d2sigxi',`d2musig' \ `d2sigxi'',`d2xixi',`d2muxi' \ `d2musig'',`d2muxi'', `d2mumu' } } } end