<s> """<STR_LIT>""" <EOL> __all__ = [ "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , "<STR_LIT>" , <EOL> "<STR_LIT>" , "<STR_LIT>" ] <EOL> import re <EOL> import numpy as np <EOL> from . formatnum import format_uncertainty <EOL> class VarStats ( object ) : <EOL> def __init__ ( self , ** kw ) : <EOL> self . __dict__ = kw <EOL> def var_stats ( draw , vars = None ) : <EOL> if vars is None : <EOL> vars = range ( draw . points . shape [ <NUM_LIT:1> ] ) <EOL> return [ _var_stats_one ( draw , v ) for v in vars ] <EOL> ONE_SIGMA = <NUM_LIT:1> - <NUM_LIT:2> * <NUM_LIT> <EOL> def _var_stats_one ( draw , var ) : <EOL> weights , values = draw . weights , draw . points [ : , var ] . flatten ( ) <EOL> best_idx = np . argmax ( draw . logp ) <EOL> best = values [ best_idx ] <EOL> p95 , p68 , p0 = credible_intervals ( x = values , weights = weights , <EOL> ci = [ <NUM_LIT> , ONE_SIGMA , <NUM_LIT:0.0> ] ) <EOL> mean , std = stats ( x = values , weights = weights ) <EOL> vstats = VarStats ( label = draw . labels [ var ] , index = var + <NUM_LIT:1> , <EOL> p95 = p95 , p68 = p68 , <EOL> median = p0 [ <NUM_LIT:0> ] , mean = mean , std = std , best = best ) <EOL> return vstats <EOL> def format_num ( x , place ) : <EOL> precision = <NUM_LIT:10> ** place <EOL> digits_after_decimal = abs ( place ) if place < <NUM_LIT:0> else <NUM_LIT:0> <EOL> return "<STR_LIT>" % ( digits_after_decimal , np . round ( x / precision ) * precision ) <EOL> def format_vars ( all_vstats ) : <EOL> v = dict ( parameter = "<STR_LIT>" , <EOL> mean = "<STR_LIT>" , median = "<STR_LIT>" , best = "<STR_LIT>" , <EOL> interval68 = "<STR_LIT>" , <EOL> interval95 = "<STR_LIT>" ) <EOL> s = [ "<STR_LIT>" <EOL> "<STR_LIT>" % v ] <EOL> for v in all_vstats : <EOL> place = int ( np . log10 ( v . p95 [ <NUM_LIT:1> ] - v . p95 [ <NUM_LIT:0> ] ) ) - <NUM_LIT:2> <EOL> summary = dict ( mean = format_uncertainty ( v . mean , v . std ) , <EOL> median = format_num ( v . median , place - <NUM_LIT:1> ) , <EOL> best = format_num ( v . best , place - <NUM_LIT:1> ) , <EOL> lo68 = format_num ( v . p68 [ <NUM_LIT:0> ] , place ) , <EOL> hi68 = format_num ( v . p68 [ <NUM_LIT:1> ] , place ) , <EOL> loci = format_num ( v . p95 [ <NUM_LIT:0> ] , place ) , <EOL> hici = format_num ( v . p95 [ <NUM_LIT:1> ] , place ) , <EOL> parameter = v . label , <EOL> index = v . index ) <EOL> s . append ( "<STR_LIT>" <EOL> "<STR_LIT>" % summary ) <EOL> return "<STR_LIT:\n>" . join ( s ) <EOL> VAR_PATTERN = re . compile ( r"""<STR_LIT>""" , re . VERBOSE ) <EOL> def parse_var ( line ) : <EOL> """<STR_LIT>""" <EOL> m = VAR_PATTERN . match ( line ) <EOL> if m : <EOL> exp = int ( m . group ( '<STR_LIT>' ) ) if m . group ( '<STR_LIT>' ) else <NUM_LIT:0> <EOL> return VarStats ( index = int ( m . group ( '<STR_LIT>' ) ) , <EOL> name = m . group ( '<STR_LIT>' ) , <EOL> mean = float ( m . group ( '<STR_LIT>' ) ) * <NUM_LIT:10> ** exp , <EOL> median = float ( m . group ( '<STR_LIT>' ) ) , <EOL> best = float ( m . group ( '<STR_LIT>' ) ) , <EOL> p68 = ( float ( m . group ( '<STR_LIT>' ) ) , float ( m . group ( '<STR_LIT>' ) ) ) , <EOL> p95 = ( float ( m . group ( '<STR_LIT>' ) ) , float ( m . group ( '<STR_LIT>' ) ) ) , <EOL> ) <EOL> else : <EOL> return None <EOL> def stats ( x , weights = None ) : <EOL> """<STR_LIT>""" <EOL> if weights is None : <EOL> x = np . sort ( x ) <EOL> mean , std = np . mean ( x ) , np . std ( x , ddof = <NUM_LIT:1> ) <EOL> else : <EOL> mean = np . mean ( x * weights ) / np . sum ( weights ) <EOL> var = np . sum ( ( weights * ( x - mean ) ) ** <NUM_LIT:2> ) / np . sum ( weights ) <EOL> std = np . sqrt ( var ) <EOL> return mean , std <EOL> def credible_intervals ( x , ci , weights = None ) : <EOL> """<STR_LIT>""" <EOL> from numpy import asarray , vstack , sort , cumsum , searchsorted , round , clip <EOL> ci = asarray ( ci , '<STR_LIT:d>' ) <EOL> target = ( <NUM_LIT:1> + vstack ( ( - ci , + ci ) ) ) . T / <NUM_LIT:2> <EOL> if weights is None : <EOL> idx = clip ( round ( target * ( x . size - <NUM_LIT:1> ) ) , <NUM_LIT:0> , x . size - <NUM_LIT:1> ) . astype ( '<STR_LIT:i>' ) <EOL> return sort ( x ) [ idx ] <EOL> else : <EOL> idx = np . argsort ( x ) <EOL> x , weights = x [ idx ] , weights [ idx ] <EOL> w = cumsum ( weights / sum ( weights ) ) <EOL> return x [ searchsorted ( w , target ) ] <EOL> def shortest_credible_interval ( x , ci = <NUM_LIT> , weights = None ) : <EOL> """<STR_LIT>""" <EOL> sorted = np . all ( x [ <NUM_LIT:1> : ] >= x [ : - <NUM_LIT:1> ] ) <EOL> if not sorted : <EOL> idx = np . argsort ( x ) <EOL> x = x [ idx ] <EOL> if weights is not None : <EOL> weights = weights [ idx ] <EOL> if weights is not None : <EOL> w = np . cumsum ( weights / sum ( weights ) ) <EOL> idx = np . searchsorted ( w , np . arange ( <NUM_LIT:0> , <NUM_LIT:1> , <NUM_LIT> ) ) <EOL> x = x [ idx ] <EOL> if np . isscalar ( ci ) : <EOL> return _find_interval ( x , ci ) <EOL> else : <EOL> return [ _find_interval ( x , i ) for i in ci ] <EOL> def _find_interval ( x , ci ) : <EOL> """<STR_LIT>""" <EOL> n = len ( x ) <EOL> size = int ( ci * n + np . sqrt ( <NUM_LIT:1> - ci ) * np . log ( n ) ) <EOL> if size >= n : <EOL> return x [ <NUM_LIT:0> ] , x [ - <NUM_LIT:1> ] <EOL> else : <EOL> width = x [ size : ] - x [ : - size ] <EOL> idx = np . argmin ( width ) <EOL> return x [ idx ] , x [ idx + size ] </s>
