1次元混合ガウスモデルのプログラムの出力値がわからない

夜分投稿失礼いたします。
以下のプログラムは1次元の混合ガウスモデルのプログラムです。
何個のガウス分布を重ね合わせるのが最適かをBICというもので計算しています。

混合ガウスモデルのアルゴリズムは理解しているのですがこのプログラムの出力値が何を示しているか自分の知識不足のせいでわかりません。
自分で考えた中では入力データの負担率だと思うのですがあってますか？
このプログラムにおけるy_trueとy_estimatedは何を表していますか？

負担率とは、あるデータxにおいて、k番目のクラスの正規分布からの観測される確率のことです。別の言い方をすると、あるデータxにおける混合正規分布の割合とも考えられそうです。

３つのクラスA,B,Cに分かれると想定した場合１つのデータを入力すると
[クラスAに属する確率,クラスBに属する確率,クラスCに属する確率] = [0.5, 0.3, 0.2]と出るはずです．

加えて負担率はクラスごとに計算されるものだと思うのですがこのプログラムで負担率はクラスごとに計算されますか？

ちなみに25個のデータを入力したら重ね合わせる最適なガウス分布の個数は1で出力データは入力データの個数と同じ個数出てきました。

入力データ：
[1. 1. 1. 2. 1. 1. 1. 2. 1. 2. 3. 1. 3. 1. 2. 4. 1. 2. 4. 1. 1. 2. 5. 1. 1. ]

出力データ：
<y_true>
[0.41558912 0.57036429 0.52914185 0.48112792 0.60293698 0.46554157
0.52676959 0.54020389 0.37723733 0.44722378 0.56863281 0.58951637
0.51391066 0.51027414 0.58159079 0.45768141 0.45611032 0.43729612
0.55446468 0.48241021 0.45572774 0.55418018 0.51759982 0.44739818
0.45398047]

<y_estimated>
[0.49164472 0.49164472 0.49164472 0.50693811 0.49164472 0.49164472
0.49164472 0.50693811 0.49164472 0.50693811 0.52245777 0.49164472
0.52245777 0.49164472 0.50693811 0.5224578 0.49164472 0.50693811
0.5224578 0.49164472 0.49164472 0.50693811 0.5224578 0.49164472
0.49164472]

<図1>

<図2>

python
```
コード
```
import numpy as np
from scipy.optimize import leastsq
import matplotlib.pyplot as plt
from tqdm import tqdm
import math
import csv
 
"""Setting up test data"""
def gaussian_func(x, A, mu, sigma): 
    return A * np.exp( - (x - mu)**2 / (2 * sigma**2))
 
num_peak_true = 3
A = np.array([0.8, 1, 0.5]) #np.random.rand(num_peak) 
mu = np.array([0.2, 0.4, 0.8])
sigma = np.array([0.05, 0.05, 0.03])
 


list = []
with open('data/src/sample.csv') as f:#ファイル名は架空のものを記載
    reader = csv.reader(f)
    for row in reader:
        list.append(row[3])#sample.csvの4列目の値を取得
float = [float(s) for s in list]#文字列をfloat型に変換
x = np.array(float)#numpy型に変換　xは入力したデータ
 
y_true = np.random.normal(0.5, 0.05, len(x))
for i in range(num_peak_true):
    y_true += gaussian_func(x, A[i], mu[i], sigma[i])
 
 
"""Solving"""
max_num_peak = 6
 
max_iter = 100
bic = 1e5
bic_arr = np.zeros(max_num_peak)
for num_peak in range(1, max_num_peak+1):
    def residual(p, y, x):
      y_fit = p[0]
      for i in range(num_peak):
          y_fit += gaussian_func(x, p[3*i+1], p[3*i+2], p[3*(i+1)])
      error = y - y_fit
      return error
     
    sse = 1e5 # set initial loss large  
    num_param = num_peak*3 + 1
    for i in tqdm(range(max_iter)):
        p = np.random.rand(num_param) # initial guesses for leastsq
        plsq, ier = leastsq(residual, p, args = (y_true, x))
        sse_i = np.sum(residual(plsq, y_true, x)**2)
        if sse_i < sse:
            sse = sse_i
            plsq_opt = plsq
     
    #BIC = n* ln(sse/n) + k*ln(n)

    bic_i = num_sample*math.log(sse/num_sample) + num_param*math.log(num_sample)
    print("num of peak: ", num_peak, "BIC:", bic_i)
    bic_arr[num_peak-1] = bic_i
    if bic_i < bic:
        bic = bic_i
        num_peak_estimated = num_peak
        plsq_global_opt = plsq_opt
  
"""plot results"""
y_estimated = plsq_global_opt[0]
for i in range(num_peak_estimated):    
    y_estimated += gaussian_func(x, plsq_global_opt[3*i+1], plsq_global_opt[3*i+2],plsq_global_opt[3*(i+1)])

#出力データ
print(y_true) 
print(y_estimated)

# gaussian fit
fig, ax = plt.subplots(figsize=(6, 4))
plt.scatter(x, y_true, label='true data')
plt.plot(x, y_estimated, label="estimated curve", color="r")
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_ylim(0, 1.75)
plt.legend()
plt.tight_layout()
plt.savefig("mix_gauss_fit.png")
plt.show()#図１
plt.close()
 
# BIC
x_bic = np.arange(1,max_num_peak+1, 1)
fig, ax = plt.subplots(figsize=(6, 4))
plt.plot(x_bic, bic_arr, color="r")
ax.set_xlabel('Number of Gaussian fit')
ax.set_ylabel('BIC')
plt.xticks(x_bic)
plt.tight_layout()
plt.savefig("bic.png")
plt.show()#図2
plt.close()
```[リンク内容](https://omedstu.jimdofree.com/2018/12/01/scipyによる1次元混合ガウス回帰/)