diff --git a/mle_estimation.py b/mle_estimation.py deleted file mode 100644 index ee6e87eb000e4718bb8a7235ced6d65a8f5ddd16..0000000000000000000000000000000000000000 --- a/mle_estimation.py +++ /dev/null @@ -1,64 +0,0 @@ -''' -Created on Nov 16, 2015 -Class which does MLE estimation for perfect Data. -For real Binning see POPRES-Analysis -@author: Harald -''' -l0 = 0.05 # Minimum block length which is reported - - -import numpy as np -import math -from statsmodels.base.model import GenericLikelihoodModel -from scipy.special import kv as kv - - -def single_pair(l_vec, r, C, sigma): - '''Estimates the likelihood for a single pair. Assert len l_vec>0''' - return np.sum(np.log([C * r ** 2 / (2 * l * sigma ** 2) * kv(2, np.sqrt(2 * l) * r / sigma) for l in l_vec])) - -def pairwise_ll(l, r, C, sigma): - '''Full Pairwise Likelihood function for data.''' - print("C: %.5f" % C) - print("Sigma: %.4f" % sigma) - - if C <= 0 or sigma <= 0: - return np.zeros_like(l) # If Parameters do not make sense. - else: - pr_noshare = -C * r / (np.sqrt(2 * l0) * sigma) * kv(1, np.sqrt(2 * l0) * r / sigma) # Standard vector of no-sharing - f_share = np.array([single_pair(l[i], r[i], C, sigma) if l[i] != 0 else 0.0 for i in range(0, len(r))]) # For vector send it to single_pair - res = pr_noshare[:, 0] + f_share - return res.astype(np.float) - -def pairwise_ll01(l, r, C, sigma): - '''Pairwise Likelihood function only caring about block or not''' - print("C: %.5f" % C) - print("Sigma: %.4f" % sigma) - if C <= 0 or sigma <= 0: - return np.zeros_like(l) # If Parameters do not make sense. - else: - lambd = (C * r / (np.sqrt(2 * l0) * sigma) * kv(1, np.sqrt(2 * l0) * r / sigma)) # Probability of sharing, vectorized. - pr_noshare = np.exp(-lambd) # Probabilities of no sharing - l = [len(l[i]) if l[i] != 0 else 0 for i in range(0, len(r))] # Number of shared blocks - pr_share = np.array([(lambd[i] ** l[i]) / math.factorial(l[i]) if l[i] != 0 else 1 for i in range(0, len(r))]) - res = pr_noshare[:, 0] * pr_share # Bring together the two terms - return res.astype(np.float) - -class MLE_estimation(GenericLikelihoodModel): - def __init__(self, endog, exog, **kwds): - super(MLE_estimation, self).__init__(endog, exog, **kwds) - - def nloglikeobs(self, params): - C = params[0] - sigma = params[1] - p_ll = pairwise_ll(self.endog, self.exog, C, sigma) - nll = -p_ll # First is length of shared block, second is pairwise distance - return nll - - def fit(self, start_params=None, maxiter=10000, maxfun=5000, **kwds): - # we have one additional parameter and we need to add it for summary - if start_params == None: - start_params = np.array([0.02, 2]) - return super(MLE_estimation, self).fit(start_params=start_params, - maxiter=maxiter, maxfun=maxfun, - **kwds)