<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Load the required modules
Step1:  Define a procedure for estimating the Hurst exponent.
Step2:  The procedure below loads data from the specified files and tests the pooled subcrossing size distribution against teh conjectured one.
Step3:  Load the simulated reults.
Step4:  Set the levels of the crossing tree where the self-similarity seems to kick in.
Step6:  Plot the estimated distribution
Step7:  Demonstrate that there are not that little crossings of sizes greater than 4.
Step8:  Show the estimated Hurst exponents
Step9:  Test the self similarity hypothesis.
<ASSISTANT_TASK:>
Python Code:
import os, re
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import warnings
warnings.filterwarnings("ignore")
import self_similarity as ss
from main import sim_load, list_files
Explanation: Load the required modules
End of explanation
def xing_hurst( file_name, p = 5, q = 7 ) :
    levels = np.arange( p, q + 1 ) - 1
    H, Njn, Djnk, Vjnde = sim_load( file_name )
    Cjn = np.reshape( np.sum( Djnk, axis = ( 2, ), dtype = np.float ), Djnk.shape[ :-1 ] + ( 1, ) )
    Zk = 2 * np.reshape( 1 + np.arange( Djnk.shape[ -1 ] ), ( 1, 1, Djnk.shape[ -1 ] ) )
    Pjnk = np.where( Cjn > 0, Djnk / Cjn, 0 )
    Mjn = np.reshape( np.sum( Pjnk * Zk, axis = ( 2, ) ), Cjn.shape )
    Hjn = np.log( 2 ) / np.log( Mjn )
    return ( ( 1 + np.arange( Djnk.shape[ -2 ] ) )[levels], H,
        np.average( Hjn, axis = (0, ) ).flatten( )[levels],
        np.std( Hjn, axis = (0, ) ).flatten( )[levels] )
Explanation: Define a procedure for estimating the Hurst exponent.
End of explanation
def xing_probs( file_name, p = 5, q = 7 ) :
    levels = np.arange( p, q + 1 ) - 1
    H, Njn, Djnk, Vjnde = sim_load( file_name )
    Cj = np.reshape( np.sum( Djnk[:, levels ], axis = ( 1, 2, ), dtype = np.float ), ( Djnk.shape[ 0 ], 1, ) )
    Pjk = np.where( Cj > 0, np.sum( Djnk[:, levels ], axis = ( 1, ) ) / Cj, 0 )
    Zk = np.reshape( 1 + np.arange( Pjk.shape[ -1 ] ), ( 1, Pjk.shape[ -1 ] ) )
    theta = 2.0**( 1.0 - 1.0 / H )
    Pk = theta * ( 1 - theta )**( Zk - 1 )
    return H, 2*Zk, Pk, np.average( Pjk, axis = ( 0, ) ), np.std( Pjk, axis = ( 0, ) )
Explanation: The procedure below loads data from the specified files and tests the pooled subcrossing size distribution against teh conjectured one.
End of explanation
# files = list_files( './output/25/std/', pattern = r'\.npz$' )
# files = list_files( '/Users/user/Dropbox/study_notes/year_14_15/course_project/code/output/HRM-2_19-32/med/', pattern = r'\.npz$' )
files = list_files( './output/HRM/local/HRM-4_18-16/med/', pattern = r'\.npz$' )
# files = list_files( './output/20150527/HRM-3_20-16/med/', pattern = r'\.npz$' )
Explanation: Load the simulated reults.
End of explanation
p, q = 5, 7
Explanation: Set the levels of the crossing tree where the self-similarity seems to kick in.
End of explanation
def draw_probs( files, ax, theor_color = 'gray' ) :
    ax.set_xticks( np.arange( 2, 42, 2) )
    ax.grid()
    ax.set_color_cycle( plt.cm.rainbow( np.linspace( 0, 1, max( len( files ), 10 ) )[::-1] ) )
    for f in files[1:] :
        H, Zk, Pk, Phat, Pstd = xing_probs( f, p = p, q = q )
        ax.plot( Zk[0], Pk[0], linestyle = '-', color = theor_color )
        mask = Phat > 0
        lwr = np.where( Phat > Pstd, Pstd, 0 )
        ax.errorbar( Zk[0][mask], Phat[mask], yerr = [lwr[mask], Pstd[mask]],
            fmt = '-o', label = 'HRM-4 H=%.2f' % ( H, ) )
