R Under development (unstable) (2024-09-21 r87186 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> #the cumulant for binomial coded in R
> cumR<-function(eta,ntrials){
+    if(eta <= 0){
+ 	out <- ntrials * log1p(exp(eta))
+ 	}
+ 	
+ 	if(eta>0){
+ 	out <- ntrials * (eta + log1p(exp(-eta)))
+ 	}
+ 	
+ 	out
+ }
> 
> 
> library(glmm)
Loading required package: trust
Loading required package: mvtnorm
Loading required package: Matrix
Loading required package: parallel
Loading required package: doParallel
Loading required package: foreach
Loading required package: iterators
> 
> #choose a positive value for eta and compare to 
> #calc of cum for bernoulli
> eta <- 1
> ntrials <- 1
> this <- cumR(eta, ntrials=ntrials)
> that <- bernoulli.glmm()$cum(eta)
> all.equal(this, that)
[1] TRUE
> 
> #choose a negative value for eta and compare to 
> #calc of cum for bernoulli
> eta <- -1
> this <- cumR(eta, ntrials=ntrials)
> that <- bernoulli.glmm()$cum(eta)
> all.equal(this, that)
[1] TRUE
> 
> #########
> #the first deriv of the cumulant for binom in R
> cpR<-function(eta,ntrials){
+ 
+ 	denom <- 1 + exp(-eta)
+ 	out <- ntrials/denom
+ 	out
+ }
> 
> #choose a couple values for eta and compare to 
> #calc of cp for bernoulli
> eta <- 1
> this <- cpR(eta, ntrials=ntrials)
> that <- bernoulli.glmm()$cp(eta)
> all.equal(this, that)
[1] TRUE
> 
> eta <- -1
> this <- cpR(eta, ntrials=ntrials)
> that <- bernoulli.glmm()$cp(eta)
> all.equal(this, that)
[1] TRUE
> 
> #########
> #the second deriv of the cumulant for binom in R
> cppR<-function(eta,ntrials){
+ 
+ 	first <- 1+exp(-eta)
+ 	second <-  1+exp(eta)
+ 	out <- ntrials/(first*second)
+ 	out
+ }
> 
> #choose a couple values for eta and compare to 
> #calc of cpp for bernoulli
> eta <- 1
> this <- cppR(eta, ntrials=ntrials)
> that <- bernoulli.glmm()$cpp(eta)
> all.equal(this,that)
[1] TRUE
> 
> eta <- -1
> this <- cppR(eta, ntrials=ntrials)
> that <- bernoulli.glmm()$cpp(eta)
> all.equal(this,that)
[1] TRUE
> 
> #########
> # check cumR for ntrials = 5 for both pos, neg eta
> ntrials <- 5
> eta <- -1
> this <- cumR(eta, ntrials)
> that <- ntrials * log(1+exp(eta))
> all.equal(this, that)
[1] TRUE
> eta <- 1
> this <- cumR(eta,ntrials)
> that <- ntrials * (eta+log(1+exp(-eta)))
> all.equal(this, that)
[1] TRUE
> 
> # check cpR for ntrials = 5 for both pos, neg eta
> ntrials <- 5
> eta <- -1
> this <- cpR(eta,ntrials)
> that <- ntrials/(1+exp(-eta))
> all.equal(this, that)
[1] TRUE
> eta <- 1
> this <- cpR(eta,ntrials)
> that <- ntrials/(1+exp(-eta))
> all.equal(this, that)
[1] TRUE
> 
> # check cppR for ntrials = 5 for both pos, neg eta
> ntrials <- 5
> eta <- -1
> this <- cppR(eta, ntrials)
> that <- ntrials/((1+exp(-eta))*(1+exp(eta)))
> all.equal(this, that)
[1] TRUE
> eta <- 1
> this <- cppR(eta, ntrials)
> that <- ntrials/((1+exp(-eta))*(1+exp(eta)))
> all.equal(this, that)
[1] TRUE
> 
> #########
> # check central finite differences for cumR, cpR, cppR
> ntrials <- 5
> eta <- 1
> delta <- .0001
> 
> #going from value to deriv
> this <- cpR(eta, ntrials) * delta
> that <- cumR(eta + delta/2, ntrials) - cumR(eta - delta/2, ntrials)
> all.equal(this, that)
[1] TRUE
> 
> #going from first to second deriv
> this <- cppR(eta, ntrials) * delta
> that <- cpR(eta + delta/2, ntrials) - cpR(eta - delta/2, ntrials)
> all.equal(this, that)
[1] TRUE
> 
> 
> # ntrials*cumulant for bern = cumulant for binomial in R
> eta <- .8
> ntrials <-2
> this <- cumR(eta, ntrials)
> that <- ntrials * bernoulli.glmm()$cum(eta)
> all.equal(this, that)
[1] TRUE
> 
> eta <- -.8
> this <- cumR(eta, ntrials)
> that <- ntrials * bernoulli.glmm()$cum(eta)
> all.equal(this, that)
[1] TRUE
> 
> ntrials <-5
> this <- cumR(eta, ntrials)
> that <- ntrials * bernoulli.glmm()$cum(eta)
> all.equal(this, that)
[1] TRUE
> 
> 
> 
> ##################  
> #I am now confident in the R calcs for binom's cum and its 2 derivs
> ##################  
> 
> 
> ##################  
> #Now I'll check the C calcs for bern against C calcs for bin
> #void cum3(double *etain, int *neta, int *typein, int *ntrials, double *wts, double *cumout)
> #bin = bern for cumulant, ntrials = 1
> eta <- 1
> ntrials <- 1
> wts <- 1
> this <- .C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(3), ntrials=as.integer(ntrials), wts=as.double(wts), out=as.double(1.1))
> that <- .C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(1), ntrials=as.integer(ntrials), wts=as.double(wts), out=as.double(1.1))
> all.equal(this$out, that$out)
[1] TRUE
> 
> #bin = 2 bern for cumulant
> eta <- 1
> ntrials <- 2
> bincum <- .C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(3), ntrials=as.integer(ntrials), wts=as.double(wts), out=as.double(0.0))$out
> berncum <- ntrials* .C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(1), ntrials=as.integer(ntrials), wts=as.double(wts), out=as.double(0.0))$out
> all.equal(bincum, berncum)
[1] TRUE
> doublecheck <- cumR(eta, ntrials)
> all.equal(bincum, doublecheck)
[1] TRUE
> 
> eta <- -1
> ntrials <-2
> bincum <- .C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(3), ntrials=as.integer(ntrials), wts=as.double(wts), out=as.double(0.0))$out
> berncum <- ntrials* .C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(1), ntrials=as.integer(ntrials), wts=as.double(wts), out=as.double(0.0))$out
> all.equal(bincum, berncum)
[1] TRUE
> doublecheck <- cumR(eta, ntrials)
> all.equal(bincum, doublecheck)
[1] TRUE
> 
> 
> #compare cumulant's derivative (R vs C)
> eta <- 1
> ntrials <- 4
> this <- .C(glmm:::C_cp3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(3), ntrials=as.integer(ntrials), out=as.double(0.0))$out
> that <- cpR(eta, ntrials)
> all.equal(this, that)
[1] TRUE
> 
> #compare cumulant's 2nd derivative  (R vs C)
> eta <- 1
> ntrials <- 4
> this <- .C(glmm:::C_cpp3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(3), ntrials=as.integer(ntrials), out=as.double(0.0))$out
> that <- cppR(eta, ntrials)
> all.equal(this, that)
[1] TRUE
> 
> #make sure cumulant still works for eta length greater than 1
> neta<- 6
> eta <- 1:neta
> wts <- rep(1,neta)
> ntrials <- rep(4, times=neta)
> length(eta) == length(ntrials)
[1] TRUE
> length(eta) == neta
[1] TRUE
> length(wts) == neta
[1] TRUE
> this <- .C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(neta), typein=as.integer(3), ntrials=as.integer(ntrials), wts=as.double(wts), out=as.double(0.0))$out
> that <- 0
> for(i in 1:neta){
+ 	that <- that + cumR(eta[i], ntrials[i])
+ }
> all.equal(this, that)
[1] TRUE
> #no segfault errors and the same answer. Yay!
> 
> 
> # one last thing to check about the cumulant calcs
> eta<- 1
> ntrials<- 1
> wts <- 1
> length(eta) == length(ntrials)
[1] TRUE
> length(wts) == length(eta)
[1] TRUE
> that <- bernoulli.glmm()$cum(eta)
> this <- binomial.glmm()$cum(eta, ntrials)
> all.equal(this, that)
[1] TRUE
> that <-  .C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(3), ntrials=as.integer(ntrials), wts=as.double(wts), out=as.double(0.0))$out
> all.equal(this, that)
[1] TRUE
> 
> eta<--2
> ntrials<-1
> that <- bernoulli.glmm()$cum(eta)
> this <- binomial.glmm()$cum(eta, ntrials)
> all.equal(this, that)
[1] TRUE
> 
> eta<--2
> ntrials<-2
> that <- 2*bernoulli.glmm()$cum(eta)
> this <- binomial.glmm()$cum(eta, ntrials)
> all.equal(this, that)
[1] TRUE
> that <-  .C(glmm:::C_cum3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(3), ntrials=as.integer(ntrials), wts=as.double(wts), out=as.double(0.0))$out
> all.equal(this, that)
[1] TRUE
> 
> eta<-1
> ntrials<-1
> that <- bernoulli.glmm()$cp(eta)
> this <- binomial.glmm()$cp(eta, ntrials)
> all.equal(this, that)
[1] TRUE
> 
> eta<-1
> ntrials<-3
> that <- 3*bernoulli.glmm()$cp(eta)
> this <- binomial.glmm()$cp(eta, ntrials)
> all.equal(this, that)
[1] TRUE
> that <-  .C(glmm:::C_cp3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(3), ntrials=as.integer(ntrials), out=as.double(0.0))$out
> all.equal(this, that)
[1] TRUE
> 
> eta<-1
> ntrials<-1
> that <- bernoulli.glmm()$cpp(eta)
> this <- binomial.glmm()$cpp(eta, ntrials)
> all.equal(this, that)
[1] TRUE
> 
> eta<-1
> ntrials<-3
> that <- 3*bernoulli.glmm()$cpp(eta)
> this <- binomial.glmm()$cpp(eta, ntrials)
> all.equal(this, that)
[1] TRUE
> that <-  .C(glmm:::C_cpp3, eta=as.double(eta), neta=as.integer(1), typein=as.integer(3), ntrials=as.integer(ntrials), out=as.double(0.0))$out
> all.equal(this, that)
[1] TRUE
> 
> eta<--1
> ntrials<-1
> that <- bernoulli.glmm()$cpp(eta)
> this <- binomial.glmm()$cpp(eta, ntrials)
> all.equal(this, that)
[1] TRUE
> 
> 
> 
> 
> 
> proc.time()
   user  system elapsed 
   1.25    0.25    1.48