Hồi quy nhị thức âm Mô hình thống kê Python

Các ước tính tham số MLE chung của statsmodels và R thống nhất đến chữ số thập phân thứ tư. Tuy nhiên, các lỗi tiêu chuẩn chỉ đồng ý đến chữ số thập phân thứ hai. Sự khác biệt này là kết quả của sự thiếu chính xác trong ước tính số Hessian của chúng tôi. Trong ngữ cảnh hiện tại, sự khác biệt giữa ước tính sai số chuẩn của MASS và statsmodels về cơ bản là không liên quan, nhưng nó làm nổi bật thực tế là những người dùng cần ước tính rất chính xác có thể không phải lúc nào cũng muốn dựa vào cài đặt mặc định khi sử dụng đạo hàm số. Trong những trường hợp như vậy, tốt hơn là sử dụng các công cụ phái sinh phân tích với lớp LikelihoodModel

Trong ví dụ này, chúng tôi muốn sử dụng AlgoPy để giúp tính toán ước tính khả năng tối đa và sai số chuẩn của các tham số của mô hình phi tuyến tính. Điều này tuân theo ví dụ về khả năng tối đa chung của mô hình thống kê sử dụng bộ dữ liệu medpar

\[\begin{split}\mathcal{L}[\beta_j; y, \alpha] = \sum_{i=1}^n y_i \log & \left[ \frac{\alpha \exp[X'_i \

Đây là mã trăn

"""
Negative Binomial Regression

This is an algopy implementation of the statsmodels example
//statsmodels.sourceforge.net/devel/examples/generated/example_gmle.html
.
"""

import functools

import numpy
import scipy.optimize
import algopy
import numdifftools
import pandas
import patsy

g_url = '//vincentarelbundock.github.com/Rdatasets/csv/COUNT/medpar.csv'

def get_aic[y, X, theta]:
    return 2*len[theta] + 2*get_neg_ll[y, X, theta]

def get_neg_ll[y, X, theta]:
    alpha = theta[-1]
    beta = theta[:-1]
    a = alpha * algopy.exp[algopy.dot[X, beta]]
    ll = algopy.sum[
        -y*algopy.log1p[1/a] +
        -algopy.log1p[a] / alpha +
        algopy.special.gammaln[y + 1/alpha] +
        -algopy.special.gammaln[y + 1] +
        -algopy.special.gammaln[1/alpha]]
    neg_ll = -ll
    return neg_ll

def eval_grad[f, theta]:
    theta = algopy.UTPM.init_jacobian[theta]
    return algopy.UTPM.extract_jacobian[f[theta]]

def eval_hess[f, theta]:
    theta = algopy.UTPM.init_hessian[theta]
    return algopy.UTPM.extract_hessian[len[theta], f[theta]]

def main[]:

    # read the data from the internet into numpy arrays
    medpar = pandas.read_csv[g_url]
    y_patsy, X_patsy = patsy.dmatrices['los~type2+type3+hmo+white', medpar]
    y = numpy.array[y_patsy].flatten[]
    X = numpy.array[X_patsy]

    # define the objective function and the autodiff gradient and hessian
    f = functools.partial[get_neg_ll, y, X]
    g = functools.partial[eval_grad, f]
    h = functools.partial[eval_hess, f]

    # init the search for max likelihood parameters
    theta0 = numpy.array[[
        numpy.log[numpy.mean[y]],
        0, 0, 0, 0,
        0.5,
        ], dtype=float]

    # do the max likelihood search
    results = scipy.optimize.fmin_ncg[
            f,
            theta0,
            fprime=g,
            fhess=h,
            avextol=1e-6,
            ]

    # compute the hessian a couple of different ways
    algopy_hessian = h[results]
    num_hessian = numdifftools.Hessian[f][results]

    # report the results of the search including aic and standard error
    print['search results:']
    print[results]
    print[]
    print['aic:']
    print[get_aic[y, X, results]]
    print[]
    print['standard error using observed fisher information,']
    print['with hessian computed using algopy:']
    print[numpy.sqrt[numpy.diag[scipy.linalg.inv[algopy_hessian]]]]
    print[]
    print['standard error using observed fisher information,']
    print['with hessian computed using numdifftools:']
    print[numpy.sqrt[numpy.diag[scipy.linalg.inv[num_hessian]]]]
    print[]


if __name__ == '__main__':
    main[]

Đây là đầu ra của nó

Optimization terminated successfully.
         Current function value: 4797.476603
         Iterations: 10
         Function evaluations: 11
         Gradient evaluations: 10
         Hessian evaluations: 10
search results:
[ 2.31027893  0.22124897  0.70615882 -0.06795522 -0.12906544  0.44575671]

aic:
9606.95320507

standard error using observed fisher information,
with hessian computed using algopy:
[ 0.06794736  0.05059255  0.07613111  0.05326133  0.06854179  0.01981577]

standard error using observed fisher information,
with hessian computed using numdifftools:
[ 0.06794736  0.05059255  0.07613111  0.05326133  0.06854179  0.01981577]

Sự thống nhất giữa kết quả đầu ra này và các kết quả ví dụ về mô hình thống kê cho thấy rằng các mô hình thống kê báo trước về độ chính xác bằng số có thể không phải là kết quả của các vấn đề về số với đạo hàm, mà thay vào đó, việc triển khai R MASS có thể đưa ra các ước tính sai số chuẩn kém chính xác hơn hoặc có thể không sử dụng