## Create an appropriately sized figure
plt.figure( figsize = ( 16, 9 ) )
ax = plt.subplot( 111 )
## Setup view of the polt's main area
draw_probs( files, ax, theor_color = 'black' )
ax.set_yscale( 'log' ) ; ax.set_ylim( 1e-8, 1 ) ; ax.set_xlim( 2, 40 )
## Add legend
legend = ax.legend( loc = 'lower right', frameon = 1 )
frame = legend.get_frame( )
frame.set_facecolor( 'whitesmoke' )
legend.set_title( r'Empirical' )
## Add proper labels and titles
ax.set_title(rSubcrossing size distribution: empirical against conjectured $\mathbb{P}(Z=2k) = \theta \cdot (1-\theta)^{k-1}$ with $\theta = 2^{1-H^{-1}}$)
ax.set_ylabel( r'Probability (log)' )
ax.set_xlabel( r'Subcrossing size' )
## Add a zoomed view of the small size subcrossings region
ax2 = plt.axes( [.64, .54, .25, .35 ], axisbg = 'w' )
draw_probs( files, ax2, theor_color = 'gray' )
ax2.set_yscale( 'log' ) ; ax2.set_ylim( .25e-2, 1 ) ; ax2.set_xlim( 1, 7 )
## Output
plt.savefig('./output/HRM-4_18-16-log_probs.png')
plt.show()
Explanation: Plot the estimated distribution
End of explanation
# H, Njn, Djnk, Vjnde = sim_load( './output/25/std/fbm_std_20150517-075535_25_0.5000_10.npz' )
H, Njn, Djnk, Vjnde = sim_load( '/Users/user/Dropbox/study_notes/year_14_15/course_project/code/output/HRM-4_19-32/med/HRM-4_med_20150526-073757_19-32_0.5000_1000.npz' )
print "Average number of observed subcrossings of specified size per each sample path replication"
print np.average( Djnk, axis = ( 0, ) )[p:q+1, np.arange( 6, 14, 2)//2-1]
Explanation: Demonstrate that there are not that little crossings of sizes greater than 4.
End of explanation
p, q = 1, 19
plt.figure( figsize = ( 16, 9 ) )
ax = plt.subplot( 111 )
ax.set_title('Crossing tree estimates of the Hurst exponent')
ax.set_xlabel( 'Tree level' ) ; ax.set_ylabel( '$H$', rotation = 0 )
ax.set_xlim( p, q ) ; ax.set_ylim( 0.45, 1.0 )
ax.set_color_cycle( plt.cm.rainbow( np.linspace( 0, 1, max( len( files ), 10 ) )[::-1] ) )
for f in files[1:] :
    Ln, H, Hhat, Hstd = xing_hurst( f, p = p, q = q )
    ax.axhline( y = H, color = 'black', linestyle = '-.' )
    lwr = np.where( Hhat > Hstd, Hstd, 0 )
    mask = Hhat > H*0.75
    ax.errorbar( Ln[mask], Hhat[mask], yerr = [ lwr[mask], Hstd[mask] ],
        fmt = '-o', label = 'fBM H=%.2f' % ( H, ) )
## Add a legend with white opaque background.
legend = ax.legend( loc = 'lower right', frameon = 1 )
frame = legend.get_frame( )
frame.set_facecolor( 'whitesmoke' )
Explanation: Show the estimated Hurst exponents
End of explanation
p, q, h = 7, 8, 5
plt.figure( figsize = ( 16, 9 ) )
ax = plt.subplot( 111 )
prc = np.arange( 1, 100 )/100.0
ax.set_color_cycle( plt.cm.rainbow( np.linspace( 0, 1, len( files ) ) ) )
labels, data = [], []
for f in files :
    H, Njn, Djnk, Vjnde = sim_load( f, return_durations = False )
    stat, dof, pv = ss.pooled( Djnk, Njn, p = p, q = q, h = h )
    labels.append( 'fBM H=%.2f' % ( H, ) )
    data.append( pv )
    ax.hist( pv, bins = prc, label = 'fBM H=%.2f' % ( H, ), histtype = 'barstacked', stacked= True )
# ax.hist( data, bins = 10, label = labels, histtype = 'bar', stacked = False )
#     ax.plot( np.percentile( pv, prc ), prc, '-o', label = 'fBM H=%.2f' % ( H, ) )
legend = ax.legend( loc = 'upper right', frameon = 1 )
frame = legend.get_frame( )
frame.set_facecolor( 'whitesmoke' )
H, Njn, Djnk, Vjnde, Wjnp, Wbarjn, Wstdjn = sim_load(
#     './output/25/std/fbm_std_20150517-075535_25_0.5000_10.npz', return_durations = True )
    './output/HRM-2-15-4/med/HRM-2_med_20150525-144147_15-16_0.5000_100.npz', return_durations = True )
plt.figure( figsize = ( 16, 9 ) )
ax = plt.subplot( 111 )
# ax.set_yscale( 'log' )
# ax.errorbar( x = np.arange( Wbarjn.shape[1] ), y = np.average( Wjnp, axis = ( 0, ) ).T )
plt.plot( np.log( np.average( Wjnp, axis = ( 0, ) ).T ) )
ax = plt.subplot( 111 )
ax.set_color_cycle( plt.cm.rainbow( np.linspace( 0, 1, Wjnp.shape[1] ) ) )
prc = np.array( [ 0.5, 1.0, 2.5, 5.0, 10, 25, 50, 75, 90, 95, 97.5, 99, 99.5 ] )
ax.plot( np.log( np.average( Wjnp, axis = ( 0, ) )[:14].T ), prc )
# H, Njn, Djnk, Vjnde = sim_load( './output/25/std/fbm_std_20150517-080403_25_0.6500_10.npz' )
basepath = '/Users/user/Dropbox/study_notes/year_14_15/course_project/code/output/HRM-3_19-32/med/'
H, Njn, Djnk, Vjnde = sim_load( basepath+'HRM-3_med_20150526-042928_19-32_0.5000_1000.npz' )
p, q = 6, 10
levels = np.arange( p, 1 + q ) - 1
Cjnd = np.reshape( np.sum( Vjnde, axis = ( 3, ), dtype = np.float ), ( Vjnde.shape[ :-1 ] + ( 1, ) ) )
V_hat = Vjnde[:,levels] / Cjnd[:,levels]
Vpool = np.sum( Vjnde[:,levels], axis = ( 1, ) ) / np.sum( Cjnd[:,levels], axis = ( 1, ) )
print np.average( Vpool, axis = 0 )
print np.std( Vpool, axis = 0 )
theo = 2**(1-1.0/H)
# print theo
pv = 1.0 / np.sqrt( 2.0 / theo )
print pv
print np.array( np.sum( Djnk > 0, axis = ( 0, ) ) / 50, np.int )#[p:q+1,:]
Explanation: Test the self similarity hypothesis.
End of explanation
<END_TASK>