如何在Scipy中计算位置-尺度t分布的联合似然？

现有代码可计算已知σ的正态分布下3个变量的联合似然，现在需要扩展到包含位置（loc）、尺度（scale）、自由度（df）三个未知参数的位置-尺度t分布，但手动指定参数域时遇到困难，想知道有没有可用的辅助函数。

已知σ的正态分布似然计算代码 #

n = 3; x = stats.norm(loc=0,scale=1).rvs(n, random_state = 10); # set seed
parameter_domain = np.linspace(-10,10,1001);
likelihood = stats.norm.pdf((x[:,np.newaxis])*np.ones(parameter_domain.shape),
                            loc=parameter_domain1, scale=1)#.prod(axis=0) # sigma=1 known

尝试的位置-尺度t分布似然计算代码 #

n = 3
df_true = 15, loc_true = 10, scale_true = 2 
x = stats.t(df = df_true, loc = loc_true, scale = scale_true).rvs(n, random_state = 10) # Data from location-scale-t distn.
df_parameter_domain = np.linspace(1, 100, 100)#np.linspace(-10,10,1001);
loc_parameter_domain = np.linspace(-10,10,1001)
scale_parameter_domain = np.linspace(0, 10, 1001)

likelihood = stats.t.pdf((x[:,np.newaxis])*np.ones(parameter_domain.shape),
                            loc=df_parameter_domain, scale=scale_parameter_domain, df = df_parameter_domain)#.prod(axis=0) # no variables known

解决建议 #

1. 合理设置参数域的方法 #

• 自由度（df）：从1开始，上限可根据数据分布调整——如果数据接近正态，设到30-100即可（df越大越接近正态分布），也可以用对数刻度来捕捉小df的细节：

  df_domain = np.logspace(0, 2, 50)  # 生成1到100的对数刻度参数点
  

• 位置参数（loc）：不要盲目设固定范围，可根据数据的均值和标准差来确定，比如取数据均值±3倍标准差的区间：

  loc_domain = np.linspace(x.mean() - 3*x.std(), x.mean() + 3*x.std(), 1001)
  

• 尺度参数（scale）：必须大于0，从极小正数（比如0.01）开始，到数据标准差的2-3倍即可：

  scale_domain = np.linspace(0.01, 3*x.std(), 1001)
  

2. 修正联合似然计算的维度问题 #

原代码存在参数维度不匹配的问题，需要生成所有参数的组合网格，再计算每个组合的联合似然：

import numpy as np
from scipy import stats

n = 3
df_true, loc_true, scale_true = 15, 10, 2 
x = stats.t(df=df_true, loc=loc_true, scale=scale_true).rvs(n, random_state=10)

# 合理设置参数域
df_domain = np.linspace(1, 50, 50)
loc_domain = np.linspace(x.mean()-3*x.std(), x.mean()+3*x.std(), 100)
scale_domain = np.linspace(0.01, 3*x.std(), 100)

# 生成三维参数网格
df_mesh, loc_mesh, scale_mesh = np.meshgrid(df_domain, loc_domain, scale_domain, indexing='ij')

# 计算每个数据点在所有参数组合下的似然，再求乘积得到联合似然
likelihoods = stats.t.pdf(x[:, np.newaxis, np.newaxis, np.newaxis], 
                          df=df_mesh, loc=loc_mesh, scale=scale_mesh)
joint_likelihood = np.prod(likelihoods, axis=0)

3. 利用参数估计辅助设置参数域 #

可以先用scipy.stats.t.fit()得到参数的初始估计值，再围绕估计值设置参数域，避免盲目指定范围：

# 获取参数初始估计
df_est, loc_est, scale_est = stats.t.fit(x)

# 围绕估计值设置合理的参数域
df_domain = np.linspace(max(1, df_est - 10), df_est + 20, 50)
loc_domain = np.linspace(loc_est - 5, loc_est + 5, 100)
scale_domain = np.linspace(max(0.01, scale_est - 2), scale_est + 3, 100)

内容的提问来源于stack exchange，提问作者Ator