# Rayleigh Test circular stats

**URL:** <https://discourse.julialang.org/t/rayleigh-test-circular-stats/60882>\
**Category:** General Usage\
**Tags:** statistics\
**Created:** [May 10, 2021, 12:03pm UTC](https://discourse.julialang.org/t/rayleigh-test-circular-stats/60882 "2021-05-10T12:03:26Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![katdeane](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/katdeane/32/15883_2.png) [@katdeane](https://discourse.julialang.org/u/katdeane)\
**Post date:** [May 10, 2021, 12:03pm UTC](https://discourse.julialang.org/t/rayleigh-test-circular-stats/60882/1 "2021-05-10T12:03:26Z")

</div>

Dear Julia peeps,

I am trying to relate brain activity to an amplitude modulated tone stimulus. As in, does the peak of brain activity occur in synchrony with the phase of the modulated tone? HypothesisTests.jl has this function for a Rayleigh Test: ( [HypothesisTests.jl/circular.jl at master · JuliaStats/HypothesisTests.jl (github.com)](https://github.com/JuliaStats/HypothesisTests.jl/blob/master/src/circular.jl)) but I am not sure from the documentation how to use it. Below is the code I’m using to find instantaneous phase of a 2 Hz amplitude modulation wave and relate a series of mock peak latencies to it. I then take the phase as orientation for a series of unit vectors and find the sum of vectors for the length of the resultant vector (vstrength). I want to compare that resultant vector length between conditions and I believe I should be able to do that with the Rayleigh Test. Any help would be appreciated on how to use this function properly.

(I’m using this publication as reference: [Auditory Cortex Phase Locking to Amplitude-Modulated Cochlear Implant Pulse Trains (physiology.org)](https://journals.physiology.org/doi/pdf/10.1152/jn.01109.2007))

```julia
using DSP, FFTW, HypothesisTests

# prep the data which will fluctuate
peaklat1 = Int.(ceil.(rand(100)*100)) # list of fake latencies
peaklat1[peaklat1 .> 50] = peaklat1[peaklat1 .> 50].+150 # create more variance
peaklat2 = Int.(ceil.(rand(100)*100)) # list of fake latencies
clickfreq = 2 # stimulus frequency
sr = 1000 # sampling rate

function getAMwave(clickfreq=2,sr=1000)
    t = 0:1/sr:1 # time course
    wave = sin.((2*π*clickfreq*t).+3π/2) # AM wave starting at the bottom
    return wave
end

# generate the AM wave used and get instantaneous phase
wave = getAMwave(clickfreq,sr) # function above
wavephase = angle.(hilbert(wave)) # Compute phase of wave

# get the phase at each time point of signal peak amplitude
orientation1 = wavephase[peaklat1]

# break vector into components i and j to compute sums, averages, and then the resultant vector 
i1 = cos.(orientation1) 
j1 = sin.(orientation1)
sumi1 = sum(i1) 
sumj1 = sum(j1)
avgi1 = sumi1/length(i1) # average of component i
avgj1 = sumj1/length(j1) # average of component j

vstrength1 = sqrt((avgi1^2)+(avgj1^2)) # find the vector based on i and j components (Pythagorean)

# now the same for the second latency series
orientation2 = wavephase[peaklat2]
i2 = cos.(orientation2) 
j2 = sin.(orientation2)
sumi2 = sum(i2) 
sumj2 = sum(j2)
avgi2 = sumi2/length(i2) 
avgj2 = sumj2/length(j2) 
vstrength2 = sqrt((avgi2^2)+(avgj2^2)) 

# what I feel I should be able to do: 
RayleighTest(vstrength1,vstrength2)
# or 
RayleighTest([vstrength1,vstrength2],100) # 100 observations per condition

```

---

<div class="post-metadata">

**Author:** ![halleysfifthinc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/halleysfifthinc/32/206280_2.png) [@halleysfifthinc](https://discourse.julialang.org/u/halleysfifthinc)\
**Post date:** [May 10, 2021, 2:00pm UTC](https://discourse.julialang.org/t/rayleigh-test-circular-stats/60882/2 "2021-05-10T14:00:29Z")

</div>

It looks like the `RayleighTest` constructor is expecting the `i` and `j` components as a single complex vector. Since you are already calculating the resultant length (`vstrength`), you can use the default constructor as `RayleighTest(vstrength1, 100)`.

However, as implemented in HypothesisTests, it does not appear that you can make the kind of comparison that you wish to. I’m not familiar with the Rayleigh Test, so I don’t know if that is a limitation of the test itself or the implementation.

---

<div class="post-metadata">

**Author:** ![katdeane](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/katdeane/32/15883_2.png) [@katdeane](https://discourse.julialang.org/u/katdeane)\
**Post date:** [May 19, 2021, 12:08pm UTC](https://discourse.julialang.org/t/rayleigh-test-circular-stats/60882/3 "2021-05-19T12:08:26Z")

</div>

@ararslan, Sorry for the direct tag but I’m working on implementing the HypothesisTests.jl - Rayleigh test in my analysis. I have been able to figure out what I was misunderstanding in the original post here but now I’m comparing results with R and python versions of this test and I noticed that the p value is inconsistent between those two to Julia.

R and python only take `exp(-Z)` but it’s written here to only perform this calculation if the sample size is larger than `1e6`. Could you tell me why there is are is the extra equation in the p value calculations if sample size is less than a million?

```julia
function pvalue(x::RayleighTest)
    Z = x.Rbar^2 * x.n
    x.n > 1e6 ? exp(-Z) :
        exp(-Z)*(1+(2*Z-Z^2)/(4*x.n)-(24*Z - 132*Z^2 + 76*Z^3 - 9*Z^4)/(288*x.n^2))
end

```

---

<div class="post-metadata">

**Author:** ![katdeane](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/katdeane/32/15883_2.png) [@katdeane](https://discourse.julialang.org/u/katdeane)\
**Post date:** [May 19, 2021, 12:10pm UTC](https://discourse.julialang.org/t/rayleigh-test-circular-stats/60882/4 "2021-05-19T12:10:29Z")

</div>

Thank you! It took me a while but I understand now that it is indeed just the way the test works. 🙂

---

<div class="post-metadata">

**Author:** ![halleysfifthinc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/halleysfifthinc/32/206280_2.png) [@halleysfifthinc](https://discourse.julialang.org/u/halleysfifthinc)\
**Post date:** [May 19, 2021, 3:22pm UTC](https://discourse.julialang.org/t/rayleigh-test-circular-stats/60882/5 "2021-05-19T15:22:35Z")

</div>

That equation is found [here](http://wexler.free.fr/library/files/wilkie%20(1983)%20rayleigh%20test%20for%20randomness%20of%20circular%20data.pdf). My guess is that the latter method is more accurate, albeit more expensive, and that for `n > 1e6`, `exp(-Z)` is likely sufficiently accurate (the exponential is going to dominate) and computationally cheaper to make sense.

---

<div class="post-metadata">

**Author:** ![katdeane](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/katdeane/32/15883_2.png) [@katdeane](https://discourse.julialang.org/u/katdeane)\
**Post date:** [May 20, 2021, 6:12pm UTC](https://discourse.julialang.org/t/rayleigh-test-circular-stats/60882/6 "2021-05-20T18:12:30Z")

</div>

The longer computation put me a couple of orders of magnitude off of the `exp(-Z)` computation so I was unsure. But it’s great to see the source here, thank you for finding that!