<jupyter_start><jupyter_text>Summary: Statistical effects on Edelweiss C (corrected function version)This notebook examines the error on the paramters C and m in the case where nothing else is allowed to vary.<jupyter_code>import numpy as np
import pandas as pd
from joblib import Memory
import matplotlib.pyplot as plt
%matplotlib inline 
import h5py
from scipy import optimize
import lmfit as lmf
import emcee

import sys
sys.path.append('../python/')
from EdwRes import *
from prob_dist import *
# GGA3 parameters from Edelweiss tables
ion_center_0keV = 1.3
ion_guard_0keV = 1.5
heat_0keV = 0.4
ion_122keV = 3.1 
heat_122keV = 2.7
aH = 0.0157
pars = {
    'V' : 4.0,
    'eps_eV' : 3.0
    }
from edelweiss_fit import *

# import data from Edelweiss
resNR_data = pd.read_csv("data/edelweiss_NRwidth_GGA3_data.txt", skiprows=1, \
                       names=['E_recoil', 'sig_NR', 'E_recoil_err', 'sig_NR_err'], \
                       delim_whitespace=True)

resER_data = pd.read_csv("data/edelweiss_ERwidth_GGA3_data.txt", skiprows=1, \
                         names=['E_recoil', 'sig_ER', 'sig_ER_err'], \
                         delim_whitespace=True)

# the sorting is necessary!
# otherwise the mask defined below will select the wrong data
resER_data = resER_data.sort_values(by='E_recoil')

#print (res_data.head(4))

# set the data up for the fits
# Edelweiss discards ER points near peaks
# and first two NR points since they're affected by the threshold
mask = [True, True, False, False, True, True, True, True, True]
ER_data = {'Erecoil': resER_data["E_recoil"][mask], 'sigma': resER_data["sig_ER"][mask], 'sigma_err': resER_data["sig_ER_err"][mask]}
NR_data = {'Erecoil': resNR_data["E_recoil"][2::], 'sigma': resNR_data["sig_NR"][2::], 'sigma_err': resNR_data["sig_NR_err"][2::]}<jupyter_output><empty_output><jupyter_text>Fit with an MCMC method<jupyter_code>def log_likelihood(theta, ER_data, NR_data, pars):
    C, m = theta
    aH, scale, A, B = pars
    V = np.abs(scale)*4.0
    
    # extract the data
    x_ER, y_ER, yerr_ER = ER_data['Erecoil'], ER_data['sigma'], ER_data['sigma_err']
    x_NR, y_NR, yerr_NR = NR_data['Erecoil'], NR_data['sigma'], NR_data['sigma_err']
    
    # expected parameter values and widths
    # uncertainty on aH is the uncertainty on the parameter aH when fitting only the ER band
    # scale width estimated by assuming a 10 mV error on V and 0.5 eV error on epsilon
    # information for A and B from Astroparticle Physics 14 (2001) 329±337
    # Anthony's function use the FWHM version of aH
    exp_aH = 0.016*2*np.sqrt(2*np.log(2))
    exp_aH_sig = exp_aH*0.046
    exp_scale = 1
    exp_scale_sig = 0.17
    exp_A = 0.16
    exp_A_sig = 0.07
    exp_B = 0.18
    exp_B_sig = 0.1
    
    ## get the ER prediction for the input parameters
    model_ER = [np.sqrt(series_NRQ_var(x,V=V,aH=aH,A=1.0,B=0.0,alpha=0.0000001)) for x in x_ER]
    sigma2_ER = yerr_ER**2
    
    ## get the NR prediction for the input parameters
    # series_NRQ_var_corr1(Er=10.0,F=0.0,V=4.0,aH=0.0381,alpha=(1/18.0),A=0.16,B=0.18,label='GGA3',corr1file='data/sigdiff_test.h5')
    # series_NRQ_var_corr1 returns the *variance*
    model_NR_0 = [series_NRQ_var_corr1(x, 0, V, aH, 1/18.0, A, B, 'GGA3') for x in x_NR] 
    model_NR = np.sqrt(np.power(C + m*x_NR, 2) + model_NR_0)
    sigma2_NR = yerr_NR**2
    
    return -0.5*(np.sum((y_NR-model_NR)**2/sigma2_NR + np.log(2*np.pi*sigma2_NR)) \
                 + np.sum((y_ER-model_ER)**2/sigma2_ER + np.log(2*np.pi*sigma2_ER)) \
                 + (aH - exp_aH)**2/exp_aH_sig**2 \
                 + (scale - exp_scale)**2/exp_scale_sig**2 \
                 + (A - exp_A)**2/exp_A_sig**2 \
                 + (B - exp_B)**2/exp_B_sig**2)<jupyter_output><empty_output><jupyter_text>Find the minimum likelihoodWhen all six parameters are allowed to vary, the solution finds a small value for A and a larger value for B:aH, C, m, scale, A, B = x: array([1.42450916e-02, 4.23384048e-02, 1.45366242e-05, 1.17278141e+00,       5.26524413e-02, 2.84524357e-01])<jupyter_code># start at the minimum determined by the six-parameter fit
