2

polyfit を使用して一連のデータに最適な直線を見つけようとしていますが、パラメーターの不確実性も知る必要があるため、共分散行列も必要です。オンライン ドキュメントでは、次のように書くことを提案しています。

polyfit(x, y, 2, cov=真)

しかし、これはエラーを与えます:

TypeError: polyfit() が予期しないキーワード引数 'cov' を取得しました

確かに、help(polyfit) はキーワード引数 'cov' を表示しません。

オンライン ドキュメントは numpy の以前のリリースを参照していますか? （私は最新の1.6.1を持っています）。独自のポリフィット関数を作成できますが、ポリフィットに共分散オプションがない理由について誰か提案がありますか?

ありがとう

4

1 に答える 1

1

ライブラリからのソリューションの場合、 scikits.statsmodelsを使用するのが便利な選択であることがわかりました。statsmodels では、回帰オブジェクトには、パラメーターと標準エラーを返す呼び出し可能な属性があります。これがどのように機能するかの例を以下に示します。

# Imports, I assume NumPy for forming your data.
import numpy as np
import scikits.statsmodels.api as sm

# Form the data here
(X, Y) = ....

reg_x_data =  np.ones(X.shape);                      # 0th degree term. 
for ii in range(1,deg+1):
    reg_x_data = np.hstack(( reg_x_data, X**(ii) )); # Append the ii^th degree term.


# Store OLS regression results into `result`
result = sm.OLS(Y,reg_x_data).fit()


# Print the estimated coefficients
print result.params

# Print the basic OLS standard error in the coefficients
print result.bse

# Print the estimated basic OLS covariance matrix
print result.cov_params()   # <-- Notice, this one is a function by convention.

# Print the heteroskedasticity-consistent standard error
print result.HC0_se

# Print the heteroskedasticity-consistent covariance matrix
print result.cov_HC0

オブジェクトには追加の堅牢な共分散属性resultもあります。印刷してご覧いただけますdir(result)。また、慣例により、異分散一貫性推定量の共分散行列は、対応する標準誤差を既に呼び出した後にのみ使用できます。たとえば、最初の参照によって 2 番目の参照が計算されて保存されるため、result.HC0_se前に呼び出す必要があります。result.cov_HC0

Pandasは、おそらくこれらの操作に対してより高度なサポートを提供する別のライブラリです。

非ライブラリ関数

これは、追加のライブラリ関数に依存したくない/依存できない場合に役立つ場合があります。

以下は、OLS 回帰係数を返すために作成した関数と、さまざまなものです。残差、回帰分散と標準誤差 (残差の二乗の標準誤差)、大規模サンプル分散の漸近式、OLS 共分散行列、不均一分散一貫性のある "ロバスト" 共分散推定 (OLS 共分散) を返します。ただし、残差に応じて重み付けされます)、および「ホワイト」または「バイアス補正された」異分散一貫性共分散。

import numpy as np

###
# Regression and standard error estimation functions
###
def ols_linreg(X, Y):
    """ ols_linreg(X,Y) 

        Ordinary least squares regression estimator given explanatory variables 
        matrix X and observations matrix Y.The length of the first dimension of 
        X and Y must be the same (equal to the number of samples in the data set).

        Note: these methods should be adapted if you need to use this for large data.
        This is mostly for illustrating what to do for calculating the different
        classicial standard errors. You would never really want to compute the inverse
        matrices for large problems.

        This was developed with NumPy 1.5.1.
    """
    (N, K) = X.shape
    t1 = np.linalg.inv( (np.transpose(X)).dot(X) )
    t2 = (np.transpose(X)).dot(Y)

    beta = t1.dot(t2)
    residuals = Y - X.dot(beta)
    sig_hat = (1.0/(N-K))*np.sum(residuals**2)
    sig_hat_asymptotic_variance = 2*sig_hat**2/N
    conv_st_err = np.sqrt(sig_hat)

    sum1 = 0.0
    for ii in range(N):
        sum1 = sum1 + np.outer(X[ii,:],X[ii,:])

    sum1 = (1.0/N)*sum1
    ols_cov = (sig_hat/N)*np.linalg.inv(sum1)

    PX = X.dot(  np.linalg.inv(np.transpose(X).dot(X)).dot(np.transpose(X))   )
    robust_se_mat1 = np.linalg.inv(np.transpose(X).dot(X))
    robust_se_mat2 = np.transpose(X).dot(np.diag(residuals[:,0]**(2.0)).dot(X))
    robust_se_mat3 = np.transpose(X).dot(np.diag(residuals[:,0]**(2.0)/(1.0-np.diag(PX))).dot(X))

    v_robust = robust_se_mat1.dot(robust_se_mat2.dot(robust_se_mat1))
    v_modified_robust = robust_se_mat1.dot(robust_se_mat3.dot(robust_se_mat1))

    """ Returns:
        beta -- The vector of coefficient estimates, ordered on the columns on X.
        residuals -- The vector of residuals, Y - X.beta
        sig_hat -- The sample variance of the residuals.
        conv_st_error -- The 'standard error of the regression', sqrt(sig_hat).
        sig_hat_asymptotic_variance -- The analytic formula for the large sample variance
        ols_cov -- The covariance matrix under the basic OLS assumptions.
        v_robust -- The "robust" covariance matrix, weighted to account for the residuals and heteroskedasticity.
        v_modified_robust -- The bias-corrected and heteroskedasticity-consistent covariance matrix.
    """
    return beta, residuals, sig_hat, conv_st_err, sig_hat_asymptotic_variance, ols_cov, v_robust, v_modified_robust

あなたの問題では、次のように使用します。

import numpy as np

# Define or load your data:
(Y, X) = ....

# Desired polynomial degree
deg = 2;

reg_x_data =  np.ones(X.shape); # 0th degree term. 
for ii in range(1,deg+1):
    reg_x_data = np.hstack(( reg_x_data, X**(ii) )); # Append the ii^th degree term.


# Get all of the regression data.
beta, residuals, sig_hat, conv_st_error, sig_hat_asymptotic_variance, ols_cov, v_robust, v_modified_robust = ols_linreg(reg_x_data,Y)

# Print the covariance matrix:
print ols_cov

私の計算 (特に不等分散整合性推定器) にバグを見つけた場合は、お知らせください。できるだけ早く修正します。

于 2012-04-21T19:24:01.770 に答える