#  Kernel density estimation
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from sklearn.neighbors import KernelDensity
# Code reference: http://scikit-learn.org/stable/auto_examples/neighbors/
# plot_kde_1d.html
N = 200
np.random.seed(1)
# Create 2 normal distributed data set
norm_data_1 = np.random.normal(0, 1, int(0.3 * N))
norm_data_2 = np.random.normal(5, 1, int(0.7 * N))
norm_data = np.concatenate((norm_data_1, norm_data_2))

X_plot = np.linspace(-5, 10, 1000) # Create x axis range
# Create linear combination of 2 normal distributed random variable
norm_linear = (0.3 * norm(0, 1).pdf(X_plot) + 0.7 * norm(5, 1).pdf(X_plot))

# figure
fig, ax = plt.subplots()
# Plot the real distribution
ax.fill(X_plot, norm_linear, fc='black', alpha=0.2,
        label='Linearcombination')
# Use 3 different kernels to estimate
for kernel in ['gaussian', 'tophat', 'epanechnikov']:
    # Initial an object to use kernel function to fit data,
    # bandwidth will affect the result
    kde = KernelDensity(kernel=kernel, bandwidth=0.5).fit(norm_data.reshape(-1, 1))
    # Evaluate the density model on the data
    log_dens = kde.score_samples(X_plot.reshape(-1, 1))
    ax.plot(X_plot, np.exp(log_dens), '-',
            label="kernel ='{0}'".format(kernel))

# Add text on the plot, position argument can be arbitrary
ax.text(6, 0.38, "N={0} points".format(N))
ax.legend(loc='upper left')
# Plot the random points, squeeze them into narrow space
ax.plot(norm_data, -0.005 - 0.01 *
        np.random.random(norm_data.shape[0]), '+k')

# Set x-axis y-axis limit to adjust the figure
ax.set_xlim(-4, 9)
ax.set_ylim(-0.03, 0.4)
fig.savefig('kernel_estimation.png', dpi=300)
plt.show()

二维散点图：

# Using the Box-Mueller Method to generate 2-dim normally distributed variables
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(100) # Set seed from comparability
# For mu = (0,0), covariance matrix Sigma = identity matrix
n = 500 # Number of random numbers
msize = 0.1 # determines the size of the plotted points

# a good size might be msize=5 for n=500 pts and msize=0.1 for n>50K
a = np.random.exponential(scale=1, size=n)
phi = np.random.uniform(low=0, high=2 * np.pi, size=n)
# change to cartesian coordinates
x = a * np.cos(phi)
y = a * np.sin(phi)
plt.figure(figsize=(4, 4))
plt.plot(x, y, 'ro', markersize=msize)

# for covariance matrix Sigma = A: Y = X/sqrt(Sigma) ~ N(0,I) => Y*sqrt(Sigma)
# Calculate sqrt(A) with Jordan decomposition
A = [[3, 1], [1, 1]]
A_eig = np.linalg.eig(A)

E_val = A_eig[0]
Gamma = A_eig[1]
Lambda = np.diag(E_val)
np.sqrt(Lambda)
Lambda12 = np.sqrt(Lambda)

A12 = np.dot(np.dot(Gamma, Lambda12), np.transpose(Gamma))

# Solve with matrix multiplication
c = [x, y]
tfxy = np.dot(A12, c)

# print(N)
plt.figure(2, figsize=(6, 4))
plt.plot(tfxy[0], tfxy[1], 'ro', markersize=msize)