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

def MCMC_pdf_samples(self, fNew, num_samples=1000, starting_loc=None, stepsize=0.1, burn_in=1000, Y_metadata=None):
        """
        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
        """
        print("Warning, using MCMC for sampling y*, needs to be tuned!")
        if starting_loc is None:
            starting_loc = fNew
        from functools import partial
        logpdf = partial(self.logpdf, f=fNew, Y_metadata=Y_metadata)
        pdf = lambda y_star: np.exp(logpdf(y=y_star[:, None]))
        #Should be the link function of f is a good starting point
        #(i.e. the point before you corrupt it with the likelihood)
        par_chains = starting_loc.shape[0]
        chain_values = np.zeros((par_chains, num_samples))
        chain_values[:, 0][:,None] = starting_loc
        #Use same stepsize for all par_chains
        stepsize = np.ones(par_chains)*stepsize
        accepted = np.zeros((par_chains, num_samples+burn_in))
        accept_ratio = np.zeros(num_samples+burn_in)
        #Whilst burning in, only need to keep the previous lot
        burnin_cache = np.zeros(par_chains)
        burnin_cache[:] = starting_loc.flatten()
        burning_in = True
        for i in range(burn_in+num_samples):
            next_ind = i-burn_in
            if burning_in:
                old_y = burnin_cache
            else:
                old_y = chain_values[:,next_ind-1]

            old_lik = pdf(old_y)
            #Propose new y from Gaussian proposal
            new_y = np.random.normal(loc=old_y, scale=stepsize)
            new_lik = pdf(new_y)
            #Accept using Metropolis (not hastings) acceptance
            #Always accepts if new_lik > old_lik
            accept_probability = np.minimum(1, new_lik/old_lik)
            u = np.random.uniform(0,1,par_chains)
            #print "Accept prob: ", accept_probability
            accepts = u < accept_probability
            if burning_in:
                burnin_cache[accepts] = new_y[accepts]
                burnin_cache[~accepts] = old_y[~accepts]
                if i == burn_in:
                    burning_in = False
                    chain_values[:,0] = burnin_cache
            else:
                #If it was accepted then new_y becomes the latest sample
                chain_values[accepts, next_ind] = new_y[accepts]
                #Otherwise use old y as the sample
                chain_values[~accepts, next_ind] = old_y[~accepts]

            accepted[~accepts, i] = 0
            accepted[accepts, i] = 1
            accept_ratio[i] = np.sum(accepted[:,i])/float(par_chains)

            #Show progress
            if i % int((burn_in+num_samples)*0.1) == 0:
                print("{}% of samples taken ({})".format((i/int((burn_in+num_samples)*0.1)*10), i))
                print("Last run accept ratio: ", accept_ratio[i])

        print("Average accept ratio: ", np.mean(accept_ratio))
        return chain_values