Commit 167e55db by Olaf

### Error estimation

parent 109c44f7
errorest.py 0 → 100644
 # import import numpy as np from scipy import optimize def chi(A,x): """ function to determine chi_A(t) last modified: 20/03/19 Parameters ---------- A : (1,n_max)-array input data x : (1,n_max-1)-array x data Returns ------- chi : (1,x)-array autocorrelation function """ # Process input n_max = np.size(A) A_tile = np.tile(A,[n_max,1]) # (n_max,n_max)-array # Calculate scalefactor N N = 1/(n_max-x) # (1,x)-array # Calculate chi # # chi_1(t) = 1/(n_max-t) * sum_{n=0}^{n_max-t} A(n)A(n+t) AA = np.convolve(A,np.flip(A,0))[0:n_max][::-1] chi_1 = np.multiply(N,AA) # (1,n_max)-array # # chi_2(t) = 1/(n_max-t) * sum_{n=0}^{n_max-t} A(n) A_triu_2 = np.triu(np.flip(A_tile,1)) A_2 = np.sum(A_triu_2,1) chi_2 = np.multiply(N,A_2) # (1,n_max)-array # # chi_3(t) = 1/(n_max-t) * sum_{n=0}^{n_max-t} A(n+1) A_triu_3 = np.triu(A_tile) A_3 = np.sum(A_triu_3,1) chi_3 = np.multiply(N,A_3) # (1,n_max)-array # # chi = chi_1 + chi_2*chi_3 chi = chi_1 - np.multiply(chi_2,chi_3) # (1,n_max)-array # return return chi def test_func(x, a, b): """ exponential function to fit chi last modified: 26/03/2019 Parameters ---------- x : (1,N) array x_data a : (float) amplitude b : (float) tau Returns ------- a*np.exp(-x/b) : (1,N) array function """ return a*np.exp(-x/b) def tau(y_data): """ function to find correlation time of the simulation (tau) by fitting chi to exponential test function last modified: 26/03/2019 Parameters ---------- y_data : (1,N) array data to which function is fitted Returns ------- tau : (float) correlation time of the simulation """ # data x_data = np.array(range(len(y_data))) # fit and determine tau [amplitude, tau], params_covariance = optimize.curve_fit(test_func, x_data,y_data) # return return tau def calc_std_mean(tau,data): """ function to calculate standard deviation (std) and mean (mean) of data using data blocking last modified: 20/03/19 Parameters ---------- tau : float correlation time of the simulation data : (1,N)-array data from which std and mean are calculated Returns ------- data_std : float standard deviation of data data_mean : float mean of data """ # define length of data block m = 2*np.abs(np.ceil(tau)) # select multiple of m datasets of data data_select = data[:int(len(data)-len(data)%m)] # reshape data to (n,m)-array n = len(data_select)/m data_reshape= data_select.reshape(int(n),int(m)) # average in data_blocks data_j = np.mean(data_reshape,1) # determine std and mean from average data_blocks data_std = np.std(data_j) data_mean= np.mean(data_j) # return return data_std, data_mean def error_combine(data): """ overhead function to calculate standard deviation and mean using datablocking last modified: 20/03/2019 Parameters ---------- data : (1,N)-array data from which std and mean are calculated Returns ------- data_std : (float) standard deviation of data data_mean : (float) mean of data """ # process input n_max = np.size(data) x_data= np.linspace(0,n_max-1,n_max) # calculate chi chi_output = chi(data,x_data) # calculate tau tau_output = tau(chi_output) # calculate std and mean [data_std, data_mean] = calc_std_mean(tau_output,data) # return return data_std, data_mean \ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment