def MCMC_pdf_samples(self, fNew, num_samples=<NUM_LIT:1000>, starting_loc=None, stepsize=<NUM_LIT:0.1>, burn_in=<NUM_LIT:1000>, Y_metadata=None):
print("<STR_LIT>")<EOL>if starting_loc is None:<EOL><INDENT>starting_loc = fNew<EOL><DEDENT>from functools import partial<EOL>logpdf = partial(self.logpdf, f=fNew, Y_metadata=Y_metadata)<EOL>pdf = lambda y_star: np.exp(logpdf(y=y_star[:, None]))<EOL>par_chains = starting_loc.shape[<NUM_LIT:0>]<EOL>chain_values = np.zeros((par_chains, num_samples))<EOL>chain_values[:, <NUM_LIT:0>][:,None] = starting_loc<EOL>stepsize = np.ones(par_chains)*stepsize<EOL>accepted = np.zeros((par_chains, num_samples+burn_in))<EOL>accept_ratio = np.zeros(num_samples+burn_in)<EOL>burnin_cache = np.zeros(par_chains)<EOL>burnin_cache[:] = starting_loc.flatten()<EOL>burning_in = True<EOL>for i in range(burn_in+num_samples):<EOL><INDENT>next_ind = i-burn_in<EOL>if burning_in:<EOL><INDENT>old_y = burnin_cache<EOL><DEDENT>else:<EOL><INDENT>old_y = chain_values[:,next_ind-<NUM_LIT:1>]<EOL><DEDENT>old_lik = pdf(old_y)<EOL>new_y = np.random.normal(loc=old_y, scale=stepsize)<EOL>new_lik = pdf(new_y)<EOL>accept_probability = np.minimum(<NUM_LIT:1>, new_lik/old_lik)<EOL>u = np.random.uniform(<NUM_LIT:0>,<NUM_LIT:1>,par_chains)<EOL>accepts = u < accept_probability<EOL>if burning_in:<EOL><INDENT>burnin_cache[accepts] = new_y[accepts]<EOL>burnin_cache[~accepts] = old_y[~accepts]<EOL>if i == burn_in:<EOL><INDENT>burning_in = False<EOL>chain_values[:,<NUM_LIT:0>] = burnin_cache<EOL><DEDENT><DEDENT>else:<EOL><INDENT>chain_values[accepts, next_ind] = new_y[accepts]<EOL>chain_values[~accepts, next_ind] = old_y[~accepts]<EOL><DEDENT>accepted[~accepts, i] = <NUM_LIT:0><EOL>accepted[accepts, i] = <NUM_LIT:1><EOL>accept_ratio[i] = np.sum(accepted[:,i])/float(par_chains)<EOL>if i % int((burn_in+num_samples)*<NUM_LIT:0.1>) == <NUM_LIT:0>:<EOL><INDENT>print("<STR_LIT>".format((i/int((burn_in+num_samples)*<NUM_LIT:0.1>)*<NUM_LIT:10>), i))<EOL>print("<STR_LIT>", accept_ratio[i])<EOL><DEDENT><DEDENT>print("<STR_LIT>", np.mean(accept_ratio))<EOL>return chain_values<EOL>
Simple implementation of Metropolis sampling algorithm

Will run a parallel chain for each input dimension (treats each f independently)
Thus assumes f*_1 independant of f*_2 etc.

:param num_samples: Number of samples to take
:param fNew: f at which to sample around
:param starting_loc: Starting locations of the independant chains (usually will be conditional_mean of likelihood), often link_f
:param stepsize: Stepsize for the normal proposal distribution (will need modifying)
:param burnin: number of samples to use for burnin (will need modifying)
:param Y_metadata: Y_metadata for pdf