NumPy网格参数化函数评估：解决广播错误与高效计算

import numpy as np

def log_likelihood_function(A, nu_0, alpha, nu, x_i, sigma):
    # 注意：nu和x_i在这里被重新转换为np.array，但通常它们已经是数组
    nu = np.array(nu)
    x_i = np.array(x_i)

    # 核心计算部分，包含了A, nu_0, alpha与nu, x_i的混合运算
    return -len(nu)/2 * np.log(2 * np.pi * sigma**2) - \
           1 / (2 * sigma**2) * np.sum((x_i - A * (nu / nu_0)**alpha * (1 + nu / nu_0)**(-4 * alpha))**2)

def model(A, nu_0, alpha, nu):
    return A * (nu / nu_0)**alpha * (1 + nu / nu_0)**(-4 * alpha)

# 模拟数据生成
nu_data = np.linspace(0.05, 1.0, 500)
x_i_data = model(4.5, 1, 2/3, nu_data) + np.random.normal(0, 0.05, len(nu_data))

# 定义参数范围并创建3D网格
A_range = np.arange(0.0, 10.0, 0.1)  # 100个点
nu_0_range = np.arange(0.0, 5.0, 0.1) # 50个点
alpha_range = np.arange(0.0, 5.0, 0.1) # 50个点
sigma = 0.05

AA, NU_0_MESH, Alpha_MESH = np.meshgrid(A_range, nu_0_range, alpha_range, indexing="ij")

# 尝试直接计算
# L = log_likelihood_function(AA, NU_0_MESH, Alpha_MESH, nu_data, x_i_data, sigma)

ValueError: operands could not be broadcast together with shapes (500,) (100,50,50)

# log_likelihood_function 保持不变，它期望A, nu_0, alpha是标量
def log_likelihood_function(A, nu_0, alpha, nu, x_i, sigma):
    # 确保nu和x_i是数组，这里保留原始代码的转换，但实际中通常已经是
    nu = np.asarray(nu) # 使用asarray更高效，避免不必要的复制
    x_i = np.asarray(x_i)

    # 核心计算：A, nu_0, alpha现在将是标量，与nu, x_i进行广播
    # nu / nu_0 将是 (500,) / scalar，可以广播
    model_output = A * (nu / nu_0)**alpha * (1 + nu / nu_0)**(-4 * alpha)

    return -len(nu)/2 * np.log(2 * np.pi * sigma**2) - \
           1 / (2 * sigma**2) * np.sum((x_i - model_output)**2)

# 定义模型函数
def model(A, nu_0, alpha, nu):
    return A * (nu / nu_0)**alpha * (1 + nu / nu_0)**(-4 * alpha)

# 定义对数似然函数
def log_likelihood_function(A, nu_0, alpha, nu_data_fixed, x_i_data_fixed, sigma_fixed):
    # 注意：这里的nu_data_fixed, x_i_data_fixed, sigma_fixed将是固定值/数组
    # A, nu_0, alpha 将是来自网格的单个标量值

    # 计算模型预测
    model_prediction = model(A, nu_0, alpha, nu_data_fixed)

    # 计算残差平方和
    sum_of_squares = np.sum((x_i_data_fixed - model_prediction)**2)

    # 计算对数似然
    term1 = -len(nu_data_fixed) / 2 * np.log(2 * np.pi * sigma_fixed**2)
    term2 = -1 / (2 * sigma_fixed**2) * sum_of_squares

    return term1 + term2

# 模拟数据生成（与之前相同）
nu_data = np.linspace(0.05, 1.0, 500)
x_i_data = model(4.5, 1, 2/3, nu_data) + np.random.normal(0, 0.05, len(nu_data))

# 定义参数范围并创建3D网格（与之前相同）
A_range = np.arange(0.0, 10.0, 0.1)
nu_0_range = np.arange(0.0, 5.0, 0.1)
alpha_range = np.arange(0.0, 5.0, 0.1)
sigma = 0.05

AA, NU_0_MESH, Alpha_MESH = np.meshgrid(A_range, nu_0_range, alpha_range, indexing="ij")

# 使用 np.vectorize 包装函数
# 'excluded' 参数指定哪些参数不应该被向量化，即它们将作为完整的对象传递给每次函数调用
# 在本例中，nu_data, x_i_data, sigma 应该是固定不变的，而 A, nu_0, alpha 才是网格上的变量
vectorized_log_likelihood = np.vectorize(log_likelihood_function, 
                                         excluded=['nu_data_fixed', 'x_i_data_fixed', 'sigma_fixed'],
                                         otypes=[float]) # 指定输出类型为浮点数

# 调用向量化函数
# 将网格参数作为前三个参数传入，固定数据作为被排除的参数传入
L = vectorized_log_likelihood(AA, NU_0_MESH, Alpha_MESH, 
                              nu_data_fixed=nu_data, 
                              x_i_data_fixed=x_i_data, 
                              sigma_fixed=sigma)

print("计算完成，似然函数结果形状:", L.shape)
# 预期输出: 计算完成，似然函数结果形状: (100, 50, 50)

# 这是一个概念性的尝试，实际操作复杂且易错
# 假设我们想让 nu_data 广播到 (100, 50, 50, 500) 这样的形状
# 然后与 (100, 50, 50) 的参数进行运算。
# 这需要对 nu_data 进行 reshape 或 expand_dims，使其在与参数运算时能够兼容。
# 例如，如果 nu_data 是 (500,)，而 A, nu_0, alpha 是 (100, 50, 50)，
# 那么 model_output 应该是一个 (100, 50, 50, 500) 的数组。
# 这意味着 A, nu_0, alpha 需要在最后一个维度上广播到 500，
# 而 nu_data 需要在前面三个维度上广播到 (100, 50, 50)。
# 这种复杂性正是 np.vectorize 旨在简化的地方。