# aH, C, m, scale, A, B = np.array([1.59555935e-02, 4.01182258e-02, 4.10754863e-05, 1.03491851e+00, 1.16818488e-01, 1.51974313e-01])
aH, C, m, scale, A, B = np.array([3.75915355e-02, 3.91363112e-02, 5.25852969e-05, 1.02252631e+00, 9.51212967e-02, 1.45393161e-01])
pars = (aH, scale, A, B)

from scipy import optimize
#np.random.seed(42)

nll = lambda *args: -log_likelihood(*args)
initial = np.array([C, m]) + np.multiply([0.01*C, 0.01*m], np.random.randn(2))
#soln = optimize.minimize(nll, initial, method='Powell', args=(ER_data, NR_data, pars)) #, options={'adaptive':True})
soln = optimize.brute(nll, ([0.03, 0.048], [-0.001, 0.002]), args=(ER_data, NR_data, pars)) #, options={'adaptive':True})

#C_fit, m_fit = soln.x
C_fit, m_fit = soln
print (soln)
# [3.91580254e-02 5.23003616e-05]
Test = True
filename = 'data/edelweiss_C_systematicErrors_sampler_nll_Candm_gausPrior_corr.h5'

# random walks starting near the minimum
ndim, nwalkers = 2, 100
nburn = 500
nsteps = 5000
    
if not Test:
    pos0 = [[C_fit, m_fit] + 1e-8*np.random.randn(ndim) for i in range(nwalkers)]
    sampler = emcee.EnsembleSampler(nwalkers, ndim, log_likelihood, args=(ER_data, NR_data, pars))
    #pos, prob, state = sampler.run_mcmc(pos0, nburn, storechain=False)
    pos, prob, state = sampler.run_mcmc(pos0, nburn)
    sampler.reset()
    sampler.run_mcmc(pos, nsteps)
    samples = sampler.chain
    
else:
    f = h5py.File(filename,'r')

    # need to store data in an array:
    # The sampler will now have a chains attribute 
    # which is an array with shape (Nwalker,N,Ndim) 
    # where N is the number of interations (500 in our inital run)
    # and Ndim is the number of fit parameters
    path='{}/{}/'.format('mcmc','sampler')

    C = np.asarray(f[path+'C'])
    m = np.asarray(f[path+'m'])
    samples = np.asarray(f[path+'samples'])

    f.close()
import os.path

