| Title: | An Empirical Bayes Method for Chi-Squared Data |
|---|---|
| Description: | We provide the main R functions to compute the posterior interval for the noncentrality parameter of the chi-squared distribution. The skewness estimate of the posterior distribution is also available to improve the coverage rate of posterior intervals. Details can be found in Du and Hu (2022) <doi:10.1080/01621459.2020.1777137>. |
| Authors: | Lilun Du [aut, cre], Inchi Hu [aut] |
| Maintainer: | Lilun Du <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.1.1 |
| Built: | 2026-05-15 05:44:22 UTC |
| Source: | https://github.com/dulilun/ebchs |
The l_1 to l_4 derivative from the g-modeling method
density_g_model(x, k, pi_0, lambda_set, g_prior)density_g_model(x, k, pi_0, lambda_set, g_prior)
x |
a sequence of chi-squared test statistics |
k |
degrees of freedom |
pi_0 |
the proportion of the null |
lambda_set |
the set of noncentrality values |
g_prior |
the prior probability for the noncentrality values |
a list: the margianl density, and its first-to-fourth derivatives
Assuming the log density of the chi-squared statistics admits a parametric form, this function estimates up to the fourth order log-density derivatives.
density_LS(x)density_LS(x)
x |
a sequence of chi-squared test statistics |
a list: the first-to-fourth log density derivatives
p = 1000 k = 7 # the prior distribution for lambda alpha = 2 beta = 10 # lambda lambda = rep(0, p) pi_0 = 0.8 p_0 = floor(p*pi_0) p_1 = p-p_0 lambda[(p_0+1):p] = stats::rgamma(p_1, shape = alpha, rate=1/beta) # Generate a Poisson RV J = sapply(1:p, function(x){rpois(1, lambda[x]/2)}) X = sapply(1:p, function(x){rchisq(1, k+2*J[x])}) out = density_LS(X)p = 1000 k = 7 # the prior distribution for lambda alpha = 2 beta = 10 # lambda lambda = rep(0, p) pi_0 = 0.8 p_0 = floor(p*pi_0) p_1 = p-p_0 lambda[(p_0+1):p] = stats::rgamma(p_1, shape = alpha, rate=1/beta) # Generate a Poisson RV J = sapply(1:p, function(x){rpois(1, lambda[x]/2)}) X = sapply(1:p, function(x){rchisq(1, k+2*J[x])}) out = density_LS(X)
The semiparametric model is employed to estimate the log density derivatives of the chi-squared statistics.
density_PLS(x, qq)density_PLS(x, qq)
x |
a sequence of chi-squared test statistics |
qq |
the quantiles used for splines |
a list: the first and second density derivatives
p = 1000 k = 7 # the prior distribution for lambda alpha = 2 beta = 10 # lambda lambda = rep(0, p) pi_0 = 0.5 p_0 = floor(p*pi_0) p_1 = p-p_0 lambda[(p_0+1):p] = stats::rgamma(p_1, shape = alpha, rate=1/beta) # Generate a Poisson RV J = sapply(1:p, function(x){rpois(1, lambda[x]/2)}) X = sapply(1:p, function(x){rchisq(1, k+2*J[x])}) qq = c(0.2, 0.4, 0.6, 0.8) out = density_PLS(X, qq)p = 1000 k = 7 # the prior distribution for lambda alpha = 2 beta = 10 # lambda lambda = rep(0, p) pi_0 = 0.5 p_0 = floor(p*pi_0) p_1 = p-p_0 lambda[(p_0+1):p] = stats::rgamma(p_1, shape = alpha, rate=1/beta) # Generate a Poisson RV J = sapply(1:p, function(x){rpois(1, lambda[x]/2)}) X = sapply(1:p, function(x){rchisq(1, k+2*J[x])}) qq = c(0.2, 0.4, 0.6, 0.8) out = density_PLS(X, qq)
Give a sequence of chi-squared statistic values, the function computes the posterior mean, variance, and skewness of the non-centrality parameter given the data.
EB_CS( x, df, qq = c(0.2, 0.4, 0.6, 0.8), method = c("LS", "PLS", "g_model"), mixture = FALSE )EB_CS( x, df, qq = c(0.2, 0.4, 0.6, 0.8), method = c("LS", "PLS", "g_model"), mixture = FALSE )
x |
a sequence of chi-squared test statistics |
df |
the degrees of freedom |
qq |
the quantiles used in spline basis |
method |
LS: parametric least-squares; PLS: penalized least-squares; g-model: g-modeling |
mixture |
default is FALSE: there is no point mass at zero. |
a list: posterior mean, variance, and skewness estimates
Du and Hu (2022), An Empirical Bayes Method for Chi-Squared Data, Journal of American Statistical Association, forthcoming.
p = 1000 k = 7 # the prior distribution for lambda alpha = 2 beta = 10 # lambda lambda = rep(0, p) pi_0 = 0.8 p_0 = floor(p*pi_0) p_1 = p-p_0 lambda[(p_0+1):p] = rgamma(p_1, shape = alpha, rate=1/beta) # Generate a Poisson RV J = sapply(1:p, function(x){rpois(1, lambda[x]/2)}) X = sapply(1:p, function(x){rchisq(1, k+2*J[x])}) qq_set = seq(0.01, 0.99, 0.01) out = EB_CS(X, k, qq=qq_set, method='LS', mixture = TRUE) E = out$E_lambda V = out$V_lambda S = out$S_lambdap = 1000 k = 7 # the prior distribution for lambda alpha = 2 beta = 10 # lambda lambda = rep(0, p) pi_0 = 0.8 p_0 = floor(p*pi_0) p_1 = p-p_0 lambda[(p_0+1):p] = rgamma(p_1, shape = alpha, rate=1/beta) # Generate a Poisson RV J = sapply(1:p, function(x){rpois(1, lambda[x]/2)}) X = sapply(1:p, function(x){rchisq(1, k+2*J[x])}) qq_set = seq(0.01, 0.99, 0.01) out = EB_CS(X, k, qq=qq_set, method='LS', mixture = TRUE) E = out$E_lambda V = out$V_lambda S = out$S_lambda
Predictive recursion by Newton (2002)
predictive_recursion(x, k)predictive_recursion(x, k)
x |
a sequence of chi-squared test statistics |
k |
degrees of freedom |
a list: null proportion, prior probability, and lambda-mesh values
set.seed(2021) p = 1000 k = 7 # the prior distribution for lambda alpha = 2 beta = 10 # lambda lambda = rep(0, p) pi_0 = 0 p_0 = floor(p*pi_0) p_1 = p-p_0 lambda[(p_0+1):p] = stats::rgamma(p_1, shape = alpha, rate=1/beta) # Generate a Poisson RV J = sapply(1:p, function(x){rpois(1, lambda[x]/2)}) X = sapply(1:p, function(x){rchisq(1, k+2*J[x])}) out = predictive_recursion(X, k)set.seed(2021) p = 1000 k = 7 # the prior distribution for lambda alpha = 2 beta = 10 # lambda lambda = rep(0, p) pi_0 = 0 p_0 = floor(p*pi_0) p_1 = p-p_0 lambda[(p_0+1):p] = stats::rgamma(p_1, shape = alpha, rate=1/beta) # Generate a Poisson RV J = sapply(1:p, function(x){rpois(1, lambda[x]/2)}) X = sapply(1:p, function(x){rchisq(1, k+2*J[x])}) out = predictive_recursion(X, k)