как получить асимтотическую ковариационную матрицу 2-параметрического распределения Вейбулла в питоне - PullRequest
Новая
       51

как получить асимтотическую ковариационную матрицу 2-параметрического распределения Вейбулла в питоне

0 голосов
/ 26 июня 2018

Я пытаюсь реализовать код Matlab с помощью Python.

В Matlab доверительный интервал распределения Вейбулла получается следующим образом: 

CI = 0.95 % 
p = [scale_hat shape_hat];
[nlogl,avar]=wbllike(p,y_ordered);
[qmid,qlo,qup]=wblinv(1-1/nCount,p(1),p(2),avar,1-CI);

В Scipy я не смог найтиэквивалент функции wbllike.Похоже, функция не реализована.Вместо этого я попытался реализовать функцию wblike самостоятельно.до сих пор отрицательное логарифмическое правдоподобие можно было легко получить следующим образом: 

def weibull_neg_loglike(x, shape, scale):
    return -stats.weibull_min.logpdf(x, shape_param, scale=scale_param).sum()

Я обнаружил, что в python есть пакет statsmodels, который предоставляет GeneralLikelyhoodModel.но не разбирайтесь, как заставить вещи работать.Я немного изучил сайт, http://rlhick.people.wm.edu/posts/estimating-custom-mle.html#statsmodels, но все же тщетно.

Любой совет будет полезен.

1 Ответ

0 голосов
/ 26 июня 2018

Ниже приведена часть записной книжки, которую я написал некоторое время назад.Он предназначен для цензурной регрессии Вейбулла, поэтому для него требуется явно указанный экзог, например, массив из них.

Он также включает обходной путь для logsf, который не работал в моей старой версии scipy. 

    def weibull_min_logsf(x, c, scale=1):
        x = x / scale
        return -np.power(x, c)

    from statsmodels.base.model import GenericLikelihoodModel

    class WeibullModel(GenericLikelihoodModel):

        def __init__(self, endog, exog, censored=None, **kwds):
            super(WeibullModel, self).__init__(endog, exog, **kwds)
            if censored is not None:
                self.censored = censored.astype(int)
            else:
                self.censored = np.zeros(len(self.endog), np.int)

            self.k_params = self.exog.shape[1] + 1


        def loglike(self, params):
            params_ex = params[:-1]
            params_shape = params[-1]
            m = self.exog.dot(params_ex)
            llf = (1 - self.censored) * stats.weibull_min.logpdf(self.endog, params_shape, scale=m)
            #stats.weibull_min.logsf overflows in my older version of scipy, changed in newer versions
            #llf2 = self.censored * stats.weibull_min.logsf(self.endog, params_shape, scale=m)
            llf += self.censored * weibull_min_logsf(self.endog, params_shape, scale=m)
            return llf.sum()

        def _get_distribution(self, params, exog):
            """similar to a predict method
            """
            params_ex = params[:-1]
            params_shape = params[-1]
            m = exog.dot(params_ex)
            return stats.weibull_min(params_shape, scale=m)

    mod = WeibullModel(data['duration'], np.ones(len(data)))
    res = mod.fit(start_params=np.ones(2)) 
    print(res.summary())
...