<jupyter_start><jupyter_text>$Like = exp(-\sum \frac{(y(x) - y_{exp})^2}{2\sigma^2})$ $like = e^{-(y(x)-y_{exp})(YcovM)(y(x)-y_{exp}).T}$$YcovM = Ycov - Y_{experror}I$<jupyter_code>#creating the likelihood function
Y_exp = df3.loc[:, features].values
y_exp_std =df4.loc[:,error_features].values
covM = Ycov - y_exp_std*np.eye(33)
Z = Y-Y_exp
H = np.dot(Z,covM) #matmul?
J = np.dot(H,Z.T) #J = (Y-Y_exp)*covM*(Y-Y_exp).T, Should be inverse covM
like = np.exp(-sum(J))
#like = np.exp(-sum(((Y-Y_exp)**2)/(2*y_exp_std**2)))
plt.hist(like,5)
print(like.shape)<jupyter_output>(47,)<jupyter_text>$Posterior = \frac{Like*S0*L*ms*mv}{\sum(like*(S0+L+ms+mv)}\\\\posterior \propto Like*S0*L*ms*mv$<jupyter_code>N = 10000 #number of iterations

#size of the steps
met_step = 1
met_step2 = 0.05
met_step3 = 0.01
met_accept = 0

#empty arrays to store the data from the mcmc
S0 = np.zeros(N)
L = np.zeros(N)
ms = np.zeros(N)
mv = np.zeros(N)
lik = np.zeros(N)
mu = np.zeros(N)

#this is the likelihood function that I derived in the above line
data = like

#starting points
S0[0] = 30
L[0] = 70
mu[0] = 0.8
lik[0] = 0.9
mv[0] = 1

#Here is the random walk
for i in range(N-1):
    if np.random.rand() > 0.5:
        S0_c = S0[i] + uniform(0,met_step).rvs()
        L_c = L[i] + uniform(0,met_step).rvs()
        mv_c = mv[i] + uniform(0,met_step2).rvs()
    else:
        S0_c = S0[i] - uniform(0,met_step).rvs()
        L_c = L[i] - uniform(0,met_step).rvs()
        mv_c = mv[i] - uniform(0,met_step2).rvs()
       
    mu_c = norm(mu[i],met_step3).rvs()    
    ms_i = norm(0.7, 0.05).pdf(mu[i])
    ms_c = norm(0.7, 0.05).pdf(mu_c)
    lik_i = norm(mu[i],1.4).pdf(data).prod()
    lik_c = norm(mu_c,1.4).pdf(data).prod()
    
    #Here we have set up the parameter range
    constraints = [0.6<=ms_c<=1,
                 0.6<mv_c<1.2,
                 25.7<S0_c<36,
                 32<=L_c<=120]
    if all(constraints):
        #the metroplois ratio
        met_r = lik_c*S0_c*L_c*ms_c*mv_c/lik_i*S0[i]*L[i]*ms_i*mv[i]
        if np.random.rand() < min(1,met_r):
            S0[i+1] = S0_c
            L[i+1] = L_c
            mu[i+1] = mu_c
            mv[i+1] = mv_c
            ms[i] = ms_c
            lik[i] = lik_c
            met_accept = met_accept +1
        else:
            S0[i+1] = S0[i]
            L[i+1] = L[i]
            mu[i+1] = mu[i]
            mv[i+1] = mv[i]
            ms[i] = ms_i
            lik[i] = lik_i
    else:
        S0[i+1] = S0[i]
        L[i+1] = L[i]
        mu[i+1] = mu[i]
        mv[i+1] = mv[i]
        ms[i] = ms_i
        lik[i] = lik_i
        
        
M = met_accept/N
print(M)
#Using the MCMC sampling figure 2 can be created using the plots below
#plt.hist(S0,30, range = (10,40))
#plt.hist(L,30, range=(30,120))
#plt.hist(ms,30,range =(0.6,1))
#plt.hist(mv,100, range = (0.6,1.2))
#plt.hist2d(ms,mv)
#plt.hist2d(ms,L)
#plt.hist2d(ms,S0)
#plt.hist2d(mv,ms)
#plt.hist2d(mv,L)
#plt.hist2d(mv,S0)
#plt.hist2d(L,ms)
#plt.hist2d(L,mv)
#plt.hist2d(L,S0)
#plt.hist2d(S0,mv)
#plt.hist2d(S0,L)
#plt.hist2d(S0,ms)
#using the ms and mv sampling from mcmc, figure 3 can be created
fI = 1/ms - 1/mv
plt.hist(fI, 100, range=(-1,1))
print(fI.mean())
#reshaping the sampling from mcmc to look like the inputs
S0 = S0.reshape(N,1)
L = L.reshape(N,1)
ms = ms.reshape(N,1)
mv = mv.reshape(N,1)
SL = np.concatenate((S0,L),axis=1)
msv = np.concatenate((ms,mv),axis=1)
x_mcmc = np.concatenate((SL,msv),axis=1)
x_MCMC = pd.DataFrame(data = x_mcmc, columns =  ['S0', 'L', 'ms', 'mv'])
x_MCMC.head()
#using the trained gaussian process to predict a new y based on the mcmc sampling. This new y will be the posterior
y_mcmc = gp.predict(x_mcmc, return_std=True)
#Reverse PCA on y_mcmc to look like y
Y_mcmc = np.dot(y_mcmc, pca.components_) + pca.mean_
Y_mcmcdf= pd.DataFrame(data = Y_mcmc, columns =  features)
Y_mcmcdf.head()
#posterior mean
Y_mcmc_mean = Y_mcmcdf[features].mean()
#here figure 1 is created
#plotting the Y_mcmc based on energy E_c.m. for example the column SN112EK42.5 represent energy 42.5 
Yy_mcmc = np.array([Y_mcmc[:,0],Y_mcmc[:,1],Y_mcmc[:,2],Y_mcmc[:,3],Y_mcmc[:,4],Y_mcmc[:,5],Y_mcmc[:,6],Y_mcmc[:,7],Y_mcmc[:,8],Y_mcmc[:,9],Y_mcmc[:,10],Y_mcmc[:,11],Y_mcmc[:,12],Y_mcmc[:,13],Y_mcmc[:,14],Y_mcmc[:,15],Y_mcmc[:,16],Y_mcmc[:,17],Y_mcmc[:,18],Y_mcmc[:,19],Y_mcmc[:,20],Y_mcmc[:,21],Y_mcmc[:,22],Y_mcmc[:,23],Y_mcmc[:,24],Y_mcmc[:,25],Y_mcmc[:,26],Y_mcmc[:,27],Y_mcmc[:,28],Y_mcmc[:,29],Y_mcmc[:,30],Y_mcmc[:,31],Y_mcmc[:,32]])
plt.plot(features,Yy_mcmc, 'g')
plt.plot(features, Y_mcmc_mean, 'k')
plt.plot(features,Yy, 'r')
plt.errorbar(features,Y_ex, yerr = y_ex_std)
plt.xlabel('$E_{c.m.}$ (MeV)')
plt.ylabel('$R_{n/p}$ & $DR_{n/p}$')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20530785562632697, 'fraction_numerical': 0.052441613588110406, 'mean_word_length': 3.521113243761996, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 9, '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': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '250498', 'license': 'CC0-1.0', 'n_tokens_mistral': 2397, 'n_tokens_neox': 2303, 'n_words': 385, 'path': 'doc/Programs/DimRed/Bayes.ipynb', 'repo_name': 'anacost/MachineLearning', 'size': '4710'}