python中模型的协方差矩阵

2024-09-29 20:17:41 发布

您现在位置:Python中文网/ 问答频道 /正文

在python中求拟合模型的协方差矩阵(相当于python中的vcov()(R function))

lmfit <- lm(formula = Y ~ X, data=Data_df)
lmpred <- predict(lmfit, newdata=Data_df, se.fit=TRUE, interval = "prediction")
std_er <- sqrt(((X0) %*% vcov(lmfit)) %*% t(X0))

尝试用python转换上面的代码。为此,我需要找到拟合模型的协方差矩阵,即vcov。 我不能使用np.cov公司()试图找到模型的协方差矩阵。在

我已经用过了statsmodels.回归.线性模型.OLSResults.cov_参数(),但我得到的值与R中的值不同


Tags: 模型dfdatafunction矩阵covlmlmfit
2条回答

scipy ODR代码可以独立计算参数协方差矩阵,这里是从my的源代码中提取的一个示例尊尊网在线曲线拟合员:

from scipy.optimize import curve_fit
import numpy as np
import scipy.odr
import scipy.stats

x = np.array([5.357, 5.797, 5.936, 6.161, 6.697, 6.731, 6.775, 8.442, 9.861])
y = np.array([0.376, 0.874, 1.049, 1.327, 2.054, 2.077, 2.138, 4.744, 7.104])

def f(x,b0,b1):
    return b0 + (b1 * x)


def f_wrapper_for_odr(beta, x): # parameter order for odr
    return f(x, *beta)

parameters, cov= curve_fit(f, x, y)

model = scipy.odr.odrpack.Model(f_wrapper_for_odr)
data = scipy.odr.odrpack.Data(x,y)
myodr = scipy.odr.odrpack.ODR(data, model, beta0=parameters,  maxit=0)
myodr.set_job(fit_type=2)
parameterStatistics = myodr.run()
df_e = len(x) - len(parameters) # degrees of freedom, error
cov_beta = parameterStatistics.cov_beta # parameter covariance matrix from ODR
sd_beta = parameterStatistics.sd_beta * parameterStatistics.sd_beta
ci = []
t_df = scipy.stats.t.ppf(0.975, df_e)
ci = []
for i in range(len(parameters)):
    ci.append([parameters[i] - t_df * parameterStatistics.sd_beta[i], parameters[i] + t_df * parameterStatistics.sd_beta[i]])

tstat_beta = parameters / parameterStatistics.sd_beta # coeff t-statistics
pstat_beta = (1.0 - scipy.stats.t.cdf(np.abs(tstat_beta), df_e)) * 2.0    # coef. p-values

for i in range(len(parameters)):
    print('parameter:', parameters[i])
    print('   conf interval:', ci[i][0], ci[i][1])
    print('   tstat:', tstat_beta[i])
    print('   pstat:', pstat_beta[i])
    print()

print('Covariance matrix:')    
print(cov_beta)

请提供您正在使用的具体细节。在

假设您对数据使用numpy数组,这里有^{}估计器

相关问题 更多 >

    热门问题