5

One can create a multivariate kernel density estimate (KDE) with scikitlearn (https://scikit-learn.org/stable/modules/generated/sklearn.neighbors.KernelDensity.html#sklearn.neighbors.KernelDensity) and scipy (https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.gaussian_kde.html)

Both allow for random sampling from the estimated distribution. Is there a way to do conditional sampling in either of the two libraries (or any other library)? In the 2-variable (x,y) case this would mean sample from P(x|y) (or P(y|x)), thus from a cross-section of the probability function (and that cross-section has to be rescaled to unit area under its curve).

x = np.random.random(100)
y =np.random.random(100)
kde = stats.gaussian_kde([x,y])
# sampling from the whole pdf:
kde.resample()

I am looking for something like

# sampling y, conditional on x
kde.sample_conditional(x=1.5) #does not exist
4
  • Have you solved this?
    – vps
    Sep 10, 2020 at 10:53
  • also interested in this.
    – Vlad
    Apr 16, 2021 at 15:36
  • any development? Sep 15, 2021 at 9:12
  • Check my answer Oct 29, 2021 at 9:41

2 Answers 2

1

This function samples conditionally one column based on all other columns.

def conditional_sample(kde, training_data, columns, values, samples = 100):
    if len(values) - len(columns) != 0:
        raise ValueError("length of columns and values should be equal")
    if training_data.shape[1] - len(columns) != 1:
        raise ValueError(f'Expected {training_data.shape[1] - 1} columns/values but {len(columns)} have be given')  
    cols_values_dict = dict(zip(columns, values))
    
    #create array to sample from
    steps = 10000  
    X_test = np.zeros((steps,training_data.shape[1]))
    for i in range(training_data.shape[1]):
        col_data = training_data.iloc[:, i]
        X_test[:, i] = np.linspace(col_data.min(),col_data.max(),steps)
    for col in columns: 
        idx_col_training_data = list(training_data.columns).index(col)
        idx_in_values_list = list(cols_values_dict.keys()).index(col)
        X_test[:, idx_col_training_data] = values[idx_in_values_list] 
        
    #compute probability of each sample
    prob_dist = np.exp(kde.score_samples(X_test))/np.exp(kde.score_samples(X_test)).sum()
    
    #sample using np.random.choice
    return X_test[np.random.choice(range(len(prob_dist)), p=prob_dist, size = (samples))]
2
  • 1
    nice workaround for part of the problem, but not for everything. If I understand it correctly, it works with probing the KDE with a discrete grid, which will introduce errors. Also, a KDE can actually extrapolate beyond the smallest and largest training point along an axis, which this method wont be able to do.
    – Sip
    Nov 2, 2021 at 9:15
  • Indeed for more efficient conditional sampling, we need to construct the pdf i think. Sadly no implementation of pdf in sklearn.kde and no sampling implementation of sampling in statsmodels.kde Nov 2, 2021 at 12:10
0

This function works for use cases with only 1-dimensional conditions. It's built in a way that the condition doesn't need to be fulfilled exactly, but with some tolerance.

from scipy import stats
import numpy as np

x = np.random.random(1000)
y = np.random.random(1000)
kde = stats.gaussian_kde([x,y])

def conditional_sample(kde, condition, precision = 0.5,  sample_size = 1000):

    block_size = int(sample_size / (2 * precision / 100))
    sample = kde.resample(block_size)

    lower_limit = np.percentile(sample[-1], stats.percentileofscore(sample[-1], condition) - precision)
    upper_limit = np.percentile(sample[-1], stats.percentileofscore(sample[-1], condition) + precision)

    conditional_sample = sample[:, (upper_limit > sample[-1]) * (sample[-1] > lower_limit)]

    return conditional_sample[0:-1]

conditional_sample(kde=kde, condition=0.6)

Not a universal solution but works well for my use case.

1
  • 1
    Your answer could be improved with additional supporting information. Please edit to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.
    – Community Bot
    Jan 4, 2022 at 0:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.