1""" 2Module of kernels that are able to handle continuous as well as categorical 3variables (both ordered and unordered). 4 5This is a slight deviation from the current approach in 6statsmodels.nonparametric.kernels where each kernel is a class object. 7 8Having kernel functions rather than classes makes extension to a multivariate 9kernel density estimation much easier. 10 11NOTE: As it is, this module does not interact with the existing API 12""" 13 14import numpy as np 15from scipy.special import erf 16 17 18#TODO: 19# - make sure we only receive int input for wang-ryzin and aitchison-aitken 20# - Check for the scalar Xi case everywhere 21 22 23def aitchison_aitken(h, Xi, x, num_levels=None): 24 r""" 25 The Aitchison-Aitken kernel, used for unordered discrete random variables. 26 27 Parameters 28 ---------- 29 h : 1-D ndarray, shape (K,) 30 The bandwidths used to estimate the value of the kernel function. 31 Xi : 2-D ndarray of ints, shape (nobs, K) 32 The value of the training set. 33 x : 1-D ndarray, shape (K,) 34 The value at which the kernel density is being estimated. 35 num_levels : bool, optional 36 Gives the user the option to specify the number of levels for the 37 random variable. If False, the number of levels is calculated from 38 the data. 39 40 Returns 41 ------- 42 kernel_value : ndarray, shape (nobs, K) 43 The value of the kernel function at each training point for each var. 44 45 Notes 46 ----- 47 See p.18 of [2]_ for details. The value of the kernel L if :math:`X_{i}=x` 48 is :math:`1-\lambda`, otherwise it is :math:`\frac{\lambda}{c-1}`. 49 Here :math:`c` is the number of levels plus one of the RV. 50 51 References 52 ---------- 53 .. [*] J. Aitchison and C.G.G. Aitken, "Multivariate binary discrimination 54 by the kernel method", Biometrika, vol. 63, pp. 413-420, 1976. 55 .. [*] Racine, Jeff. "Nonparametric Econometrics: A Primer," Foundation 56 and Trends in Econometrics: Vol 3: No 1, pp1-88., 2008. 57 """ 58 Xi = Xi.reshape(Xi.size) # seems needed in case Xi is scalar 59 if num_levels is None: 60 num_levels = np.asarray(np.unique(Xi).size) 61 62 kernel_value = np.ones(Xi.size) * h / (num_levels - 1) 63 idx = Xi == x 64 kernel_value[idx] = (idx * (1 - h))[idx] 65 return kernel_value 66 67 68def wang_ryzin(h, Xi, x): 69 r""" 70 The Wang-Ryzin kernel, used for ordered discrete random variables. 71 72 Parameters 73 ---------- 74 h : scalar or 1-D ndarray, shape (K,) 75 The bandwidths used to estimate the value of the kernel function. 76 Xi : ndarray of ints, shape (nobs, K) 77 The value of the training set. 78 x : scalar or 1-D ndarray of shape (K,) 79 The value at which the kernel density is being estimated. 80 81 Returns 82 ------- 83 kernel_value : ndarray, shape (nobs, K) 84 The value of the kernel function at each training point for each var. 85 86 Notes 87 ----- 88 See p. 19 in [1]_ for details. The value of the kernel L if 89 :math:`X_{i}=x` is :math:`1-\lambda`, otherwise it is 90 :math:`\frac{1-\lambda}{2}\lambda^{|X_{i}-x|}`, where :math:`\lambda` is 91 the bandwidth. 92 93 References 94 ---------- 95 .. [*] Racine, Jeff. "Nonparametric Econometrics: A Primer," Foundation 96 and Trends in Econometrics: Vol 3: No 1, pp1-88., 2008. 97 http://dx.doi.org/10.1561/0800000009 98 .. [*] M.-C. Wang and J. van Ryzin, "A class of smooth estimators for 99 discrete distributions", Biometrika, vol. 68, pp. 301-309, 1981. 100 """ 101 Xi = Xi.reshape(Xi.size) # seems needed in case Xi is scalar 102 kernel_value = 0.5 * (1 - h) * (h ** abs(Xi - x)) 103 idx = Xi == x 104 kernel_value[idx] = (idx * (1 - h))[idx] 105 return kernel_value 106 107 108def gaussian(h, Xi, x): 109 """ 110 Gaussian Kernel for continuous variables 111 Parameters 112 ---------- 113 h : 1-D ndarray, shape (K,) 114 The bandwidths used to estimate the value of the kernel function. 115 Xi : 1-D ndarray, shape (K,) 116 The value of the training set. 117 x : 1-D ndarray, shape (K,) 118 The value at which the kernel density is being estimated. 119 120 Returns 121 ------- 122 kernel_value : ndarray, shape (nobs, K) 123 The value of the kernel function at each training point for each var. 124 """ 125 return (1. / np.sqrt(2 * np.pi)) * np.exp(-(Xi - x)**2 / (h**2 * 2.)) 126 127 128def tricube(h, Xi, x): 129 """ 130 Tricube Kernel for continuous variables 131 Parameters 132 ---------- 133 h : 1-D ndarray, shape (K,) 134 The bandwidths used to estimate the value of the kernel function. 135 Xi : 1-D ndarray, shape (K,) 136 The value of the training set. 137 x : 1-D ndarray, shape (K,) 138 The value at which the kernel density is being estimated. 139 140 Returns 141 ------- 142 kernel_value : ndarray, shape (nobs, K) 143 The value of the kernel function at each training point for each var. 144 """ 145 u = (Xi - x) / h 146 u[np.abs(u) > 1] = 0 147 return (70. / 81) * (1 - np.abs(u)**3)**3 148 149 150def gaussian_convolution(h, Xi, x): 151 """ Calculates the Gaussian Convolution Kernel """ 152 return (1. / np.sqrt(4 * np.pi)) * np.exp(- (Xi - x)**2 / (h**2 * 4.)) 153 154 155def wang_ryzin_convolution(h, Xi, Xj): 156 # This is the equivalent of the convolution case with the Gaussian Kernel 157 # However it is not exactly convolution. Think of a better name 158 # References 159 ordered = np.zeros(Xi.size) 160 for x in np.unique(Xi): 161 ordered += wang_ryzin(h, Xi, x) * wang_ryzin(h, Xj, x) 162 163 return ordered 164 165 166def aitchison_aitken_convolution(h, Xi, Xj): 167 Xi_vals = np.unique(Xi) 168 ordered = np.zeros(Xi.size) 169 num_levels = Xi_vals.size 170 for x in Xi_vals: 171 ordered += aitchison_aitken(h, Xi, x, num_levels=num_levels) * \ 172 aitchison_aitken(h, Xj, x, num_levels=num_levels) 173 174 return ordered 175 176 177def gaussian_cdf(h, Xi, x): 178 return 0.5 * h * (1 + erf((x - Xi) / (h * np.sqrt(2)))) 179 180 181def aitchison_aitken_cdf(h, Xi, x_u): 182 x_u = int(x_u) 183 Xi_vals = np.unique(Xi) 184 ordered = np.zeros(Xi.size) 185 num_levels = Xi_vals.size 186 for x in Xi_vals: 187 if x <= x_u: #FIXME: why a comparison for unordered variables? 188 ordered += aitchison_aitken(h, Xi, x, num_levels=num_levels) 189 190 return ordered 191 192 193def wang_ryzin_cdf(h, Xi, x_u): 194 ordered = np.zeros(Xi.size) 195 for x in np.unique(Xi): 196 if x <= x_u: 197 ordered += wang_ryzin(h, Xi, x) 198 199 return ordered 200 201 202def d_gaussian(h, Xi, x): 203 # The derivative of the Gaussian Kernel 204 return 2 * (Xi - x) * gaussian(h, Xi, x) / h**2 205 206 207def aitchison_aitken_reg(h, Xi, x): 208 """ 209 A version for the Aitchison-Aitken kernel for nonparametric regression. 210 211 Suggested by Li and Racine. 212 """ 213 kernel_value = np.ones(Xi.size) 214 ix = Xi != x 215 inDom = ix * h 216 kernel_value[ix] = inDom[ix] 217 return kernel_value 218 219 220def wang_ryzin_reg(h, Xi, x): 221 """ 222 A version for the Wang-Ryzin kernel for nonparametric regression. 223 224 Suggested by Li and Racine in [1] ch.4 225 """ 226 return h ** abs(Xi - x) 227