if not Test:
    #save the MCMC results
    path='{}/{}/'.format('mcmc','sampler')

    #remove vars
    f = h5py.File(filename,'a')
    exC = path+'C' in f
    exm = path+'m' in f
    
    exsamp = path+'samples' in f

    if exC:
      del f[path+'C']
    if exm:
      del f[path+'m']
    if exsamp:
      del f[path+'samples']

    dset = f.create_dataset(path+'C',np.shape(C),dtype=np.dtype('float64').type)
    dset[...] = C
    dset = f.create_dataset(path+'m',np.shape(m),dtype=np.dtype('float64').type)
    dset[...] = m
    dset = f.create_dataset(path+'samples',np.shape(samples),dtype=np.dtype('float64').type, \
    compression="gzip",compression_opts=9)
    dset[...] = samples

    f.close()
fig, axes = plt.subplots(ndim, figsize=(10, 7), sharex=True)

print (np.shape(samples))
print (np.shape(samples[:, :, 1]))
labels = ["C", "m"]
for i in range(ndim):
    ax = axes[i]
    # when you pass ax.plot a pair, the columns represent seperate data sets
    ax.plot(np.transpose(samples[:, :, i]), "k", alpha=0.03)
    ax.set_xlim(0, np.shape(samples)[1])
    ax.set_ylabel(labels[i])
    ax.yaxis.set_label_coords(-0.1, 0.5)

axes[-1].set_xlabel("step number");
import corner
print(np.shape(samples))
samples = samples[:, 300:, :].reshape((-1, ndim))
print(np.shape(samples))
fig = corner.corner(samples, labels=["$C$", "$m$"],
                      truths=[C_fit, m_fit])
fig, axs = plt.subplots(1,2,figsize=(15,6))
# the ER-band fit
ax1 = axs[0]
ax2 = axs[1]

Erecoil = np.linspace(10, 200, 50)
   
# GGA3 parameters from Edelweiss tables
V = scale*4.0 
    
## get the ER prediction for the input parameters
## note that the ER prediction is unchanged for changing C, m
model_ER = [np.sqrt(series_NRQ_var(x,V=V,aH=aH,A=1.0,B=0.0,alpha=0.0000001)) for x in Erecoil]
model_NR_0 = [series_NRQ_var_corr1(x, 0, V, aH, 1/18.0, A, B, 'GGA3') for x in Erecoil]

## get the NR fit prediction for the input parameters
# series_NRQ_var_corr1(Er=10.0,F=0.0,V=4.0,aH=0.0381,alpha=(1/18.0),A=0.16,B=0.18,label='GGA3',corr1file='data/sigdiff_test.h5')
# series_NRQ_var_corr1 returns the *variance*
model_NR_0 = [series_NRQ_var_corr1(x, 0, V, aH, 1/18.0, A, B, 'GGA3') for x in Erecoil] 
model_NR_fit = np.sqrt(np.power(C_fit + m_fit*Erecoil, 2) + model_NR_0)

for C, m in samples[np.random.randint(len(samples), size=100)]:
    ## get the NR prediction for the input parameters
    # series_NRQ_var_corr1(Er=10.0,F=0.0,V=4.0,aH=0.0381,alpha=(1/18.0),A=0.16,B=0.18,label='GGA3',corr1file='data/sigdiff_test.h5')
    # series_NRQ_var_corr1 returns the *variance*
    model_NR = np.sqrt(np.power(C + m*Erecoil, 2) + model_NR_0)
    
    ax1.plot(Erecoil, model_ER, color="k", alpha=0.05)
    ax2.plot(Erecoil, model_NR, color="k", alpha=0.05)

# plot the best fit
ax1.plot(Erecoil, model_ER, color="r", lw=2, alpha=0.8)
ax2.plot(Erecoil, model_NR_fit, color="r", lw=2, alpha=0.8)

# plot the data
ax1.errorbar(ER_data['Erecoil'], ER_data['sigma'], yerr=ER_data['sigma_err'], fmt="o", label="data")
ax2.errorbar(NR_data['Erecoil'], NR_data['sigma'], yerr=NR_data['sigma_err'], fmt="o", label="data")

