# Work to undnerstand how to use the kappa distribution in pymc3

``````

In [44]:

import itertools

import matplotlib.pyplot as plt
import matplotlib as mpl
from pymc3 import Model, Normal, Slice
from pymc3 import sample
from pymc3 import traceplot
from pymc3.distributions import Interpolated
from theano import as_op
import theano.tensor as tt
import numpy as np
from scipy import stats
import tqdm

%matplotlib inline

%version_information pymc3, scipy

``````
``````

Out[44]:

SoftwareVersionPython3.6.2 64bit [GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)]IPython6.1.0OSDarwin 15.6.0 x86_64 i386 64bitpymc33.1scipy0.19.1Mon Sep 18 12:13:55 2017 MDT

``````
``````

In [3]:

stats.kappa3?

``````
``````

In [14]:

k3_rvs = stats.kappa3.rvs(a=1, loc=90, scale=45, size=10000)
plt.hist(k3_rvs, 50);

``````
``````

``````
``````

In [111]:

# stats.kappa4?
h=100
k=0
k4_pdf = stats.kappa4.rvs(h=k, k=k, loc=90, scale=20, size=10000)
plt.hist(k4_pdf, 50);
mean, var, skew, kurt = stats.kappa4.stats(h, k, moments='mvsk')
print(mean, var, skew, kurt)

``````
``````

4.621500393020947 0.023720538535741298 17.440372233406986 431.65646930084574

``````
``````

In [24]:

import scipy.optimize

``````
``````

In [92]:

def moments(x):
h, k = x
mean, var, skew, kurt = stats.kappa4.stats(h, k, moments='mvsk')
return skew

minme = lambda x: np.abs(moments(x))

x0 = [.2, 4]

minme(x0)

``````
``````

Out[92]:

7.5529025101903082

``````
``````

In [97]:

# optimization seems slow, lets try a brute force to get a feel
N=7
h = np.linspace(0, 5, N)
k = np.linspace(0, 5, N)

H, K = np.meshgrid(h, k)

ans = np.zeros((N,N), dtype=float)

``````
``````

In [98]:

for i, j in tqdm.tqdm(itertools.product(range(N), range(N)), total=N*N):
ans[i,j] = minme((H[i,j], K[i,j]))

``````
``````

86%|████████▌ | 42/49 [00:52<00:02,  2.66it/s]/Users/balarsen/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/integrate/quadpack.py:364: IntegrationWarning: The algorithm does not converge.  Roundoff error is detected
in the extrapolation table.  It is assumed that the requested tolerance
cannot be achieved, and that the returned result (if full_output = 1) is
the best which can be obtained.
warnings.warn(msg, IntegrationWarning)
100%|██████████| 49/49 [00:58<00:00,  2.84it/s]

``````
``````

In [100]:

plt.contourf(H,K,ans)
plt.colorbar()
plt.contour(H,K,ans, levels=[0,.1,.2], colors='w')

print(np.unravel_index(np.argmin(ans), ans.shape))
print(ans.min())
min = np.unravel_index(np.argmin(ans), ans.shape)
plt.scatter(H[min], K[min], marker='o', c='r')
print(H[min], K[min])

``````
``````

(1, 1)
0.00946777071186
0.833333333333 0.833333333333

``````
``````

In [91]:

niter = 0
def callbackF(x):
global niter
print("{} : {}, {}".format(niter, x[0], x[1]))
niter += 1

scipy.optimize.minimize(minme, x0, method='TNC', callback=callbackF, bounds=[[0,1], [0,1]])

``````
``````

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-91-6a675168301b> in <module>()
6
7
----> 8 scipy.optimize.minimize(minme, x0, method='TNC', callback=callbackF, bounds=[[0,1], [0,1]])

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/optimize/_minimize.py in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
451     elif meth == 'tnc':
452         return _minimize_tnc(fun, x0, args, jac, bounds, callback=callback,
--> 453                              **options)
454     elif meth == 'cobyla':
455         return _minimize_cobyla(fun, x0, args, constraints, **options)

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/optimize/tnc.py in _minimize_tnc(fun, x0, args, jac, bounds, eps, scale, offset, mesg_num, maxCGit, maxiter, eta, stepmx, accuracy, minfev, ftol, xtol, gtol, rescale, disp, callback, **unknown_options)
407                                         offset, messages, maxCGit, maxfun,
408                                         eta, stepmx, accuracy, fmin, ftol,
--> 409                                         xtol, pgtol, rescale, callback)
410

