Compute parameters of a PDF (probability density function) for which no closed form expression is available
Your problem could be solved either by direct numeric integration or by MCMC.

Numeric integration can be performed most easily by scipy:

import numpy as np
import scipy.stats
import scipy.integrate

def weird_density(x):
    """ The function I want to sample """
    return scipy.stats.lognorm.pdf(x, 1.0)

def quad_quantile(fun, q, precision=1e-10, minus_inf=-10e10, lower=-10e10, upper=10e10):
    """ Use bisection to evaluate a quantile """
    while upper - lower &gt; precision:
        med = (upper + lower) * 0.5
        cdf_med = scipy.integrate.quad(fun, minus_inf, med)[0]
        if cdf_med &lt; q:
            lower = med
        else:
            upper = med
    return (upper + lower) * 0.5

print('true mean  :', scipy.stats.lognorm(1).mean())
print('integrated mean:', scipy.integrate.quad(lambda x: weird_density(x) * x, 0, 100)[0])
print('true median  :', scipy.stats.lognorm(1).ppf(0.5))
print('integrated median:', quad_quantile(weird_density, 0.5, 0.001, 0, 0, 100))


The output is 

true mean  : 1.6487212707
integrated mean: 1.6484641126903046
true median  : 1.0
integrated median: 0.9998321533203125


However, if the distribution is multidimensional, direct integration may be intractable. In this case,  MCMC methods (Markov chain Monte Carlo) can help.

What they do is just make a correct (in some sense) sample from the distribution for which you know the PDF. When you have a sample, you can calculate all your parameters from it as classical sample statistics, just like from any observed data.

MCMC algorithms may be implemented manually, like in the example below. Another option is to use the PyMC3 library or its analogs.

One of the best known MCMC algorithms, Metropolis-Hastings, works as follows:

def sample_next(x):
    """ Generate random next point"""
    return scipy.stats.norm(loc=x).rvs(1)[0]

def metropolis_hastings(density, generator, n_samples, starting_point=1, random_state=None):
    """ Generate sample from the given density function """
    np.random.seed(random_state)
    result = [starting_point]
    for i in range(n_samples):
        current_point = result[-1]
        next_candidate = generator(current_point)
        acceptance_ratio = density(next_candidate) / density(current_point)
        if np.random.uniform() &lt;= acceptance_ratio:
            result.append(next_candidate)
        else:
            result.append(current_point)
    return np.array(result)

sample = metropolis_hastings(weird_density, sample_next, 1000, random_state=1)

print('sample mean:', np.mean(sample))
print('sample median:', np.median(sample))


This code prints

sample mean: 1.2872857165
sample median: 0.932826883898


You see that the results are not very close, but with longer sampling (1000 points is too few) they will converge. And the shape of distribution is already well matched:

import matplotlib.pyplot as plt
plt.hist(sample, bins=30, normed=True)
idx = np.linspace(0, sample.max())
plt.plot(idx, scipy.stats.lognorm(1.0).pdf(idx))
plt.legend(['true density', 'sample density'])
plt.title('convergence of MCMC sample')
plt.show()