#ax1.plot(E_recoil_arr, sig_er_func(E_recoil_arr), label="Edelweiss ER resolution, C=0.0313")
ax1.set(xlabel="Recoil energy (keV)", ylabel="1-sigma width of ER band")
ax2.set(xlabel="Recoil energy (keV)", ylabel="1-sigma width of NR band")
    
plt.xlim(10, 200)
ax2.set_ylim(0.04, 0.1)
sampsize=10000
Erecoil_keV = np.linspace(10, 200, 50)
evec = np.zeros((np.size(Erecoil_keV), sampsize))
V = scale*4.0 
model_NR_0 = [series_NRQ_var_corr1(x, 0, V, aH, 1/18.0, A, B, 'GGA3') for x in Erecoil_keV]
    
for i, (C, m) in enumerate(samples[np.random.randint(len(samples), size=sampsize)]):
    #for i, (C, m) in enumerate(samples):
    # GGA3 parameters from Edelweiss tables
    model_NR = np.sqrt(np.power(C + m*Erecoil_keV, 2) + model_NR_0)
    
    #print(np.std(v))
    #print(np.mean(v))
    evec[:,i] = model_NR
    
   
# use the best-fit function instead of the mean
# to calculate the up and down limits:
# find the height at which 68% of the function is contained
# the upper and lower are the values that cross that height
#upvec = #np.mean(evec,axis=1)+np.std(evec,axis=1)
#dnvec = #np.mean(evec,axis=1)-np.std(evec,axis=1)
upvec = []
dnvec = []
for sample in evec:
    #print (np.size(sample))
    #print(sample)
    interval = hpd(sample, 0.6827)
    #print (interval)
    upvec.append(interval[0])
    dnvec.append(interval[1])
    #plt.hist(sample)
    #break
    
#print (upvec)
V = scale*4.0
model_NR_0 = [series_NRQ_var_corr1(x, 0, V, aH, 1/18.0, A, B, 'GGA3') for x in Erecoil_keV] 
model_NR_fit = np.sqrt(np.power(C_fit + m_fit*Erecoil_keV, 2) + model_NR_0)
    
plt.plot(Erecoil_keV, model_NR_fit, color="r", lw=2, alpha=0.8)
plt.plot(Erecoil_keV, upvec, color="r",linestyle="--", lw=2, alpha=0.8)
plt.plot(Erecoil_keV, dnvec, color="r",linestyle="--", lw=2, alpha=0.8)

