The Bivariate Poisson Conditionals Distribution (BPCD)

Introduction

This vignette introduces the Bivariate Poisson Conditionals Distribution (BPCD), defined via conditional specifications, as proposed by Ghosh, Marques, and Chakraborty (2021). The BCD package provides functions to evaluate the joint and cumulative distributions, perform random sampling, and estimate parameters via maximum likelihood.

Joint Probability: dpoisBCD()

The joint probability mass function (p.m.f.) of the BBCD is given by:

\[ P(X = x, Y = y) = K \frac{\lambda_1^x \lambda_2^y \lambda_3^{xy}}{x! y!}, \]

where \(K\) is a normalizing constant ensuring the probabilities sum to 1 and .

Example

dpoisBCD(x = 1, y = 2, lambda1 = 0.5, lambda2 =  0.5, lambda3 =  0.5)
#> [1] 0.006358011
dpoisBCD(x = 0, y = 1, lambda1 = 0.5, lambda2 =  0.5, lambda3 =  1)
#> [1] 0.1839397

Cumulative Distribution: ppoisBCD()

The function ppoisBCD() computes the cumulative distribution:

\[ P(X \leq x, Y \leq y) \]

Example

ppoisBCD(x = 1, y = 1, lambda1 = 0.5, lambda2 = 0.5, lambda3 = 0.5)
#> [1] 0.8646896
ppoisBCD(x = 2, y = 2, lambda1 = 1.0, lambda2 = 1.0, lambda3 = 0.8)
#> [1] 0.8872704

Random Sampling: rpoisBCD()

Generate samples from the BPCD using:

rpoisBCD(n, lambda1, lambda2, lambda3)

Example

set.seed(123)
samples <- rpoisBCD(n = 100, lambda1 = 0.5, lambda2 = 0.5, lambda3 = 0.5)
head(samples)
#>   X Y
#> 1 0 1
#> 2 1 2
#> 3 0 1
#> 4 0 2
#> 5 1 0
#> 6 1 0
cor(samples$X, samples$Y) # Should be negative
#> [1] -0.2450123

Maximum Likelihood Estimation: MLEpoisBCD()

Estimate the parameters of the distribution from data.

Example

data <- rpoisBCD(n = 100, lambda1 = 0.5, lambda2 = 0.5, lambda3 = 0.5)
fit <- MLEpoisBCD(data)
fit
#> $lambda1
#> [1] 0.5252073
#> 
#> $lambda2
#> [1] 0.5058323
#> 
#> $lambda3
#> [1] 0.431483
#> 
#> $logLik
#> [1] -162.3016
#> 
#> $AIC
#> [1] 330.6032
#> 
#> $BIC
#> [1] 338.4187
#> 
#> $convergence
#> [1] 0
#> 
#> attr(,"class")
#> [1] "MLEpoisBCD"

Real Data Example

The dataset lensfaults records counts of surface faults \(X\) and interior faults \(Y\) observed in 100 optical lenses.

data(lensfaults)
head(lensfaults)
#>   X Y
#> 1 0 0
#> 2 0 1
#> 3 0 2
#> 4 0 2
#> 5 0 2
#> 6 0 2
plot(lensfaults$X, lensfaults$Y, xlab = "X", ylab = "Y")