366             f = fun(x, *args)
--> 367             g = approx_fprime(x, fun, epsilon, *args)
368             return f, g
369     else:

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/optimize/optimize.py in approx_fprime(xk, f, epsilon, *args)
686
687     """
--> 688     return _approx_fprime_helper(xk, f, epsilon, args=args)
689
690

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/optimize/optimize.py in _approx_fprime_helper(xk, f, epsilon, args, f0)
626         ei[k] = 1.0
627         d = epsilon * ei
--> 628         grad[k] = (f(*((xk + d,) + args)) - f0) / d[k]
629         ei[k] = 0.0

<ipython-input-31-62019fb4b4c9> in <lambda>(x)
4     return skew
5
----> 6 minme = lambda x: np.abs(moments(x))
7
8 x0 = [.2, 4]

<ipython-input-31-62019fb4b4c9> in moments(x)
1 def moments(x):
2     h, k = x
----> 3     mean, var, skew, kurt = stats.kappa4.stats(h, k, moments='mvsk')
4     return skew
5

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py in stats(self, *args, **kwds)
1025             if 'm' in moments:
1026                 if mu is None:
-> 1027                     mu = self._munp(1, *goodargs)
1028                 out0 = default.copy()
1029                 place(out0, cond, mu * scale + loc)

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py in _munp(self, n, *args)
783         # Silence floating point warnings from integration.
784         olderr = np.seterr(all='ignore')
--> 785         vals = self.generic_moment(n, *args)
786         np.seterr(**olderr)
787         return vals

~/miniconda3/envs/python3/lib/python3.6/site-packages/numpy/lib/function_base.py in __call__(self, *args, **kwargs)
2732             vargs.extend([kwargs[_n] for _n in names])
2733
-> 2734         return self._vectorize_call(func=func, args=vargs)
2735
2736     def _get_ufunc_and_otypes(self, func, args):

~/miniconda3/envs/python3/lib/python3.6/site-packages/numpy/lib/function_base.py in _vectorize_call(self, func, args)
2808                       for a in args]
2809
-> 2810             outputs = ufunc(*inputs)
2811
2812             if ufunc.nout == 1:

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py in _mom1_sc(self, m, *args)
1603
1604     def _mom1_sc(self, m, *args):
-> 1605         return integrate.quad(self._mom_integ1, 0, 1, args=(m,)+args)[0]
1606
1607     def _pdf(self, x, *args):

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/integrate/quadpack.py in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
321     if (weight is None):
322         retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
--> 323                        points)
324     else:
325         retval = _quad_weight(func, a, b, args, full_output, epsabs, epsrel,

386     if points is None:
387         if infbounds == 0:
389         else:

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py in _mom_integ1(self, q, m, *args)
1600     # moment calculated using ppf
1601     def _mom_integ1(self, q, m, *args):
-> 1602         return (self.ppf(q, *args))**m
1603
1604     def _mom1_sc(self, m, *args):

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py in ppf(self, q, *args, **kwds)
1898         q, loc, scale = map(asarray, (q, loc, scale))
1899         args = tuple(map(asarray, args))
-> 1900         cond0 = self._argcheck(*args) & (scale > 0) & (loc == loc)
1901         cond1 = (0 < q) & (q < 1)
1902         cond2 = cond0 & (q == 0)

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/stats/_continuous_distns.py in _argcheck(self, h, k)
3496                              [f0, f1, f1, f0, f1, f1],
3497                              [h, k],
-> 3498                              default=np.nan)
3499         return h == h
3500

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/_lib/_util.py in _lazyselect(condlist, choicelist, arrays, default)
95             continue
96         cond, _ = np.broadcast_arrays(cond, arrays[0])
---> 97         temp = tuple(np.extract(cond, arr) for arr in arrays)
98         np.place(out, cond, func(*temp))
99     return out

~/miniconda3/envs/python3/lib/python3.6/site-packages/scipy/_lib/_util.py in <genexpr>(.0)
95             continue
96         cond, _ = np.broadcast_arrays(cond, arrays[0])
---> 97         temp = tuple(np.extract(cond, arr) for arr in arrays)
98         np.place(out, cond, func(*temp))
99     return out

~/miniconda3/envs/python3/lib/python3.6/site-packages/numpy/lib/function_base.py in extract(condition, arr)
2324
2325     """
-> 2326     return _nx.take(ravel(arr), nonzero(ravel(condition))[0])
2327
2328

~/miniconda3/envs/python3/lib/python3.6/site-packages/numpy/core/fromnumeric.py in take(a, indices, axis, out, mode)
132            [5, 7]])
133     """
--> 134     return _wrapfunc(a, 'take', indices, axis=axis, out=out, mode=mode)
135
136

~/miniconda3/envs/python3/lib/python3.6/site-packages/numpy/core/fromnumeric.py in _wrapfunc(obj, method, *args, **kwds)
55 def _wrapfunc(obj, method, *args, **kwds):
56     try:
---> 57         return getattr(obj, method)(*args, **kwds)
58
59     # An AttributeError occurs if the object does not have

KeyboardInterrupt:

``````

# Try VonMises

``````

In [112]:

kappa = 3.99390425811
mean, var, skew, kurt = stats.vonmises.stats(kappa, moments='mvsk')
print(mean, var, skew, kurt)

``````
``````

4.995631231732412e-17 0.2988042418167408 4.572598569626176e-16 0.6905715325848036

``````
``````

In [122]:

x = np.linspace(stats.vonmises.ppf(0.01, kappa), stats.vonmises.ppf(0.99, kappa), 100)

plt.plot(x, stats.vonmises.pdf(x, kappa), 'r-', lw=5, alpha=0.6, label='vonmises pdf')

x = np.linspace(stats.norm.ppf(0.01, loc=0, scale=.4, ), stats.norm.ppf(0.99, loc=0, scale=.4, ), 100)

plt.plot(x, stats.norm.pdf(x, loc=0, scale=.4), 'k', lw=2, alpha=0.6, label='normal pdf')

``````
``````

Out[122]:

[<matplotlib.lines.Line2D at 0x1228fb240>]

``````
``````

In [ ]:

``````
``````

In [ ]:

``````
``````

In [ ]:

``````
``````

In [ ]:

``````
``````

In [ ]:

``````
``````

In [ ]:

``````
``````

In [ ]:

``````
``````

In [ ]:

``````