Kurtosis of a normal mixture distribution: error in implementation

Question

I have been using the text Fruhwirth-Shnatter (2006) for analytic formulas on normal mixture distributions. I implemented the formula for the fourth moment but, after converting it to excess Kurtosis, I am unable to get close empirical results with Monte Carlo sampling. What is the error in my code / the formula?

Here a minimal reproducible example of my code.

import numpy as np
from scipy.stats import kurtosis


get_p2 = lambda p1: 1 - p1

def simulate(p1, mu1, mu2, sigma1, sigma2, n):
    
    param_draws = np.random.choice([0, 1], p = [p1, get_p2(p1)], size = n)
    mu_param = np.where(param_draws == 0, mu1, mu2)
    sigma_param = np.where(param_draws == 0, sigma1, sigma2)
        
    return np.array([np.random.normal(loc = mu_param[i], scale = sigma_param[i]) for i in range(n)])

def mixture_mu(p1, mu1, mu2):
    return p1 * mu1 + get_p2(p1) * mu2

def mixture_variance(p1, mu1, mu2, sigma1, sigma2):
    mu = mixture_mu(p1, mu1, mu2)
    part_1 = p1 * (np.power(sigma1, 2) + np.power(mu1, 2)) - np.power(mu, 2)
    part_2 = get_p2(p1) * (np.power(sigma2, 2) + np.power(mu2, 2)) - np.power(mu, 2)
    return part_1 + part_2

def excess_kurtosis(p1, mu1, mu2, sigma1, sigma2):
    mu = mixture_mu(p1, mu1, mu2)
    variance = mixture_variance(p1, mu1, mu2, sigma1, sigma2)
    
    part_1 = np.power(mu1 - mu, 4) + (6 * np.power(mu1 - mu, 2) * np.power(sigma1, 2)) + 3 * np.power(sigma1, 4)
    part_2 = np.power(mu2 - mu, 4) + (6 * np.power(mu2 - mu, 2) * np.power(sigma2, 2)) + 3 * np.power(sigma2, 4)

    return (p1 * part_1 + get_p2(p1) * part_2) / np.power(variance, 2) - 3


p1 = 0.5; mu1 = 0; mu2 = 1; sigma1 = 2; sigma2 = 3

X = simulate(p1, mu1, mu2, sigma1, sigma2, n = 50_000)

analytic_excess_kurtosis = excess_kurtosis(p1, mu1, mu2, sigma1, sigma2)
empirical_excess_kurtosis = kurtosis(X)

print("Analytic:", analytic_excess_kurtosis)
print("Empirical:", empirical_excess_kurtosis)

The deviations from empirical get larger as mu1 and mu2 deviate more.

What is going on? Is it a simple bug somewhere in my code or something else? Thanks in advance.