pythonscipystatistics

Scipy's WrappedCauchy function wrong?

I'd like someone to check my understanding on the wrapped cauchy function in Scipy... From Wikipedia "a wrapped Cauchy distribution is a wrapped probability distribution that results from the "wrapping" of the Cauchy distribution around the unit circle." It's similar to the Von Mises distribution in that way.

I use the following bits of code to calculate a couple thousand random variates, get a histogram and plot it.

from scipy.stats import wrapcauchy, vonmises
import plotly.graph_objects as go
import numpy as np 

def plot_cauchy(c, loc = 0, scale = 1, size = 100000):
    '''
    rvs(c, loc=0, scale=1, size=1, random_state=None)
    '''
    rvses = vonmises.rvs(c, 
                        loc = loc, 
                        scale = scale,
                        size = size) 
    
    # rvses =  wrapcauchy.rvs(c, 
    #                         loc = loc, 
    #                         scale = scale,
    #                         size = size)
    
    y,x = np.histogram(rvses, 
                   bins = 200, 
                   range = [-np.pi,np.pi],
                   density = True)
    return x,y
 
fig = go.Figure()

loc = -3
x,y = plot_cauchy(0.5, loc = loc)
fig.add_trace(go.Scatter(x=x, y=y,
                    mode='lines+markers',
                    name= f'Centered on {loc}'))

loc = 1.5
x,y = plot_cauchy(0.5, loc = loc)
fig.add_trace(go.Scatter(x=x, y=y,
                    mode='lines+markers',
                    name= f'Centered on {loc}'))

loc = 0
x,y = plot_cauchy(0.5, loc = loc)
fig.add_trace(go.Scatter(x=x, y=y,
                    mode='lines+markers',
                    name=f'Centered on {loc}')) 

fig.show()

When plotting this using the Von Mises distribution I get a couple of distributions that are wrapped from -pi to pi and centered on "loc":

enter image description here

When I replace the vonmises distribution with the wrapcauchy distribution I get a "non-wrapped" result, that to my eye just looks wrong. To plot this completely I have to adjust the ranges for the histogram

This is with Scipy version '1.15.2'.

Is there a way to correctly "wrap" the outputs of a the Scipy call, or another library that correctly wraps the output from -pi to pi?

Solution

  • Is there a way to correctly "wrap" the outputs of a the Scipy call

    You can use the modulo operator. The operation number % x wraps all output to the range [0, x). If you want the range to begin at a value other than 0, you can add and subtract a constant before and after the modulo operation to center it somewhere else. If you want the range to begin at -pi, you can do (array + pi) % (2 * pi) - pi.

    For example, this is how SciPy internally wraps the vonmises result.

    return np.mod(rvs + np.pi, 2*np.pi) - np.pi

    Source.

    You could do something similar with the result of scipy.stats.wrapcauchy().

    Here is how you could modify your code to do this:

    
from scipy.stats import wrapcauchy, vonmises
import plotly.graph_objects as go
import numpy as np 

def plot_cauchy_or_vm(c, loc = 0, scale = 1, kind="vonmises", size = 100000):
    '''
    rvs(c, loc=0, scale=1, size=1, random_state=None)
    '''
    if kind == "vonmises":
        rvses = vonmises.rvs(c, 
                            loc = loc, 
                            scale = scale,
                            size = size) 
    elif kind == "cauchy":
        rvses =  wrapcauchy.rvs(c, 
                                loc = loc, 
                                scale = scale,
                                size = size)
        rvses = ((rvses + np.pi) % (2 * np.pi)) - np.pi
    else:
        raise Exception("Unknown kind")
        
    
    y,x = np.histogram(rvses, 
                   bins = 200, 
                   range = [-np.pi,np.pi],
                   density = True)
    return x,y
 
    
for kind in ["vonmises", "cauchy"]:
    fig = go.Figure()

    loc = -3
    x,y = plot_cauchy_or_vm(0.5, kind=kind, loc = loc)
    fig.add_trace(go.Scatter(x=x, y=y,
                        mode='lines+markers',
                        name= f'Centered on {loc}'))

    loc = 1.5
    x,y = plot_cauchy_or_vm(0.5, kind=kind, loc = loc)
    fig.add_trace(go.Scatter(x=x, y=y,
                        mode='lines+markers',
                        name= f'Centered on {loc}'))

    loc = 0
    x,y = plot_cauchy_or_vm(0.5, kind=kind, loc = loc)
    fig.add_trace(go.Scatter(x=x, y=y,
                        mode='lines+markers',
                        name=f'Centered on {loc}')) 

    fig.show()

    Output:

    Cauchy Plot

    Cauchy Plot