Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/statsmodels/genmod/families/links.py : 46%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1'''
2Defines the link functions to be used with GLM and GEE families.
3'''
5import numpy as np
6import scipy.stats
7FLOAT_EPS = np.finfo(float).eps
10class Link(object):
11 """
12 A generic link function for one-parameter exponential family.
14 `Link` does nothing, but lays out the methods expected of any subclass.
15 """
17 def __call__(self, p):
18 """
19 Return the value of the link function. This is just a placeholder.
21 Parameters
22 ----------
23 p : array_like
24 Probabilities
26 Returns
27 -------
28 g(p) : array_like
29 The value of the link function g(p) = z
30 """
31 return NotImplementedError
33 def inverse(self, z):
34 """
35 Inverse of the link function. Just a placeholder.
37 Parameters
38 ----------
39 z : array_like
40 `z` is usually the linear predictor of the transformed variable
41 in the IRLS algorithm for GLM.
43 Returns
44 -------
45 g^(-1)(z) : ndarray
46 The value of the inverse of the link function g^(-1)(z) = p
47 """
48 return NotImplementedError
50 def deriv(self, p):
51 """
52 Derivative of the link function g'(p). Just a placeholder.
54 Parameters
55 ----------
56 p : array_like
58 Returns
59 -------
60 g'(p) : ndarray
61 The value of the derivative of the link function g'(p)
62 """
63 return NotImplementedError
65 def deriv2(self, p):
66 """Second derivative of the link function g''(p)
68 implemented through numerical differentiation
69 """
70 from statsmodels.tools.numdiff import approx_fprime_cs
71 # TODO: workaround proplem with numdiff for 1d
72 return np.diag(approx_fprime_cs(p, self.deriv))
74 def inverse_deriv(self, z):
75 """
76 Derivative of the inverse link function g^(-1)(z).
78 Parameters
79 ----------
80 z : array_like
81 `z` is usually the linear predictor for a GLM or GEE model.
83 Returns
84 -------
85 g'^(-1)(z) : ndarray
86 The value of the derivative of the inverse of the link function
88 Notes
89 -----
90 This reference implementation gives the correct result but is
91 inefficient, so it can be overridden in subclasses.
92 """
93 return 1 / self.deriv(self.inverse(z))
95 def inverse_deriv2(self, z):
96 """
97 Second derivative of the inverse link function g^(-1)(z).
99 Parameters
100 ----------
101 z : array_like
102 `z` is usually the linear predictor for a GLM or GEE model.
104 Returns
105 -------
106 g'^(-1)(z) : ndarray
107 The value of the second derivative of the inverse of the link
108 function
110 Notes
111 -----
112 This reference implementation gives the correct result but is
113 inefficient, so it can be overridden in subclasses.
114 """
115 iz = self.inverse(z)
116 return -self.deriv2(iz) / self.deriv(iz)**3
119class Logit(Link):
120 """
121 The logit transform
123 Notes
124 -----
125 call and derivative use a private method _clean to make trim p by
126 machine epsilon so that p is in (0,1)
128 Alias of Logit:
129 logit = Logit()
130 """
132 def _clean(self, p):
133 """
134 Clip logistic values to range (eps, 1-eps)
136 Parameters
137 ----------
138 p : array_like
139 Probabilities
141 Returns
142 -------
143 pclip : ndarray
144 Clipped probabilities
145 """
146 return np.clip(p, FLOAT_EPS, 1. - FLOAT_EPS)
148 def __call__(self, p):
149 """
150 The logit transform
152 Parameters
153 ----------
154 p : array_like
155 Probabilities
157 Returns
158 -------
159 z : ndarray
160 Logit transform of `p`
162 Notes
163 -----
164 g(p) = log(p / (1 - p))
165 """
166 p = self._clean(p)
167 return np.log(p / (1. - p))
169 def inverse(self, z):
170 """
171 Inverse of the logit transform
173 Parameters
174 ----------
175 z : array_like
176 The value of the logit transform at `p`
178 Returns
179 -------
180 p : ndarray
181 Probabilities
183 Notes
184 -----
185 g^(-1)(z) = exp(z)/(1+exp(z))
186 """
187 z = np.asarray(z)
188 t = np.exp(-z)
189 return 1. / (1. + t)
191 def deriv(self, p):
192 """
193 Derivative of the logit transform
195 Parameters
196 ----------
197 p: array_like
198 Probabilities
200 Returns
201 -------
202 g'(p) : ndarray
203 Value of the derivative of logit transform at `p`
205 Notes
206 -----
207 g'(p) = 1 / (p * (1 - p))
209 Alias for `Logit`:
210 logit = Logit()
211 """
212 p = self._clean(p)
213 return 1. / (p * (1 - p))
215 def inverse_deriv(self, z):
216 """
217 Derivative of the inverse of the logit transform
219 Parameters
220 ----------
221 z : array_like
222 `z` is usually the linear predictor for a GLM or GEE model.
224 Returns
225 -------
226 g'^(-1)(z) : ndarray
227 The value of the derivative of the inverse of the logit function
228 """
229 t = np.exp(z)
230 return t/(1 + t)**2
232 def deriv2(self, p):
233 """
234 Second derivative of the logit function.
236 Parameters
237 ----------
238 p : array_like
239 probabilities
241 Returns
242 -------
243 g''(z) : ndarray
244 The value of the second derivative of the logit function
245 """
246 v = p * (1 - p)
247 return (2*p - 1) / v**2
250class logit(Logit):
251 pass
254class Power(Link):
255 """
256 The power transform
258 Parameters
259 ----------
260 power : float
261 The exponent of the power transform
263 Notes
264 -----
265 Aliases of Power:
266 inverse = Power(power=-1)
267 sqrt = Power(power=.5)
268 inverse_squared = Power(power=-2.)
269 identity = Power(power=1.)
270 """
272 def __init__(self, power=1.):
273 self.power = power
275 def __call__(self, p):
276 """
277 Power transform link function
279 Parameters
280 ----------
281 p : array_like
282 Mean parameters
284 Returns
285 -------
286 z : array_like
287 Power transform of x
289 Notes
290 -----
291 g(p) = x**self.power
292 """
293 if self.power == 1:
294 return p
295 else:
296 return np.power(p, self.power)
298 def inverse(self, z):
299 """
300 Inverse of the power transform link function
302 Parameters
303 ----------
304 `z` : array_like
305 Value of the transformed mean parameters at `p`
307 Returns
308 -------
309 `p` : ndarray
310 Mean parameters
312 Notes
313 -----
314 g^(-1)(z`) = `z`**(1/`power`)
315 """
316 if self.power == 1:
317 return z
318 else:
319 return np.power(z, 1. / self.power)
321 def deriv(self, p):
322 """
323 Derivative of the power transform
325 Parameters
326 ----------
327 p : array_like
328 Mean parameters
330 Returns
331 -------
332 g'(p) : ndarray
333 Derivative of power transform of `p`
335 Notes
336 -----
337 g'(`p`) = `power` * `p`**(`power` - 1)
338 """
339 if self.power == 1:
340 return np.ones_like(p)
341 else:
342 return self.power * np.power(p, self.power - 1)
344 def deriv2(self, p):
345 """
346 Second derivative of the power transform
348 Parameters
349 ----------
350 p : array_like
351 Mean parameters
353 Returns
354 -------
355 g''(p) : ndarray
356 Second derivative of the power transform of `p`
358 Notes
359 -----
360 g''(`p`) = `power` * (`power` - 1) * `p`**(`power` - 2)
361 """
362 if self.power == 1:
363 return np.zeros_like(p)
364 else:
365 return self.power * (self.power - 1) * np.power(p, self.power - 2)
367 def inverse_deriv(self, z):
368 """
369 Derivative of the inverse of the power transform
371 Parameters
372 ----------
373 z : array_like
374 `z` is usually the linear predictor for a GLM or GEE model.
376 Returns
377 -------
378 g^(-1)'(z) : ndarray
379 The value of the derivative of the inverse of the power transform
380 function
381 """
382 if self.power == 1:
383 return np.ones_like(z)
384 else:
385 return np.power(z, (1 - self.power)/self.power) / self.power
387 def inverse_deriv2(self, z):
388 """
389 Second derivative of the inverse of the power transform
391 Parameters
392 ----------
393 z : array_like
394 `z` is usually the linear predictor for a GLM or GEE model.
396 Returns
397 -------
398 g^(-1)'(z) : ndarray
399 The value of the derivative of the inverse of the power transform
400 function
401 """
402 if self.power == 1:
403 return np.zeros_like(z)
404 else:
405 return ((1 - self.power) *
406 np.power(z, (1 - 2*self.power)/self.power) / self.power**2)
409class inverse_power(Power):
410 """
411 The inverse transform
413 Notes
414 -----
415 g(p) = 1/p
417 Alias of statsmodels.family.links.Power(power=-1.)
418 """
419 def __init__(self):
420 super(inverse_power, self).__init__(power=-1.)
423class sqrt(Power):
424 """
425 The square-root transform
427 Notes
428 -----
429 g(`p`) = sqrt(`p`)
431 Alias of statsmodels.family.links.Power(power=.5)
432 """
433 def __init__(self):
434 super(sqrt, self).__init__(power=.5)
437class inverse_squared(Power):
438 r"""
439 The inverse squared transform
441 Notes
442 -----
443 g(`p`) = 1/(`p`\*\*2)
445 Alias of statsmodels.family.links.Power(power=2.)
446 """
447 def __init__(self):
448 super(inverse_squared, self).__init__(power=-2.)
451class identity(Power):
452 """
453 The identity transform
455 Notes
456 -----
457 g(`p`) = `p`
459 Alias of statsmodels.family.links.Power(power=1.)
460 """
461 def __init__(self):
462 super(identity, self).__init__(power=1.)
465class Log(Link):
466 """
467 The log transform
469 Notes
470 -----
471 call and derivative call a private method _clean to trim the data by
472 machine epsilon so that p is in (0,1). log is an alias of Log.
473 """
475 def _clean(self, x):
476 return np.clip(x, FLOAT_EPS, np.inf)
478 def __call__(self, p, **extra):
479 """
480 Log transform link function
482 Parameters
483 ----------
484 x : array_like
485 Mean parameters
487 Returns
488 -------
489 z : ndarray
490 log(x)
492 Notes
493 -----
494 g(p) = log(p)
495 """
496 x = self._clean(p)
497 return np.log(x)
499 def inverse(self, z):
500 """
501 Inverse of log transform link function
503 Parameters
504 ----------
505 z : ndarray
506 The inverse of the link function at `p`
508 Returns
509 -------
510 p : ndarray
511 The mean probabilities given the value of the inverse `z`
513 Notes
514 -----
515 g^{-1}(z) = exp(z)
516 """
517 return np.exp(z)
519 def deriv(self, p):
520 """
521 Derivative of log transform link function
523 Parameters
524 ----------
525 p : array_like
526 Mean parameters
528 Returns
529 -------
530 g'(p) : ndarray
531 derivative of log transform of x
533 Notes
534 -----
535 g'(x) = 1/x
536 """
537 p = self._clean(p)
538 return 1. / p
540 def deriv2(self, p):
541 """
542 Second derivative of the log transform link function
544 Parameters
545 ----------
546 p : array_like
547 Mean parameters
549 Returns
550 -------
551 g''(p) : ndarray
552 Second derivative of log transform of x
554 Notes
555 -----
556 g''(x) = -1/x^2
557 """
558 p = self._clean(p)
559 return -1. / p**2
561 def inverse_deriv(self, z):
562 """
563 Derivative of the inverse of the log transform link function
565 Parameters
566 ----------
567 z : ndarray
568 The inverse of the link function at `p`
570 Returns
571 -------
572 g^(-1)'(z) : ndarray
573 The value of the derivative of the inverse of the log function,
574 the exponential function
575 """
576 return np.exp(z)
579class log(Log):
580 """
581 The log transform
583 Notes
584 -----
585 log is a an alias of Log.
586 """
587 pass
590# TODO: the CDFLink is untested
591class CDFLink(Logit):
592 """
593 The use the CDF of a scipy.stats distribution
595 CDFLink is a subclass of logit in order to use its _clean method
596 for the link and its derivative.
598 Parameters
599 ----------
600 dbn : scipy.stats distribution
601 Default is dbn=scipy.stats.norm
603 Notes
604 -----
605 The CDF link is untested.
606 """
608 def __init__(self, dbn=scipy.stats.norm):
609 self.dbn = dbn
611 def __call__(self, p):
612 """
613 CDF link function
615 Parameters
616 ----------
617 p : array_like
618 Mean parameters
620 Returns
621 -------
622 z : ndarray
623 (ppf) inverse of CDF transform of p
625 Notes
626 -----
627 g(`p`) = `dbn`.ppf(`p`)
628 """
629 p = self._clean(p)
630 return self.dbn.ppf(p)
632 def inverse(self, z):
633 """
634 The inverse of the CDF link
636 Parameters
637 ----------
638 z : array_like
639 The value of the inverse of the link function at `p`
641 Returns
642 -------
643 p : ndarray
644 Mean probabilities. The value of the inverse of CDF link of `z`
646 Notes
647 -----
648 g^(-1)(`z`) = `dbn`.cdf(`z`)
649 """
650 return self.dbn.cdf(z)
652 def deriv(self, p):
653 """
654 Derivative of CDF link
656 Parameters
657 ----------
658 p : array_like
659 mean parameters
661 Returns
662 -------
663 g'(p) : ndarray
664 The derivative of CDF transform at `p`
666 Notes
667 -----
668 g'(`p`) = 1./ `dbn`.pdf(`dbn`.ppf(`p`))
669 """
670 p = self._clean(p)
671 return 1. / self.dbn.pdf(self.dbn.ppf(p))
673 def deriv2(self, p):
674 """
675 Second derivative of the link function g''(p)
677 implemented through numerical differentiation
678 """
679 from statsmodels.tools.numdiff import approx_fprime
680 p = np.atleast_1d(p)
681 # Note: special function for norm.ppf does not support complex
682 return np.diag(approx_fprime(p, self.deriv, centered=True))
684 def inverse_deriv(self, z):
685 """
686 Derivative of the inverse of the CDF transformation link function
688 Parameters
689 ----------
690 z : ndarray
691 The inverse of the link function at `p`
693 Returns
694 -------
695 g^(-1)'(z) : ndarray
696 The value of the derivative of the inverse of the logit function
697 """
698 return 1/self.deriv(self.inverse(z))
701class probit(CDFLink):
702 """
703 The probit (standard normal CDF) transform
705 Notes
706 -----
707 g(p) = scipy.stats.norm.ppf(p)
709 probit is an alias of CDFLink.
710 """
711 pass
714class cauchy(CDFLink):
715 """
716 The Cauchy (standard Cauchy CDF) transform
718 Notes
719 -----
720 g(p) = scipy.stats.cauchy.ppf(p)
722 cauchy is an alias of CDFLink with dbn=scipy.stats.cauchy
723 """
725 def __init__(self):
726 super(cauchy, self).__init__(dbn=scipy.stats.cauchy)
728 def deriv2(self, p):
729 """
730 Second derivative of the Cauchy link function.
732 Parameters
733 ----------
734 p: array_like
735 Probabilities
737 Returns
738 -------
739 g''(p) : ndarray
740 Value of the second derivative of Cauchy link function at `p`
741 """
742 a = np.pi * (p - 0.5)
743 d2 = 2 * np.pi**2 * np.sin(a) / np.cos(a)**3
744 return d2
747class CLogLog(Logit):
748 """
749 The complementary log-log transform
751 CLogLog inherits from Logit in order to have access to its _clean method
752 for the link and its derivative.
754 Notes
755 -----
756 CLogLog is untested.
757 """
758 def __call__(self, p):
759 """
760 C-Log-Log transform link function
762 Parameters
763 ----------
764 p : ndarray
765 Mean parameters
767 Returns
768 -------
769 z : ndarray
770 The CLogLog transform of `p`
772 Notes
773 -----
774 g(p) = log(-log(1-p))
775 """
776 p = self._clean(p)
777 return np.log(-np.log(1 - p))
779 def inverse(self, z):
780 """
781 Inverse of C-Log-Log transform link function
784 Parameters
785 ----------
786 z : array_like
787 The value of the inverse of the CLogLog link function at `p`
789 Returns
790 -------
791 p : ndarray
792 Mean parameters
794 Notes
795 -----
796 g^(-1)(`z`) = 1-exp(-exp(`z`))
797 """
798 return 1 - np.exp(-np.exp(z))
800 def deriv(self, p):
801 """
802 Derivative of C-Log-Log transform link function
804 Parameters
805 ----------
806 p : array_like
807 Mean parameters
809 Returns
810 -------
811 g'(p) : ndarray
812 The derivative of the CLogLog transform link function
814 Notes
815 -----
816 g'(p) = - 1 / ((p-1)*log(1-p))
817 """
818 p = self._clean(p)
819 return 1. / ((p - 1) * (np.log(1 - p)))
821 def deriv2(self, p):
822 """
823 Second derivative of the C-Log-Log ink function
825 Parameters
826 ----------
827 p : array_like
828 Mean parameters
830 Returns
831 -------
832 g''(p) : ndarray
833 The second derivative of the CLogLog link function
834 """
835 p = self._clean(p)
836 fl = np.log(1 - p)
837 d2 = -1 / ((1 - p)**2 * fl)
838 d2 *= 1 + 1 / fl
839 return d2
841 def inverse_deriv(self, z):
842 """
843 Derivative of the inverse of the C-Log-Log transform link function
845 Parameters
846 ----------
847 z : array_like
848 The value of the inverse of the CLogLog link function at `p`
850 Returns
851 -------
852 g^(-1)'(z) : ndarray
853 The derivative of the inverse of the CLogLog link function
854 """
855 return np.exp(z - np.exp(z))
858class cloglog(CLogLog):
859 """
860 The CLogLog transform link function.
862 Notes
863 -----
864 g(`p`) = log(-log(1-`p`))
866 cloglog is an alias for CLogLog
867 cloglog = CLogLog()
868 """
869 pass
872class NegativeBinomial(Link):
873 '''
874 The negative binomial link function
876 Parameters
877 ----------
878 alpha : float, optional
879 Alpha is the ancillary parameter of the Negative Binomial link
880 function. It is assumed to be nonstochastic. The default value is 1.
881 Permissible values are usually assumed to be in (.01, 2).
882 '''
884 def __init__(self, alpha=1.):
885 self.alpha = alpha
887 def _clean(self, x):
888 return np.clip(x, FLOAT_EPS, np.inf)
890 def __call__(self, p):
891 '''
892 Negative Binomial transform link function
894 Parameters
895 ----------
896 p : array_like
897 Mean parameters
899 Returns
900 -------
901 z : ndarray
902 The negative binomial transform of `p`
904 Notes
905 -----
906 g(p) = log(p/(p + 1/alpha))
907 '''
908 p = self._clean(p)
909 return np.log(p/(p + 1/self.alpha))
911 def inverse(self, z):
912 '''
913 Inverse of the negative binomial transform
915 Parameters
916 ----------
917 z : array_like
918 The value of the inverse of the negative binomial link at `p`.
920 Returns
921 -------
922 p : ndarray
923 Mean parameters
925 Notes
926 -----
927 g^(-1)(z) = exp(z)/(alpha*(1-exp(z)))
928 '''
929 return -1/(self.alpha * (1 - np.exp(-z)))
931 def deriv(self, p):
932 '''
933 Derivative of the negative binomial transform
935 Parameters
936 ----------
937 p : array_like
938 Mean parameters
940 Returns
941 -------
942 g'(p) : ndarray
943 The derivative of the negative binomial transform link function
945 Notes
946 -----
947 g'(x) = 1/(x+alpha*x^2)
948 '''
949 return 1/(p + self.alpha * p**2)
951 def deriv2(self, p):
952 '''
953 Second derivative of the negative binomial link function.
955 Parameters
956 ----------
957 p : array_like
958 Mean parameters
960 Returns
961 -------
962 g''(p) : ndarray
963 The second derivative of the negative binomial transform link
964 function
966 Notes
967 -----
968 g''(x) = -(1+2*alpha*x)/(x+alpha*x^2)^2
969 '''
970 numer = -(1 + 2 * self.alpha * p)
971 denom = (p + self.alpha * p**2)**2
972 return numer / denom
974 def inverse_deriv(self, z):
975 '''
976 Derivative of the inverse of the negative binomial transform
978 Parameters
979 ----------
980 z : array_like
981 Usually the linear predictor for a GLM or GEE model
983 Returns
984 -------
985 g^(-1)'(z) : ndarray
986 The value of the derivative of the inverse of the negative
987 binomial link
988 '''
989 t = np.exp(z)
990 return t / (self.alpha * (1-t)**2)
993class nbinom(NegativeBinomial):
994 """
995 The negative binomial link function.
997 Notes
998 -----
999 g(p) = log(p/(p + 1/alpha))
1001 nbinom is an alias of NegativeBinomial.
1002 nbinom = NegativeBinomial(alpha=1.)
1003 """
1004 pass