Skip to content

Commit

Permalink
fix: _circ_mean_resample
Browse files Browse the repository at this point in the history
  • Loading branch information
@huangziwei
huangziwei committed Nov 20, 2024
1 parent b41aab8 commit 41e562b
Show file tree
Hide file tree
Showing 4 changed files with 105 additions and 122 deletions.
1 change: 1 addition & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
.DS_Store
.coverage
.vscode
.circstat

*.egg-info
*/__pycache__
Expand Down
74 changes: 22 additions & 52 deletions examples/B1-Fisher-1993.ipynb
View file

Large diffs are not rendered by default.

74 changes: 35 additions & 39 deletions examples/T1-descriptive-statistics.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
"cells": [
{
"cell_type": "code",
"execution_count": 11,
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
Expand Down Expand Up @@ -48,7 +48,7 @@
},
{
"cell_type": "code",
"execution_count": 12,
"execution_count": 2,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -81,7 +81,7 @@
},
{
"cell_type": "code",
"execution_count": 13,
"execution_count": 3,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -123,16 +123,16 @@
},
{
"cell_type": "code",
"execution_count": 14,
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(0.23986529263377035, 0.23986529263377035)"
"(np.float64(0.23986529263377035), np.float64(0.23986529263377035))"
]
},
"execution_count": 14,
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
Expand All @@ -159,16 +159,16 @@
},
{
"cell_type": "code",
"execution_count": 15,
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(-0.9235490965405332, -0.9235490965405332)"
"(np.float64(-0.9235490965405332), np.float64(-0.9235490965405332))"
]
},
"execution_count": 15,
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
Expand Down Expand Up @@ -196,16 +196,16 @@
},
{
"cell_type": "code",
"execution_count": 16,
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(6.642665308562964, 6.642665308562964)"
"(np.float64(6.642665308562964), np.float64(6.642665308562964))"
]
},
"execution_count": 16,
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
Expand Down Expand Up @@ -243,7 +243,7 @@
},
{
"cell_type": "code",
"execution_count": 17,
"execution_count": 7,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -289,7 +289,7 @@
},
{
"cell_type": "code",
"execution_count": 18,
"execution_count": 8,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -332,7 +332,7 @@
},
{
"cell_type": "code",
"execution_count": 19,
"execution_count": 9,
"metadata": {},
"outputs": [
{
Expand Down Expand Up @@ -361,20 +361,20 @@
},
{
"cell_type": "code",
"execution_count": 20,
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Table 4.1 of Fisher (1993)\n",
"Table 4.1 of Fisher (1993), p.78\n",
"Sample size 10 20 30 \n",
"----------------------------------------------------------------------------------------------------\n",
"median 279.0 250.0 245.0 \n",
"mean 280.8 248.7 247.6 \n",
"95% median CI [262.9 298.9] [222.5 271.1] [229.9 263.1] \n",
"95% bootstrap mean CI [245. 315.] [229. 277.] [229. 267.] \n",
"95% median CI [245. 315.] [229. 277.] [229. 267.] \n",
"95% bootstrap mean CI [262.3 299.9] [221.8 269.4] [230.3 261.2] \n",
"95% large-sample mean CI - - [232.7 262.5] \n"
]
}
Expand All @@ -396,16 +396,16 @@
" median = circ_median(e, method='count')\n",
" table['median'].append(np.rad2deg(median).round(1))\n",
"\n",
" # CI for median\n",
" median_lb, median_ub, ci = circ_median_ci(alpha=e)\n",
" table['95% median CI'].append(np.rad2deg([median_lb, median_ub]).round(1))\n",
"\n",
" # CI for mean using bootstrap\n",
" mean_lb, mean_ub = circ_mean_ci(\n",
" alpha=e,\n",
" method=\"bootstrap\",\n",
" )\n",
" table['95% median CI'].append(np.rad2deg([mean_lb, mean_ub]).round(1))\n",
"\n",
" # CI for median\n",
" median_lb, median_ub, ci = circ_median_ci(alpha=e)\n",
" table['95% bootstrap mean CI'].append(np.rad2deg([median_lb, median_ub]).round(1))\n",
" table['95% bootstrap mean CI'].append(np.rad2deg([mean_lb, mean_ub]).round(1))\n",
"\n",
" if i == 2:\n",
" # CI for mean using dispersion\n",
Expand All @@ -417,7 +417,7 @@
" else:\n",
" table['95% large-sample mean CI'].append('-')\n",
"\n",
"print('Table 4.1 of Fisher (1993)')\n",
"print('Table 4.1 of Fisher (1993), p.78')\n",
"print(f\"{'Sample size':25} {'10':25} {'20':25} {'30':25}\")\n",
"print(\"-\"*100)\n",
"for k, v in table.items():\n",
Expand All @@ -427,24 +427,25 @@
},
{
"cell_type": "code",
"execution_count": 21,
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Last updated: 2023-04-18 11:17:25CEST\n",
"Last updated: 2024-11-20 14:08:59CET\n",
"\n",
"Python implementation: CPython\n",
"Python version : 3.10.6\n",
"IPython version : 8.5.0\n",
"Python version : 3.10.13\n",
"IPython version : 8.29.0\n",
"\n",
"pycircstat2: 0.1.0\n",
"\n",
"numpy: 1.23.5\n",
"numpy : 2.1.3\n",
"pycircstat2: 0.1.0\n",
"\n",
"Watermark: 2.3.1\n",
"Watermark: 2.5.0\n",
"\n"
]
}
Expand All @@ -457,7 +458,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "base",
"display_name": ".circstat",
"language": "python",
"name": "python3"
},
Expand All @@ -471,14 +472,9 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.6"
"version": "3.10.13"
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "4a52f57b56c27f38fc3d95daa57af6da3929fe8541a384ec0d11efb6a1b206eb"
}
}
"orig_nbformat": 4
},
"nbformat": 4,
"nbformat_minor": 2
Expand Down
78 changes: 47 additions & 31 deletions pycircstat2/descriptive.py
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
from typing import Tuple, Union

