2

我有一组 N 个样本（N~10000 到 100000）：

(y_i, P_i)

他们对未知函数进行采样：

y = f(P)

在我的例子中，P 是一组 n_p 参数，其中 n_p 通常在 10 左右。

我的目标是使用样本找到最接近未知函数的 2 阶多项式。其典型结果将是最佳拟合多项式的 (n_p+1)+(n_p+1)*n_p/2 系数。

但是，我想以下列形式获得解决方案：

f(P) = (P-mu)*(C^-1)*(P-mu)^T + K

mu 为 n_p 维向量，C 为 (n_p X n_p) 对称矩阵，K 为常数。这些 mu、C 和 K 是我想要获得的。

Python 生态系统中的某个地方是否有标准实现和/或有效的方法来做到这一点？

提前致谢 ！

编辑： 这是一个典型的设置，我创建假样本（P 和 Y），并且只使用它们，我想恢复原始的 mu、C 和 K：

import numpy as np

n_p = 3
N = 15

# Generate the "mu" we'll try to recover
mu = np.random.rand(n_p)

# Generate the "C" we'll try to recover
tmp = np.random.rand(n_p,n_p)
C = np.dot(tmp, tmp.T)

# Generate the "K" we'll try to recover
K = np.random.rand()

# Generate the samples
P = np.random.rand(n_p,N)
Y = np.array([np.dot(np.dot((P[:,i]-mu),np.linalg.inv(C)),(P[:,i]-mu))+K for i in range(N)])
4

1 回答 1

0

你在找这样的东西吗？我不是 100% 确定它对您的用例是准确的（我不知道当您拟合 2 阶多项式时您可以施加哪些额外的约束），但它应该尝试找到一个 U、A（你称为 C) 和 K 最小二乘误差。

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt

samples = 100
num_params = 20

y = np.random.rand(samples)
p = np.random.rand(samples, num_params)


def my_func(params):
    u = params[0:num_params]
    u = np.expand_dims(u, axis=1)
    a = params[num_params:-1]
    k = params[-1]
    a = a.reshape(num_params, num_params)
    a_inv = np.linalg.inv(a)
    shifted_p = p - np.transpose(u)
    mult_with_a_inv = np.dot(shifted_p, a_inv)
    mat_mult_vec = np.einsum('ij,ji->i', mult_with_a_inv, np.transpose(shifted_p))
    return_val = y - k - mat_mult_vec
    return sum(return_val**2)


guess = np.random.rand(num_params+num_params**2+1)
new_params = optimize.fmin_cg(my_func, guess)
new_u = new_params[0:num_params]
new_a = new_params[num_params:-1]
new_a = new_a.reshape(num_params, num_params)
new_k = new_params[-1]

original_y, = plt.plot(np.arange(samples), y)
new_y, = plt.plot(np.arange(samples), np.einsum('ij,ji->i', np.dot(p-np.transpose(new_u),
                                                                   np.linalg.inv(new_a)),
                                                            np.transpose(p-np.transpose(new_u))
                                 ) + new_k)
plt.legend([original_y, new_y], ['Original Y', 'Newly Fitted Y'])
plt.show()

可能有一种在数学上更合理的方法，但希望这至少会有所帮助。

编辑：刚刚注意到我弄乱了 k 的符号。我还需要增加参数的数量。不过，这似乎工作得非常好。我几乎完全恢复了原始噪音。

还添加示例输出。

样本输出

于 2017-07-26T23:09:03.957 回答