human
What is a the Python code you requested given the following:
I am working on the file "Lecture/11. 추정 및 검정/7) 베이지안 모수 추정의 예.ipynb".
The first task is:
카테고리 분포의 모수 추정
다음으로 클래스 갯수가 $K$인 카테고리 분포의 모수 $\theta$ 벡터를 베이지안 추정법으로 추정해 본다.
카테고리 분포의 모수의 각 원소는 모두 0부터 1사이의 값을 가지므로 사전 분포는 하이퍼 모수 $\alpha_i=\dfrac{1}{K}$인 디리클리 분포로 한다.
$$ P(\theta) \propto \prod_{k=1}^K \theta_i^{\alpha_i - 1} \;\;\.
(\alpha_i = 1/K , \.
\text{ for all } i) $$
데이터는 모두 독립적인 카테고리 분포의 곱이므로 우도는 다음과 같이 다항 분포가 된다.
$$ p(x_{1},\ldots,x_{N} \mid \theta) = \prod_{i=1}^N  \prod_{k=1}^K \theta_k^{x_{i,k}} $$
베이지안 규칙을 사용하여 사후 분포를 구하면 다음과 같이 갱신된 하이퍼 모수 $\alpha'_i$를 가지는 디리클리 분포가 된다.
$$ 
\begin{eqnarray}
p(\theta \mid x_{1},\ldots,x_{N})
&\propto & p(x_{1},\ldots,x_{N} \mid \theta)  P(\theta) \
&=& \prod_{i=1}^N  \prod_{k=1}^K \theta_k^{x_{i,k}} \cdot \prod_{k=1}^K \theta_i^{\alpha_i - 1}  \
&=& \prod_{k=1}^K  \theta^{\sum_{i=1}^N x_i + \alpha_i − 1}   \
&=& \prod_{k=1}^K  \theta^{N_i + \alpha_i −1}   \
&=& \prod_{k=1}^K  \theta^{\alpha'_i −1}   \
\end{eqnarray}
$$
이 경우에도 conjugated prior 임을 알 수 있다.
갱신된 하이퍼 모수의 값은 다음과 같다.
$$ \alpha'_i = N_i + \alpha_i $$
Can you write Python code for it?
assistant
#Here is the Python code you requested.

def plot_dirichlet(alpha):

    def project(x):
        n1 = np.array([1, 0, 0])
        n2 = np.array([0, 1, 0])
        n3 = np.array([0, 0, 1])
        n12 = (n1 + n2)/2
        m1 = np.array([1, -1, 0])
        m2 = n3 - n12
        m1 = m1/np.linalg.norm(m1)
        m2 = m2/np.linalg.norm(m2)
        return np.dstack([(x-n12).dot(m1), (x-n12).dot(m2)])[0]

    def project_reverse(x):
        n1 = np.array([1, 0, 0])
        n2 = np.array([0, 1, 0])
        n3 = np.array([0, 0, 1])
        n12 = (n1 + n2)/2
        m1 = np.array([1, -1, 0])
        m2 = n3 - n12
        m1 = m1/np.linalg.norm(m1)
        m2 = m2/np.linalg.norm(m2)
        return x[:,0][:, np.newaxis] * m1 + x[:,1][:, np.newaxis] * m2 + n12

    eps = np.finfo(float).eps * 10
    X = project([[1-eps,0,0], [0,1-eps,0], [0,0,1-eps]])
    
    import matplotlib.tri as mtri

    triang = mtri.Triangulation(X[:,0], X[:,1], [[0, 1, 2]])
    refiner = mtri.UniformTriRefiner(triang)
    triang2 = refiner.refine_triangulation(subdiv=6)
    XYZ = project_reverse(np.dstack([triang2.x, triang2.y, 1-triang2.x-triang2.y])[0])

    pdf = sp.stats.dirichlet(alpha).pdf(XYZ.T)
    plt.tricontourf(triang2, pdf)
    plt.axis("equal")
    plt.show()

theta0 = np.array([0.2, 0.6, 0.2])

np.random.seed(0)
x1 = np.random.choice(3, 20, p=theta0)
N1 = np.bincount(x1, minlength=3)
x2 = np.random.choice(3, 100, p=theta0)
N2 = np.bincount(x2, minlength=3)
x3 = np.random.choice(3, 1000, p=theta0)
N3 = np.bincount(x3, minlength=3)

a0 = np.ones(3) / 3
plot_dirichlet(a0)

a1 = a0 + N1
plot_dirichlet(a1)
print((a1 - 1)/(a1.sum() - 3))

a2 = a1 + N2
plot_dirichlet(a2)
print((a2 - 1)/(a2.sum() - 3))

a3 = a2 + N3
plot_dirichlet(a3)
print((a3 - 1)/(a3.sum() - 3))