import numpy as np
from scipy.stats import chi2, norm
from scipy.stats import chi2, norm, t

from .utils import angrange

Expand Down Expand Up @@ -347,8 +347,8 @@ def circ_std(


def circ_median(
alpha: np.array,
w: Union[np.array, None] = None,
alpha: np.ndarray,
w: Union[np.ndarray, None] = None,
grouped: bool = False,
method: str = "deviation",
return_average: bool = True,
Expand Down Expand Up @@ -666,10 +666,7 @@ def _circ_mean_ci_approximate(
/ R
)
else: # eq(26.25)
inner = (
np.sqrt(n**2 - (n**2 - R**2) * np.exp(chi2.isf(1 - ci, 1) / n))
/ R
)
inner = np.sqrt(n**2 - (n**2 - R**2) * np.exp(chi2.isf(1 - ci, 1) / n)) / R

d = np.arccos(inner)
lb = mean - d
Expand All @@ -684,9 +681,31 @@ def _circ_mean_ci_approximate(


def _circ_mean_ci_bootstrap(alpha, B=2000, ci=0.95, return_samples=False):
from scipy.stats import t

beta = np.array([_circ_mean_resample(alpha) for i in range(B)]).flatten()
# Precompute z0 and v0 from original data
# algo 1
X = np.cos(alpha)
Y = np.sin(alpha)
z1 = np.mean(X) # eq(8.24)
z2 = np.mean(Y)
z0 = np.array([z1, z2])

# algo 2
u11 = np.mean((X - z1) ** 2) # eq(8.25)
u22 = np.mean((Y - z2) ** 2)
u12 = u21 = np.mean((X - z1) * (Y - z2)) # eq(8.26)

β = (u11 - u22) / (2 * u12) - np.sqrt(
(u11 - u22) ** 2 / (4 * u12**2 + 1)
) # eq(8.27)
t1 = np.sqrt(β**2 * u11 + 2 * β * u12 + u22) / np.sqrt(1 + β**2) # eq(8.28)
t2 = np.sqrt(u11 - 2 * β * u12 + β**2 * u22) / np.sqrt(1 + β**2) # eq(8.29)
v11 = (β**2 * t1 + t2) / (1 + β**2) # eq(8.30)
v22 = (t1 + β**2 * t2) / (1 + β**2)
v12 = v21 = β * (t1 - t2) / (1 + β**2) # eq(8.31)
v0 = np.array([[v11, v12], [v21, v22]])

beta = np.array([_circ_mean_resample(alpha, z0, v0) for i in range(B)]).flatten()

lb, ub = t.interval(
ci,
Expand All @@ -701,40 +720,37 @@ def _circ_mean_ci_bootstrap(alpha, B=2000, ci=0.95, return_samples=False):
return lb, ub


def _circ_mean_resample(alpha):
def _circ_mean_resample(alpha, z0, v0):
"""
Implementation of Section 8.3.5 (Fisher, 1993)
Implementation of Section 8.3.5 (Fisher, 1993, P210)
"""

θ = np.random.choice(alpha, len(alpha), replace=True)
X = np.cos(θ)
Y = np.sin(θ)
z1 = np.mean(X)

# algo 1
z1 = np.mean(X) # eq(8.24)
z2 = np.mean(Y)
zB = np.array([z1, z2])

u11 = np.mean((X - z1) ** 2)
u11 = np.mean((X - z1) ** 2) # eq(8.25)
u22 = np.mean((X - z2) ** 2)
u12 = u21 = np.mean((X - z1) * (Y - z2))
u12 = u21 = np.mean((X - z1) * (Y - z2)) # eq(8.26)

β = (u11 - u22) / (2 * u12) - np.sqrt((u11 - u22) ** 2 / (4 * u12**2 + 1))
t1 = np.sqrt(β**2 * u11 + 2 * β * u12 + u22) / np.sqrt(1 + β**2)
t2 = np.sqrt(u11 - 2 * β * u12 + β**2 * u22) / np.sqrt(1 + β**2)
v11 = (β**2 * t1 + t2) / (1 + β**2)
v22 = (t1 + β**2 * t2) / (1 + β**2)
v12 = v21 = β * (t1 - t2) / (1 + β**2)

q1 = np.sqrt(1 + β**2) / np.sqrt(β**2 * u11 + 2 * β * u12 + u22)
q2 = np.sqrt(1 + β**2) / np.sqrt(u11 - 2 * β * u12 + β**2 * u22)
w11 = (β**2 * q1 + q2) / (1 + β**2)
w22 = (q1 + β**2 * q2) / (1 + β**2)
w12 = w21 = β * (t1 - q2) / (1 + β**2)
# algo 3
β = (u11 - u22) / (2 * u12) - np.sqrt(
(u11 - u22) ** 2 / (4 * u12**2 + 1)
) # eq(8.27)
t1 = np.sqrt(1 + β**2) / np.sqrt(β**2 * u11 + 2 * β * u12 + u22) # eq(8.33)
t2 = np.sqrt(1 + β**2) / np.sqrt(u11 - 2 * β * u12 + β**2 * u22) # eq(8.34)
w11 = (β**2 * t1 + t2) / (1 + β**2) # eq(8.35)
w22 = (t1 + β**2 * t2) / (1 + β**2)
w12 = w21 = β * (t1 - t2) / (1 + β**2) # eq(8.36)

z0 = np.array([z1, z2])
u0 = np.array([[u11, u12], [u21, u22]])
v0 = np.array([[v11, v12], [v21, v22]])
w0 = np.array([[w11, w12], [w21, w22]])
wB = np.array([[w11, w12], [w21, w22]])

Cbar, Sbar = z0 + v0 @ w0 @ (z0 - z0)
Cbar, Sbar = z0 + v0 @ wB @ (zB - z0)
Cbar = np.power(Cbar**2 + Sbar**2, -0.5) * Cbar
Sbar = np.power(Cbar**2 + Sbar**2, -0.5) * Sbar

Expand Down

0 comments on commit 41e562b

Please sign in to comment.