fit <- MLEpoisBCD(lensfaults)
FTtest(lensfaults, "BPCD", params = fit, num_params = 3)
#> $observed
#>    0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
#> 0  1 1 4 0 0 0 0 0 0 0  0  0  0  0  1
#> 1  3 2 6 2 5 2 0 0 0 0  0  0  0  0  0
#> 2  1 2 4 3 2 1 1 1 0 0  1  0  0  0  0
#> 3  0 5 1 2 3 2 0 0 0 0  0  0  0  0  0
#> 4  0 1 2 5 3 1 1 0 0 0  0  0  0  0  0
#> 5  0 1 2 1 2 1 0 0 0 0  1  0  0  0  0
#> 6  0 0 2 0 1 1 0 0 0 0  0  0  1  0  0
#> 7  0 1 0 1 1 0 0 0 0 0  0  0  0  0  0
#> 8  0 0 2 0 0 0 0 0 0 0  0  0  0  0  0
#> 9  0 0 0 0 0 0 0 0 0 0  0  0  0  0  0
#> 10 0 0 0 0 0 0 0 0 0 0  0  0  0  0  0
#> 11 0 1 0 0 0 0 0 0 0 0  0  0  0  0  0
#> 12 0 1 0 0 0 0 0 0 0 0  0  0  0  0  0
#> 
#> $expected
#>              0            1           2           3           4            5
#> 0  0.186393251 0.5780419270 0.896310534 0.926544756 0.718349125 0.4455482262
#> 1  0.571730986 1.7730496060 2.749286800 2.842025357 2.203419118 1.3666467257
#> 2  0.876845910 2.7192706605 4.216495076 4.358725296 3.379314905 2.0959832888
#> 3  0.896527410 2.7803068416 4.311137643 4.456560333 3.455166301 2.1430292918
#> 4  0.687488008 2.1320347711 3.305928403 3.417443517 2.649540182 1.6433484596
#> 5  0.421751532 1.3079339869 2.028079549 2.096490444 1.625406724 1.0081408295
#> 6  0.215609039 0.6686457999 1.036800699 1.071773915 0.830945132 0.5153846739
#> 7  0.094477952 0.2929946074 0.454316791 0.469641741 0.364112723 0.2258369531
#> 8  0.036224469 0.1123391637 0.174192859 0.180068708 0.139607070 0.0865897658
#> 9  0.012345852 0.0382869017 0.059367585 0.061370164 0.047580220 0.0295111139
#> 10 0.003786889 0.0117438844 0.018210041 0.018824300 0.014594459 0.0090520543
#> 11 0.001055970 0.0032747687 0.005077849 0.005249134 0.004069648 0.0025241549
#> 12 0.000269918 0.0008370684 0.001297956 0.001341739 0.001040249 0.0006452029
#>               6            7            8            9           10
#> 0  0.2302886982 0.1020243276 3.954967e-02 1.362792e-02 4.226284e-03
#> 1  0.7063731307 0.3129430329 1.213122e-01 4.180143e-02 1.296344e-02
#> 2  1.0833423516 0.4799509302 1.860527e-01 6.410954e-02 1.988162e-02
#> 3  1.1076588276 0.4907238085 1.902288e-01 6.554853e-02 2.032788e-02
#> 4  0.8493908297 0.3763038694 1.458740e-01 5.026486e-02 1.558812e-02
#> 5  0.5210736473 0.2308501845 8.948895e-02 3.083586e-02 9.562803e-03
#> 6  0.2663847787 0.1180158997 4.574880e-02 1.576400e-02 4.888724e-03
#> 7  0.1167274268 0.0517135114 2.004671e-02 6.907642e-03 2.142195e-03
#> 8  0.0447552998 0.0198278483 7.686254e-03 2.648509e-03 8.213543e-04
#> 9  0.0152532893 0.0067576333 2.619593e-03 9.026522e-04 2.799301e-04
#> 10 0.0046786984 0.0020727941 8.035175e-04 2.768739e-04 8.586401e-05
#> 11 0.0013046497 0.0005779963 2.240599e-04 7.720596e-05 2.394308e-05
#> 12 0.0003334834 0.0001477425 5.727228e-05 1.973473e-05 6.120125e-06
#>              11           12           13           14
#> 0  1.191503e-03 3.079237e-04 7.345627e-05 1.627159e-05
#> 1  3.654742e-03 9.445058e-04 2.253152e-04 4.991046e-05
#> 2  5.605163e-03 1.448559e-03 3.455588e-04 7.654611e-05
#> 3  5.730976e-03 1.481073e-03 3.533152e-04 7.826425e-05
#> 4  4.394709e-03 1.135738e-03 2.709342e-04 6.001571e-05
#> 5  2.696011e-03 6.967381e-04 1.662093e-04 3.681769e-05
#> 6  1.378263e-03 3.561885e-04 8.497001e-05 1.882204e-05
#> 7  6.039423e-04 1.560786e-04 3.723310e-05 8.247652e-06
#> 8  2.315619e-04 5.984322e-05 1.427581e-05 3.162291e-06
#> 9  7.891982e-05 2.039548e-05 4.865414e-06 1.077757e-06
#> 10 2.420737e-05 6.255982e-06 1.492387e-06 3.305845e-07
#> 11 6.750197e-06 1.744474e-06 4.161503e-07 9.218310e-08
#> 12 1.725428e-06 4.459074e-07 1.063728e-07 2.356306e-08
#> 
#> $test
#> $test$statistic
#> [1] 112.0946
#> 
#> $test$df
#> [1] 191
#> 
#> $test$p_value
#> [1] 0.999999

The dataset eplSeasonGoals is a list of data frames for five consecutive seasons (2014/15 to 2018/19) from the English Premier League. Each data frame contains the number of full-time home (\(X\)) and away (\(Y\)) goals scored in each match of the season.

