<s> from __future__ import print_function <EOL> import lmfit <EOL> import numpy as np <EOL> import emcee <EOL> try : <EOL> import matplotlib . pyplot as plt <EOL> HASPYLAB = True <EOL> except ImportError : <EOL> HASPYLAB = False <EOL> def create_prior ( params ) : <EOL> """<STR_LIT>""" <EOL> none_to_inf = lambda x , sign = <NUM_LIT:1> : sign * np . inf if x is None else x <EOL> lower_bounds = np . array ( [ none_to_inf ( i . min , - <NUM_LIT:1> ) for i in params . values ( ) if i . vary ] ) <EOL> upper_bounds = np . array ( [ none_to_inf ( i . max , <NUM_LIT:1> ) for i in params . values ( ) if i . vary ] ) <EOL> def bounds_prior ( values ) : <EOL> values = np . asarray ( values ) <EOL> is_ok = np . all ( ( lower_bounds < values ) & ( values < upper_bounds ) ) <EOL> return <NUM_LIT:0> if is_ok else - np . inf <EOL> return bounds_prior <EOL> def create_lnliklihood ( mini , sigma = None ) : <EOL> """<STR_LIT>""" <EOL> def lnprob ( vals , sigma = sigma ) : <EOL> for v , p in zip ( vals , [ p for p in mini . params . values ( ) if p . vary ] ) : <EOL> p . value = v <EOL> residuals = mini . userfcn ( mini . params ) <EOL> if not sigma : <EOL> sigma = vals [ - <NUM_LIT:1> ] <EOL> val = - <NUM_LIT:0.5> * np . sum ( np . log ( <NUM_LIT:2> * np . pi * sigma ** <NUM_LIT:2> ) + ( residuals / sigma ) ** <NUM_LIT:2> ) <EOL> return val <EOL> return lnprob <EOL> def starting_guess ( mini , estimate_sigma = True ) : <EOL> """<STR_LIT>""" <EOL> vals = [ i . value for i in mini . params . values ( ) if i . vary ] <EOL> if estimate_sigma : <EOL> vals . append ( mini . residual . std ( ) ) <EOL> return vals <EOL> def create_all ( mini , sigma = None ) : <EOL> """<STR_LIT>""" <EOL> sigma_given = not sigma is None <EOL> lnprior = create_prior ( mini . params ) <EOL> lnprob = create_lnliklihood ( mini , sigma = sigma ) <EOL> guess = starting_guess ( mini , not sigma_given ) <EOL> if sigma_given : <EOL> func = lambda x : lnprior ( x [ : ] ) + lnprob ( x ) <EOL> else : <EOL> func = lambda x : lnprior ( x [ : - <NUM_LIT:1> ] ) + lnprob ( x ) <EOL> return func , guess <EOL> params = lmfit . Parameters ( ) <EOL> params . add_many ( ( '<STR_LIT:a>' , <NUM_LIT:5> ) , <EOL> ( '<STR_LIT:b>' , - <NUM_LIT:5> ) , <EOL> ( '<STR_LIT>' , <NUM_LIT:1> , <NUM_LIT:1> , <NUM_LIT:0> ) , <EOL> ( '<STR_LIT>' , <NUM_LIT:15> , <NUM_LIT:1> , <NUM_LIT:0> ) ) <EOL> x = np . linspace ( <NUM_LIT:0> , <NUM_LIT:20> , <NUM_LIT> ) <EOL> a , b , t1 , t2 = <NUM_LIT:2> , <NUM_LIT:3> , <NUM_LIT:2> , <NUM_LIT:10> <EOL> y_true = a * np . exp ( - x / t1 ) + b * np . exp ( - x / t2 ) <EOL> sigma = <NUM_LIT> <EOL> y = y_true + np . random . randn ( x . size ) * sigma <EOL> def residuals ( paras ) : <EOL> a = paras [ '<STR_LIT:a>' ] . value <EOL> b = paras [ '<STR_LIT:b>' ] . value <EOL> t1 = paras [ '<STR_LIT>' ] . value <EOL> t2 = paras [ '<STR_LIT>' ] . value <EOL> return a * np . exp ( - x / t1 ) + b * np . exp ( - x / t2 ) - y <EOL> mini = lmfit . Minimizer ( residuals , params ) <EOL> mini . leastsq ( ) <EOL> lmfit . report_errors ( params ) <EOL> lnfunc , guess = create_all ( mini ) <EOL> nwalkers , ndim = <NUM_LIT:30> , len ( guess ) <EOL> p0 = emcee . utils . sample_ball ( guess , <NUM_LIT:0.1> * np . array ( guess ) , nwalkers ) <EOL> sampler = emcee . EnsembleSampler ( nwalkers , ndim , lnfunc ) <EOL> steps = <NUM_LIT> <EOL> sampler . run_mcmc ( p0 , steps ) <EOL> if HASPYLAB : <EOL> fig , axes = plt . subplots ( <NUM_LIT:5> , <NUM_LIT:1> , sharex = True , figsize = ( <NUM_LIT:8> , <NUM_LIT:9> ) ) <EOL> for ( i , name , rv ) in zip ( range ( <NUM_LIT:5> ) , params . keys ( ) + [ '<STR_LIT>' ] , [ a , b , t1 , t2 , sigma ] ) : <EOL> axes [ i ] . plot ( sampler . chain [ : , : , i ] . T , color = "<STR_LIT:k>" , alpha = <NUM_LIT> ) <EOL> axes [ i ] . yaxis . set_major_locator ( plt . MaxNLocator ( <NUM_LIT:5> ) ) <EOL> axes [ i ] . axhline ( rv , color = "<STR_LIT>" , lw = <NUM_LIT:2> ) <EOL> axes [ i ] . set_ylabel ( "<STR_LIT>" % name ) <EOL> axes [ - <NUM_LIT:1> ] . set_xlabel ( "<STR_LIT>" ) <EOL> plt . figure ( ) <EOL> try : <EOL> import triangle <EOL> burnin = <NUM_LIT:100> <EOL> samples = sampler . chain [ : , burnin : , : ] . reshape ( ( - <NUM_LIT:1> , ndim ) ) <EOL> triangle . corner ( samples , labels = params . keys ( ) + [ '<STR_LIT>' ] , <EOL> truths = [ a , b , t1 , t2 , sigma ] ) <EOL> except ImportError : <EOL> print ( "<STR_LIT>" ) </s>
