統計検定：ポワソン分布と条件付き分布

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
def select_most_frequent_point(lambda_mu):
    """ 
    Description
        return the tuple of the most frequent datapoint and lambda_mu
        
        Note:
          When the most frequent datapoint are more than one datapoint, it returns the smaller integer datapoint
      
    INPUT
        lambda_mu:
            poisson lambda
    
    Returns:
        x: integer
            the most frequent datapoint based on the poisson distribution with the input lambda.

        x_theory: integer
           the analytical solution of the most frequent datapoint

        x_hacked: integer
           the datapoint which is caucluated by an easy way
           max(np.int64(np.ceil(lambda_mu - eps)) - 1, 0)

        lambda: float
            the poisson parameter.
            
    """
    eps = 1e-18

    x_range = np.arange(0, 100)
    x_index = np.argmax(np.around(poisson.pmf(x_range, lambda_mu), decimals=8)) ##HACK: 有効桁は小数点以下8桁
    x_theory = np.min(x_range[x_range >= lambda_mu - 1])
    x_tips = max(np.int64(np.ceil(lambda_mu - eps)) - 1, 0)

    return x_range[x_index], x_theory, x_tips, lambda_mu

### simulation lambda_range
lambda_range = np.linspace(0, 50, 201)

### create the object which can calculate the most frequent datapoint and compute vectorization
vectorized_select_most_frequent_point = np.vectorize(select_most_frequent_point)

### caculate the simulation
x_simulated, x_theory, x_hacked, lambda_param = vectorized_select_most_frequent_point(lambda_range)

### compare x_simulated with x_theory
print(np.array_equal(x_simulated, x_theory, equal_nan=True))

### compare x_theory with x_hacked
print(np.array_equal(x_theory, x_hacked, equal_nan=True))

1
2
3
4
5
6
7
8
9
10
11
from scipy.special import lambertw

## Generating Data
y = np.array([1, 2, 3, 4, 5])
freq = np.array([15, 12, 10, 3, 2])

## Estimation
lambert_w = np.average(y, weights = freq)
estimated_lambda = lambertw(-lambert_w/np.exp(lambert_w)) + lambert_w
print(estimated_lambda)
>>> 1.813224064133259 + 0j

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
def binary_search(func, y_min, y_max, value, eps):
    left = y_min             # left seach area
    right = y_max            # right seach area
    while left <= right:
        mid = (left + right) / 2            # calculate the mid
        if abs(func(mid) - value) < eps:
            # return the numerical solution
            return mid
        elif func(mid) < value:
            # set the left at the mid beacuse the objective function is a strictly increasing function
            left = mid
        else:
            # set the right at the mid
            right = mid 
    return None            # cannot compute

def objective_function(x):
    return x/(1 - np.exp(-x))

y = np.array([1, 2, 3, 4, 5])
freq = np.array([15, 12, 10, 3, 2])
bar_y = np.average(y, weights = freq)

print(binary_search(func=objective_function, y_min=0, y_max =6, value = bar_y, eps = 1e-8))
>>> 1.8132240548729897

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
from scipy.optimize import minimize
import math
import numpy as np

class ConditionalPoissonRegression:

    def __init__(self, data):
        self.data = data

    @staticmethod
    def conditiona_poisson_pmf(x, theta):
        vec_factorial = np.vectorize(math.factorial)
        return (theta**x * np.exp(-theta) / vec_factorial(x)) / (1 - np.exp(-theta))
    
    def neg_loglike(self, theta, _data):
        return np.sum(np.log(self.conditiona_poisson_pmf(x = _data, theta = theta)))

    def fit(self):
        observed_mean = np.mean(self.data)
        f = lambda theta: -self.neg_loglike(theta = theta, _data = self.data)
        return minimize(f, observed_mean, method = 'Nelder-Mead', options={'disp': True})
## Generating Data
y = np.array([1, 2, 3, 4, 5])
freq = np.array([15, 12, 10, 3, 2])
data = np.repeat(y, freq, axis=0)

poisson_mle = ConditionalPoissonRegression(data=data)

poisson_mle.fit()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
y = np.array([1, 2, 3, 4, 5])
freq = np.array([15, 12, 10, 3, 2])

## MoM var
def calculate_conditionalvar(x):
    denominator = (1 - np.exp(-x))**2
    numerator = x - x*(x+1)*np.exp(-x)

    return numerator/denominator

theory_var = calculate_conditionalvar(x = lambert_w)

## 不偏分散
data_var = np.average((y - lambert_w)**2, weights = freq)*np.sum(freq)/(np.sum(freq) - 1)

print(data_var, theory_var)
>>> 1.3130081300813008 1.7610396317212755

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
## Library
import numpy as np
import pandas as pd
from scipy import stats
from scipy.stats import poisson
from scipy.stats import binom
from statsmodels.distributions.empirical_distribution import ECDF
from itertools import chain
import matplotlib.pyplot as plt