plt.xlim(10, 200)
plt.ylim(0.03,0.12)
plt.errorbar(NR_data['Erecoil'], NR_data['sigma'], yerr=NR_data['sigma_err'], fmt="o", label="data", capsize=0)
if not Test:
    #save the results for the MS fit
    path='{}/{}/'.format('mcmc','edwdata_sys_error')

    filename = 'data/statistical_error_fits_corr.h5'

    #remove vars
    f = h5py.File(filename,'a')
    exCms = path+'Cms' in f
    exm = path+'m' in f
    exsamp = path+'samples' in f
    exsampsize = path+'sampsize' in f
    exEr = path+'Er' in f
    exCsig_u = path+'Csig_u' in f
    exCsig_l = path+'Csig_l' in f
    exSigss = path+'Sigss' in f

    if exCms:
      del f[path+'Cms']
    if exm:
      del f[path+'m']
    if exsamp:
      del f[path+'samples']
    if exsampsize:
      del f[path+'sampsize']
    if exEr:
      del f[path+'Er']
    if exCsig_u:
      del f[path+'Csig_u']
    if exCsig_l:
      del f[path+'Csig_l']
    if exSigss:
      del f[path+'Sigss']

    # save the fit parameters
    dset = f.create_dataset(path+'Cms',np.shape(C_fit),dtype=np.dtype('float64').type)
    dset[...] = C_fit
    dset = f.create_dataset(path+'m',np.shape(m_fit),dtype=np.dtype('float64').type)
    dset[...] = m_fit
    
    # save the MCMC samples
    dset = f.create_dataset(path+'samples',np.shape(samples),dtype=np.dtype('float64').type, \
    compression="gzip",compression_opts=9)
    dset[...] = samples
    dset = f.create_dataset(path+'sampsize',np.shape(sampsize),dtype=np.dtype('float64').type)
    dset[...] = sampsize
    
    # save the with of the NR band
    dset = f.create_dataset(path+'Er',np.shape(Erecoil_keV),dtype=np.dtype('float64').type, \
    compression="gzip")
    dset[...] = Erecoil_keV
    dset = f.create_dataset(path+'Csig_u',np.shape(upvec),dtype=np.dtype('float64').type, \
    compression="gzip",compression_opts=9)
    dset[...] = upvec
    dset = f.create_dataset(path+'Csig_l',np.shape(dnvec),dtype=np.dtype('float64').type, \
    compression="gzip",compression_opts=9)
    dset[...] = dnvec
    dset = f.create_dataset(path+'Sigss',np.shape(Erecoil_keV),dtype=np.dtype('float64').type, \
    compression="gzip",compression_opts=9)
    dset[...] = model_NR_fit

    f.close()<jupyter_output><empty_output><jupyter_text>Only the systematic errorThe above calculation includes both the statistical error (from the uncertainty on the data points) and the systematic error (from the uncertainty on parameters used in the fit like the voltage across the detector and the yield model).In the interest of understanding how future analyses might improve on an estimate of the NR Fano, we examine the uncertainty introduced by the systematic error alone.The statistical error is drawn from stat_uncertainty_C.ipynb.<jupyter_code># get the data data/systematic_error_fits.h5 from edelweiss_C_systematicErrors_allParameters.ipynb
filename_stat_error = 'data/systematic_error_fits_corr.h5'
f = h5py.File(filename_stat_error,'r')

path='{}/{}/'.format('mcmc','edwdata_sys_error')

exCms = np.asarray(f[path+'Cms'])
exslope = np.asarray(f[path+'m'])
exEr = np.asarray(f[path+'Er'])
exCsig_u = np.asarray(f[path+'Csig_u'])
exCsig_l = np.asarray(f[path+'Csig_l'])
exSigss = np.asarray(f[path+'Sigss'])
model_NR_fit = exSigss

f.close()

print (np.shape(exEr))
plt.plot(Erecoil_keV, model_NR_fit, color="r", lw=2, alpha=0.8)
plt.plot(Erecoil_keV, upvec, color="r",linestyle="--", lw=2, alpha=0.8)
plt.plot(Erecoil_keV, dnvec, color="r",linestyle="--", lw=2, alpha=0.8)

# the best-fit is different for the two-parameter (C and m) fit
# than for the six-parameter fit

#plt.plot(exEr, ref_sigNR_func(exEr), color='k', lw=2, alpha=0.8)
a = model_NR_fit + exCsig_u - exSigss
b = model_NR_fit - (-exCsig_l + exSigss)
plt.plot(exEr, exCsig_u, color="k",linestyle="--", lw=2, alpha=0.8)
plt.plot(exEr, exCsig_l, color="k",linestyle="--", lw=2, alpha=0.8)
plt.plot(exEr, exSigss, color="k",linestyle="--", lw=2, alpha=0.8)

plt.xlim(10, 200)
plt.ylim(0.03,0.12)
plt.errorbar(NR_data['Erecoil'], NR_data['sigma'], yerr=NR_data['sigma_err'], fmt="o", label="data", capsize=0)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1600721828632536, 'fraction_numerical': 0.04812190883571715, 'mean_word_length': 4.05677593781683, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '535226', 'license': 'MIT', 'n_tokens_mistral': 6564, 'n_tokens_neox': 6170, 'n_words': 1310, 'path': 'analysis_notebooks/edelweiss_C_statisticalErrors_corr.ipynb', 'repo_name': 'villano-lab/nrFano_paper2019', 'size': '14962'}