生物統計論 第0回

import numpy as np
import collections
import matplotlib.pyplot as plt

N  = 100000

def BernoulliVar(q, N, seed=None):
    rng  = np.random.RandomState(seed) if seed is not None else np.random
    return np.var(rng.rand(N)<q)

qs = np.linspace(0,1,1000)

probs = [BernoulliVar(q, N) for q in qs]

plt.title("Bernoulli distribution")
plt.plot(qs, probs, label="$Var(X)$")
plt.xlabel("q")
plt.ylabel("$Var(X)$")
plt.legend()
plt.show()

def BinomialProb(n, p, N, seed=None, color=None, ax=None):
    if ax==None:
        fig,ax = plt.subplots()
    rng  = np.random.RandomState(seed) if seed is not None else np.random
    counter = collections.Counter(rng.binomial(n, p, N))
    count, freqs = zip(*sorted(counter.most_common()))
    probs = [freq/N for freq in freqs]
    ax.plot(count, probs, color=color, label=f"p={p} and n={n}")
    ax.scatter(count, probs, color=color)
    return ax

ax = BinomialProb(20, 0.5, N, seed=0, color="blue")
ax = BinomialProb(20, 0.7, N, seed=0, color="green", ax=ax)
ax = BinomialProb(40, 0.5, N, seed=0, color="red",   ax=ax)
plt.legend()
plt.show()

from scipy import stats

def NormalProb(mu, sigma):
    X = np.arange(mu - 3.5*sigma,  mu + 3.5*sigma, sigma/50)
    Y = stats.norm.pdf(X, mu, sigma)
    for i in np.arange(-3,4,1):
        idx = np.argmin(abs(X-(mu+i*sigma)))
        plt.plot([X[idx], X[idx]], [0, Y[idx]], 'k-')
        plt.scatter(X[idx], Y[idx], color="black")
    plt.plot(X,Y, color="black")

NormalProb(mu=0, sigma=0.1)

from scipy.stats import multivariate_normal

x, y = np.mgrid[-3:3.5:.01, -3:3:.01]
pos = np.dstack((x, y))
rv = multivariate_normal([0.5, -0.2], [[2.0, 2.3], [0.3, 0.8]])
plt.contourf(x, y, rv.pdf(pos))
plt.show()