data(eplSeasonGoals)
plot(eplSeasonGoals[["1415"]]$X, eplSeasonGoals[["1415"]]$Y, xlab = "X", ylab = "Y")

plot(eplSeasonGoals[["2425"]]$X, eplSeasonGoals[["2425"]]$Y, xlab = "X", ylab = "Y")

fit <- MLEpoisBCD(eplSeasonGoals[["1415"]])
FTtest(eplSeasonGoals[["1415"]], "BPCD", params = fit, num_params = 3)
#> $observed
#>    0  1  2  3 4 5 6
#> 0 31 40 14  5 1 1 0
#> 1 37 37 26 15 4 0 0
#> 2 37 35 22  6 2 0 0
#> 3 19 14  9  3 0 0 1
#> 4  4  5  2  1 0 0 0
#> 5  2  1  0  2 0 0 0
#> 6  1  2  0  0 0 0 0
#> 7  0  0  0  0 0 0 0
#> 8  1  0  0  0 0 0 0
#> 
#> $expected
#>             0           1            2           3            4            5
#> 0 28.03424236 31.78629927 18.020262649 6.810688324 1.9305549127 0.4377874416
#> 1 42.47791655 46.93740825 25.932537095 9.551677511 2.6386121498 0.5831247187
#> 2 32.18159725 34.65518704 18.659452787 6.697894466 1.8031794982 0.3883556325
#> 3 16.25401755 17.05792223  8.950793541 3.131163102 0.8215066906 0.1724274899
#> 4  6.15708454  6.29716795  3.220219234 1.097828104 0.2807013570 0.0574175513
#> 5  1.86586190  1.85974932  0.926828386 0.307930697 0.0767304785 0.0152958218
#> 6  0.47119712  0.45770146  0.222296163 0.071976443 0.0174787373 0.0033956250
#> 7  0.10199502  0.09655247  0.045700170 0.014420522 0.0034127574 0.0006461299
#> 8  0.01931805  0.01782184  0.008220755 0.002528015 0.0005830541 0.0001075791
#>              6
#> 0 8.273003e-02
#> 1 1.073906e-01
#> 2 6.970100e-02
#> 3 3.015926e-02
#> 4 9.787315e-03
#> 5 2.540952e-03
#> 6 5.497283e-04
#> 7 1.019420e-04
#> 8 1.654116e-05
#> 
#> $test
#> $test$statistic
#> [1] 43.69272
#> 
#> $test$df
#> [1] 59
#> 
#> $test$p_value
#> [1] 0.9319867
fit <- MLEpoisBCD(eplSeasonGoals[["2425"]])
FTtest(eplSeasonGoals[["2425"]], "BPCD", params = fit, num_params = 3)
#> $observed
#>    0  1  2  3 4 5 6
#> 0 16 26 28 12 5 3 1
#> 1 29 45 29 11 4 2 0
#> 2 21 33 31  5 2 1 1
#> 3 11 17 11  1 1 0 1
#> 4  7 10  6  3 0 0 0
#> 5  2  2  1  1 0 0 0
#> 6  0  0  0  0 0 0 0
#> 7  1  0  0  0 0 0 0
#> 
#> $expected
#>            0          1           2           3           4            5
#> 0 16.3518431 26.9702903 22.24203583 12.22847193 5.042328193 1.6633361060
#> 1 28.4590851 42.3127044 31.45506869 15.58903503 5.794408264 1.7230145191
#> 2 24.7653894 33.1914290 22.24214893  9.93656504 3.329331793 0.8924170594
#> 3 14.3673982 17.3576072 10.48507621  4.22242619 1.275304246 0.3081452880
#> 4  6.2513290  6.8079292  3.70704367  1.34570268 0.366380032 0.0798002886
#> 5  2.1759884  2.1361423  1.04851295  0.34310430 0.084205367 0.0165326846
#> 6  0.6311892  0.5585526  0.24713750  0.07289906 0.016127480 0.0028543094
#> 7  0.1569335  0.1251847  0.04992944  0.01327611 0.002647564 0.0004223884
#>              6
#> 0 0.4572436409
#> 1 0.4269603409
#> 2 0.1993413537
#> 3 0.0620463175
#> 4 0.0144842456
#> 5 0.0027049905
#> 6 0.0004209731
#> 7 0.0000561560
#> 
#> $test
#> $test$statistic
#> [1] 39.71118
#> 
#> $test$df
#> [1] 52
#> 
#> $test$p_value
#> [1] 0.8941378

Reference: Ghosh, I., Marques, F., and Chakraborty, S. (2021). A new bivariate Poisson distribution via conditional specification: properties and applications. , 48(16), 3025-3047.