Skip to content

Instantly share code, notes, and snippets.

@songcai
Last active December 22, 2015 18:59
Show Gist options
  • Save songcai/6516959 to your computer and use it in GitHub Desktop.
Save songcai/6516959 to your computer and use it in GitHub Desktop.
Test run for Song's drmdel package

Testing drmdel package

This is a test for my drmdel package. To run it you have to install it from CRAN in advance.

# Load the package first.
library(drmdel)

#Try the main function drmdel for density ratio modelling of multiple samples.
example(drmdel, package = "drmdel")

## drmdel> # Data generation
## drmdel> set.seed(25)
## 
## drmdel> n_samples <- c(100, 200, 180, 150, 175)  # sample sizes
## 
## drmdel> x0 <- rgamma(n_samples[1], shape=5, rate=1.8)
## 
## drmdel> x1 <- rgamma(n_samples[2], shape=12, rate=1.2)
## 
## drmdel> x2 <- rgamma(n_samples[3], shape=12, rate=1.2)
## 
## drmdel> x3 <- rgamma(n_samples[4], shape=18, rate=5)
## 
## drmdel> x4 <- rgamma(n_samples[5], shape=25, rate=2.6)
## 
## drmdel> x <- c(x0, x1, x2, x3, x4)
## 
## drmdel> # Fit a DRM with the basis function q(x) = (x, log(abs(x))), which
## drmdel> # is the basis function for gamma family.
## drmdel> 
## drmdel> # There are 11 built-in basis function in drmdel(). And q(x) = (x,
## drmdel> # log(abs(x))) is the 6th basis function, so we can fit the model
## drmdel> # by specifying basis_func=6 in drmdel() as follows:
## drmdel> drmfit <- drmdel(x=x, n_samples=n_samples, basis_func=6)
## 
## drmdel> names(drmfit)
## [1] "drm_info" "mdele"    "info_mat" "negldl"   "delr"     "df"      
## [7] "p_val"    "p_est"    "cdf_est" 
## 
## drmdel> # A brief summary of the DRM fit
## drmdel> summaryDRM(drmfit)
## Basic information about the DRM:
##     Number of samples (m+1): 5 
##     Basis function: 6 
##     Dimension of the baisi function (d): 2 
##     Sample sizes: 100 200 180 150 175 
##     Sample proportions (rho): 0.124 0.248 0.224 0.186 0.217 
## 
## Maximum dual-empirical-likelihood estimator (MDELE):
##    alpha[] beta[,1] beta[,2]
## F1  -22.14  -0.1935     13.5
## F2  -19.49   0.0243     11.4
## F3   -4.95  -4.6389     17.8
## F4  -32.88  -1.2457     22.8
## 
## Default hypothesis testing problem:
##     H0: All distribution functions, F_k's, are equal.
##     H1: One of the distribution functions is different from the others.
## Dual-empirical-likelihood ratio statistc (DELR): 1035.261 
## Degree of freedom: 8 
## p-value: 0 
## 
## Summary statistics of the estimated F_k's (mean, var -- variance, sd -- standard deviation, Q1 -- first quartile, Q3 -- third quartile, IQR -- inter-quartile range):
##     mean  var    sd   Q1    Q3  IQR
## F0  2.64 1.59 1.262 1.75  3.25 1.51
## F1 10.20 7.70 2.776 8.11 11.82 3.71
## F2 10.29 9.31 3.051 7.99 11.96 3.98
## F3  3.61 0.71 0.843 3.06  4.20 1.14
## F4  9.61 4.10 2.024 8.07 10.96 2.89
## 
## drmdel> # Another way of specifying basis function for drmdel() is to pass
## drmdel> # a user-specified R function to the basis_func argument of the
## drmdel> # drmdel() function.
## drmdel> # NOTE: If the basis function one wants to use is included in the
## drmdel> # built-in function list, one should use the built-in functions by
## drmdel> # passing an integer between 1 to 11 to the drmdel() function,
## drmdel> # because the computation will be faster with a built-in function
## drmdel> # than with a user-specified function.
## drmdel> basis_gamma <- function(x) return(c(x, log(abs(x))))
## 
## drmdel> drmfit1 <- drmdel(x=x, n_samples=n_samples,
## drmdel+                   basis_func=basis_gamma)
## 
## drmdel> # One can see the summary of this DRM fit is exactly the same as
## drmdel> # that of the previous fit with basis_func=6
## drmdel> summaryDRM(drmfit1)
## Basic information about the DRM:
##     Number of samples (m+1): 5 
##     Basis function:
## function(x) return(c(x, log(abs(x))))
##     Dimension of the baisi function (d): 2 
##     Sample sizes: 100 200 180 150 175 
##     Sample proportions (rho): 0.124 0.248 0.224 0.186 0.217 
## 
## Maximum dual-empirical-likelihood estimator (MDELE):
##    alpha[] beta[,1] beta[,2]
## F1  -22.14  -0.1935     13.5
## F2  -19.49   0.0243     11.4
## F3   -4.95  -4.6389     17.8
## F4  -32.88  -1.2457     22.8
## 
## Default hypothesis testing problem:
##     H0: All distribution functions, F_k's, are equal.
##     H1: One of the distribution functions is different from the others.
## Dual-empirical-likelihood ratio statistc (DELR): 1035.261 
## Degree of freedom: 8 
## p-value: 0 
## 
## Summary statistics of the estimated F_k's (mean, var -- variance, sd -- standard deviation, Q1 -- first quartile, Q3 -- third quartile, IQR -- inter-quartile range):
##     mean  var    sd   Q1    Q3  IQR
## F0  2.64 1.59 1.262 1.75  3.25 1.51
## F1 10.20 7.70 2.776 8.11 11.82 3.71
## F2 10.29 9.31 3.051 7.99 11.96 3.98
## F3  3.61 0.71 0.843 3.06  4.20 1.14
## F4  9.61 4.10 2.024 8.07 10.96 2.89
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment