MCMCすればいいじゃん！というツッコミはなしでオネガイシャス
ご存知の通り、scipyには多くの確率密度関数があらかじめ定義されており、それらを使うことで、サンプリングや、確率密度関数のプロットなどを容易に行うことができます。
例えば、
stats.norm.rvs(loc=50, scale=20, size=1000)
とすれば、平均50, 標準偏差20の正規分布から1000個のサンプルを得ることができますし、
x = np.linspace(0, 100, 100)
px = stats.norm.pdf(x, loc=50, scale=20)
plt.plot(x, px)
とすれば、$0<x<100$での正規分布を図示することもできます。
では、scipyにあらかじめ定義されていない自前の確率密度関数からサンプリングを行うにはどうしたらよいでしょうか。離散分布の場合は@yk-tanigawaさんの記事で紹介されているので、ここでは連続分布の場合を扱います。
例として、scipy.statsに正規分布が定義されていないとして（本当は定義されてます）、自前でこれを定義して、サンプリングする例を下に書いておきます。rv_continousから継承して、_pdf 関数の中で定義したい関数を書くだけですね。
from scipy import stats
import math
class gaussian(stats.rv_continuous):
def _pdf(self, x, mu, sigma):
normalize_factor = 1.0/(2.0*math.pi*sigma**2)**(1/2)
px = normalize_factor * math.exp(-(x-mu)**2/(2*sigma**2))
return px
gaussian = gaussian(name="gaussian", a=0.0)
sample_from_gaussian = gaussian.rvs(size=1, mu=10.0, sigma=1.0)
なお、自前で定義する確率密度関数は規格化されている必要があるので注意してください。規格化できていない場合は、棄却サンプリングを使うとよいっぽいです。
(しかし、ドキュメントをみても、a=0が必要な理由がよく分からないな…。)