## Generating Data
np.random.seed(seed=42)
mu_x, mu_y = 3, 2
N = 1000

x, y = poisson.rvs(mu_x, size=N), poisson.rvs(mu_y, size=N)
z = poisson.rvs(mu_x + mu_y, size=N)

1
2
3
4
5
6
7
8
9
fig, ax = plt.subplots(1, 1,figsize=(10, 7))
bins_range = np.arange(0, max(chain(x, y)), 1)

ax.hist(x, density = True, bins = bins_range, alpha = 0.5, label = '$\lambda = 3$')
ax.hist(y, density = True, bins = bins_range, alpha = 0.5, label = '$\lambda = 2$')
ax.set_title('random sampling from poisson', fontsize=15)
ax.set_xlabel('score', fontsize=12)
ax.set_ylabel('density', fontsize=12)
plt.legend();

1
2
3
4
5
6
7
8
9
fig, ax = plt.subplots(1, 1,figsize=(10, 7))
bins_range = np.arange(0, max(x + y), 1)

ax.hist(x + y, density = True, bins = bins_range, alpha = 0.5, label = '$x + y$')
ax.hist(z, density = True, bins = bins_range, alpha = 0.5, label = '$\lambda = 5$')
ax.set_title('random sampling from poisson', fontsize=15)
ax.set_xlabel('score', fontsize=12)
ax.set_ylabel('density', fontsize=12)
plt.legend();

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
def qqplot(x, y, title = None, xlabel = None, ylabel = None, plot = False, ax = None):
    """ 
    INPUT
        x, y : array_like, 1-Dimensional
            Two arrays of sample observations
    
    Returns:
        statistic: float
            KS statistic.
        pvalue: float
            One-tailed or two-tailed p-value.
        Figure:
            QQ-plot
            x-axis and y-axis is normalized between 0 and 1
    """
    
    sorted_x = sorted(x)
    x_xdf, y_cdf = ECDF(x), ECDF(y)
    ks_test = stats.ks_2samp(x, y)

    print(ks_test)

    if ax is None:
        fig, ax = plt.subplots(1, 1,figsize=(9, 9))
    
    if plot:
        ax.plot(x_xdf(sorted_x), y_cdf(sorted_x),lw=2, color='#0485d1')
    else:
        ax.scatter(x_xdf(sorted_x), y_cdf(sorted_x),lw=2, color='#0485d1')
    
    ax.axline([0, 0], [1, 1],lw=4, color='#A60628')
    ax.set_title(title  + '  (ks-test p-value: %.3f )'%ks_test[1] , fontsize=15)
    ax.set_xlabel(xlabel, fontsize=13)
    ax.set_ylabel(ylabel, fontsize=13)
    plt.legend()

1
2
3
4
5
6
7
8
9
10
11
12
13
## Generating Data
w = np.vstack([x, y, x+y])
x_2, y_2, z_2 = w[:, w[2,:] == 4]

## plot histogram
fig, ax = plt.subplots(1, 1,figsize=(10, 7))
bins_range = np.arange(0, max(chain(x, y)), 1)

ax.hist(x_2, density = True, bins = bins_range, alpha = 0.5, label = '$X_2\sim (X|x+y =4)$')
ax.set_title('random sampling from trunctaed poisson', fontsize=15)
ax.set_xlabel('$X_2$ sampling', fontsize=12)
ax.set_ylabel('density', fontsize=12)
plt.legend();

1
2
3
4
5
6
r_x = binom.rvs(n=4, p=mu_x/(mu_x + mu_y), size=1000)
qqplot(x_2, r_x, 
        title = 'QQ-plot between two poisson',
        xlabel = '$X_2$ from truncated X',
        ylabel = 'binomial distribution $n = 4, p = 0.6$',
        ax = None)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
## setting model parameter
N = 100
observed_lambda = 1.79

## geenrating data
x_range = np.arange(0, 8)
x = np.array([22, 23, 26, 18, 6, 4, 1, 0])

## poisson object
poisson_rv = poisson(observed_lambda)

## plot
fig, ax = plt.subplots(1, 1,figsize=(10, 7))
ax.bar(x_range, x, alpha = 0.5)
ax.plot(x_range, poisson_rv.pmf(x_range)*N)
ax.set_xlabel('outcome')
ax.set_ylabel('Frequency')
ax.set_title('compare the observed data and a theoretical distirbution')
;

1
2
3
4
5
6
7
from scipy.stats import chisquare

chisquare_poisson_pmf = lambda x: np.where(x < 7, 1, 0) * poisson_rv.pmf(x) + (1 - np.where(x < 7, 1, 0)) * (1 - poisson_rv.cdf(x-1))
theoretical_x = chisquare_poisson_pmf(x_range)*N

chisquare(x, f_exp=theoretical_x)
>>> Power_divergenceResult(statistic=4.873492084753685, pvalue=0.6753982780281518)