human
I am working on the file "Bayesian_bootstrap.ipynb".
The first task is:
The Bayesian bootstrap (with a Dirichlet prior)
(See:
http://sumsar.net/blog/2015/04/the-non-parametric-bootstrap-as-a-bayesian-model/
and
http://projecteuclid.org/euclid.aos/1176345338
)
Can you write Python code for it?
assistant
#Here is the Python code you requested.

def bayes_bstrp(data, statistic, nbstrp, samplesize):
    """Implements the Bayesian bootstrap method."""
    
    def Dirichlet_sample(m,n):
        """Returns a matrix of values drawn from a Dirichlet distribution with parameters = 1.
        'm' rows of values, with 'n' Dirichlet draws in each one."""
        # Draw from Gamma distribution
        Dirichlet_params = np.ones(m*n) # Set Dirichlet distribution parameters
        # https://en.wikipedia.org/wiki/Dirichlet_distribution#Gamma_distribution
        Dirichlet_weights = np.asarray([random.gammavariate(a,1) for a in Dirichlet_params])
        Dirichlet_weights = Dirichlet_weights.reshape(m,n) # Fold them (row by row) into a matrix
        row_sums = Dirichlet_weights.sum(axis=1)
        Dirichlet_weights = Dirichlet_weights / row_sums[:, np.newaxis] # Reweight each row to be normalised to 1
        return Dirichlet_weights
    
    Dirich_wgts_matrix = Dirichlet_sample(nbstrp, data.shape[0]) #Generate sample of Dirichlet weights
    
    # If statistic can be directly computed using the weights (such as the mean), do this since it will be faster.
    if statistic==np.mean or statistic==np.average:
        results = np.asarray([np.average(data, weights=Dirich_wgts_matrix[i]) for i in xrange(nbstrp)])
        return results
    
    # Otherwise resort to sampling according to the Dirichlet weights and computing the statistic
    else:
        results = np.zeros(nbstrp)
        for i in xrange(nbstrp): #Sample from data according to Dirichlet weights
            weighted_sample = np.random.choice(data, samplesize, replace=True, p = Dirich_wgts_matrix[i])
            results[i] = statistic(weighted_sample) #Compute the statistic for each sample
        return results