version 0.2

DESCRIPTION

@@ -1,17 +1,17 @@
 Package: IsingSampler
 Type: Package
-Title: Sampling methods and distribution functions for the Ising model
-Version: 0.1.1
-Date: 2014-01-09
+Title: Sampling Methods and Distribution Functions for the Ising Model
+Version: 0.2
+Date: 2015-03-01
 Author: Sacha Epskamp
 Maintainer: Sacha Epskamp <mail@sachaepskamp.com>
-Description: This package can be used to sample states from the Ising model and compute the probability of states. Sampling can be done for any number of nodes, but due to the intractibility of the Ising model the distribution can only be computed up to ~10 nodes.
+Description: Sample states from the Ising model and compute the probability of states. Sampling can be done for any number of nodes, but due to the intractibility of the Ising model the distribution can only be computed up to ~10 nodes.
 License: GPL-2
-Imports: plyr
-Depends: Rcpp (>= 0.10.4)
+Imports: plyr, magrittr, nnet
+Depends: Rcpp (>= 0.10.4), R (>= 3.0.0)
 LinkingTo: Rcpp
 URL: github.com/SachaEpskamp/IsingSampler
-Packaged: 2014-01-09 12:29:18 UTC; sacha
+Packaged: 2015-03-02 12:27:36 UTC; sacha
 NeedsCompilation: yes
 Repository: CRAN
-Date/Publication: 2014-01-09 13:53:20
+Date/Publication: 2015-03-02 15:51:07

MD5

@@ -1,18 +1,26 @@
-4fd2eb609887430c0f71c7b15562c0fd *DESCRIPTION
-30f81225affa398b447cf8e4076ac928 *NAMESPACE
-f90a17dabbb6e2a1df5a40bafab55e1a *R/Distribution.R
-3b6fee6061b42e322b20f7229f7bc998 *R/IsingSampler.R
-a6fdb6afd0ed6cc9fec0a3d7673ebb2a *R/LinTransform.R
-ccef8d1913c5717170ab0015c9882a63 *R/RcppExports.R
+be6bf10e5bb995d6f5e12b585cefaf28 *DESCRIPTION
+31c5003275bdf678b461ce5cc68fb4d8 *NAMESPACE
+34682408144679642ae77c1d9001a7aa *NEWS
+ee59f5a99ada1c61d659dc8da170226c *R/Distribution.R
+0d235e758b64db12096afd3c54031f49 *R/EstimateIsing.R
+d7703488891d599df35f7ca09078e017 *R/IsingSampler.R
+135a1c7333c258e6da982a501564c10a *R/LinTransform.R
+b972bac52ff7e633b41994914085c85c *R/PseudoLikelihood.R
+59a567e0e50d53b533251b6532573bcc *R/RcppExports.R
+e33860b126069542b68a2833a617bc68 *R/logisticRegressionEstimation.R
+d41d8cd98f00b204e9800998ecf8427e *R/positiveGraph.R
 ff547ababc85108875b8e8f94239d4dd *inst/COPYING
 bb280cb0bbd7ffc5f80b6cb7732c6787 *inst/COPYRIGHTS
-e52cf17250f0070efc2378143cd10eab *man/IsingLikelihood.Rd
-74ab3a37365c61269d5ca3fb7f9565e5 *man/IsingSampler-package.Rd
-120b439ca0b8b8a52907160f044ee006 *man/IsingSampler.Rd
-de4818d49ff79994ccf95cce0f8f3a87 *man/IsingStateProb.Rd
-37b56eb0b3b4e0f70da7c89429259d31 *man/IsingSumLikelihood.Rd
-94312d39f27895146ac287d8417e7b5e *man/LinTransform.Rd
-35cc87a8b0fd086aca008ed435245ea4 *src/IsingCpp_CFTP.cpp
+8f0785a149c9b554db872e7bdb24a663 *man/EstimateIsing.Rd
+d4fdbe68931f9e915e0c1cab4b6cc1f4 *man/IsingLikelihood.Rd
+8747e4758e744717dc404a67814088da *man/IsingPL.Rd
+bbdf0f5ab867bc738ee19f5aa80de865 *man/IsingSampler-package.Rd
+2e4ccc4c2c151f7be07e6cb6757bd0e2 *man/IsingSampler.Rd
+e9d55cdbd7d8faf6cc5262a63198b7ba *man/IsingStateProb.Rd
+d82f8f2d44aa9009eb9b889353826fed *man/IsingSumLikelihood.Rd
+a65816cf9b559f1261228648256d0b14 *man/LinTransform.Rd
+79956d3fd3e6c54ac806b226d1cca04d *src/IsingCpp_CFTP.cpp
 b3f86c20b81b51f5162874977d38784d *src/Makevars
-e654bda2f0cf51c4d53c8c6d713641ae *src/Makevars.win
-dc8913cce798e45dfeb464719748238f *src/RcppExports.cpp
+73505498722ad4b07b0ddb3067e92523 *src/Makevars.win
+8fd74081fe8ecc6a5dd24ad79fe3c2a8 *src/PseudoLikelihood.cpp
+8f12ee02f4a27a2a9dc43815f40302b5 *src/RcppExports.cpp

NAMESPACE

@@ -1,8 +1,17 @@
-useDynLib(IsingSampler)
-export(IsingSampler)
-export(IsingSumLikelihood)
-export(IsingLikelihood)
-export(IsingStateProb)
-import(plyr)
-import(Rcpp)
-export(LinTransform)
+useDynLib(IsingSampler)
+export(IsingSampler)
+export(IsingSumLikelihood)
+export(IsingLikelihood)
+export(IsingStateProb)
+import(plyr)
+import(Rcpp)
+export(LinTransform)
+import(magrittr)
+export(EstimateIsing)
+export(IsingPL)
+import(nnet)
+export(EstimateIsing)
+export(EstimateIsingPL)
+export(EstimateIsingUni)
+export(EstimateIsingBi)
+export(EstimateIsingLL)

NEWS

@@ -0,0 +1,2 @@
+Changes in version 0.2:
+  o Added methods for estimating the Ising model without regularization, see ?EstimateIsing

R/Distribution.R

@@ -1,45 +1,45 @@
-IsingSumLikelihood <- function(graph, thresholds, beta, responses = c(0L,1L))
-{
-  stopifnot(isSymmetric(graph))  
-  stopifnot(length(responses)==2)
-  if (any(diag(graph)!=0))
-  {
-    diag(graph) <- 0
-    warning("Diagonal set to 0")
-  }
-  N <- nrow(graph)
-  Allstates <- do.call(expand.grid,lapply(1:N,function(x)c(responses[1],responses[2])))
-  P <- exp(- beta * apply(Allstates,1,function(s)H(graph,s,thresholds)))
-  SumScores <- rowSums(1*(Allstates==1))
-
-  df <- ddply(data.frame(Sum = SumScores, P = P),"Sum",summarize,P=sum(P))
-  df$P <- df$P / sum(df$P)
-  return(df)
-}
-
-IsingLikelihood <- function(graph, thresholds, beta, responses = c(0L,1L))
-{
-  stopifnot(isSymmetric(graph))  
-  stopifnot(length(responses)==2)
-  if (any(diag(graph)!=0))
-  {
-    diag(graph) <- 0
-    warning("Diagonal set to 0")
-  }
-  N <- nrow(graph)
-  Allstates <- do.call(expand.grid,lapply(1:N,function(x)c(responses[1],responses[2])))
-  P <- exp(- beta * apply(Allstates,1,function(s)H(graph,s,thresholds)))
-  df <- cbind(Probability = P / sum(P), Allstates)
-  return(df)
-}
-
-IsingStateProb <- function(s,graph,thresholds,beta,responses=c(0L,1L))
-{
-  if (!is.list(s)) s <- list(s)
-  N <- length(s[[1]])
-  Allstates <- do.call(expand.grid,lapply(1:N,function(x)responses))
-  Dist <- exp(- beta * apply(Allstates,1,function(s)H(graph,s,thresholds)))  
-  Z <- sum(Dist)  
-
-  sapply(s, function(x)exp(-beta*H(graph,x,thresholds))/Z)
+IsingSumLikelihood <- function(graph, thresholds, beta, responses = c(0L,1L))
+{
+  stopifnot(isSymmetric(graph))  
+  stopifnot(length(responses)==2)
+  if (any(diag(graph)!=0))
+  {
+    diag(graph) <- 0
+    warning("Diagonal set to 0")
+  }
+  N <- nrow(graph)
+  Allstates <- do.call(expand.grid,lapply(1:N,function(x)c(responses[1],responses[2])))
+  P <- exp(- beta * apply(Allstates,1,function(s)H(graph,s,thresholds)))
+  SumScores <- rowSums(1*(Allstates==1))
+
+  df <- ddply(data.frame(Sum = SumScores, P = P),"Sum",summarize,P=sum(P))
+  df$P <- df$P / sum(df$P)
+  return(df)
+}
+
+IsingLikelihood <- function(graph, thresholds, beta, responses = c(0L,1L))
+{
+  stopifnot(isSymmetric(graph))  
+  stopifnot(length(responses)==2)
+  if (any(diag(graph)!=0))
+  {
+    diag(graph) <- 0
+    warning("Diagonal set to 0")
+  }
+  N <- nrow(graph)
+  Allstates <- do.call(expand.grid,lapply(1:N,function(x)c(responses[1],responses[2])))
+  P <- exp(- beta * apply(Allstates,1,function(s)H(graph,s,thresholds)))
+  df <- cbind(Probability = P / sum(P), Allstates)
+  return(df)
+}
+
+IsingStateProb <- function(s,graph,thresholds,beta,responses=c(0L,1L))
+{
+  if (!is.list(s)) s <- list(s)
+  N <- length(s[[1]])
+  Allstates <- do.call(expand.grid,lapply(1:N,function(x)responses))
+  Dist <- exp(- beta * apply(Allstates,1,function(s)H(graph,s,thresholds)))  
+  Z <- sum(Dist)  
+
+  sapply(s, function(x)exp(-beta*H(graph,x,thresholds))/Z)
 }

R/EstimateIsing.R

@@ -0,0 +1,15 @@
+# Wrapper function for:
+# - pl: pseudolikelihood
+# - uni: univariate logistic regressions
+# - bi: bivariate logistic regressions
+# - ll: Loglinear model
+
+EstimateIsing <- function(data, responses, beta = 1, method = c('pl', 'uni', 'bi', 'll'),...){
+
+  switch(method[[1]],
+        pl = EstimateIsingPL(data, responses, beta, ...)  ,
+        uni = EstimateIsingUni(data, responses, beta, ...)  ,
+        bi = EstimateIsingBi(data, responses, beta, ...),
+        ll = EstimateIsingLL(data, responses, beta, ...)
+         )
+}

R/IsingSampler.R

@@ -1,57 +1,63 @@
-# Wrapper function:
-IsingSampler <- function(n, graph, thresholds, beta = 1, nIter = 100, responses = c(0L, 1L), method = c("MH","CFTP","direct"), CFTPretry = 10)
-{
-  stopifnot(!missing(graph)|!missing(thresholds))
-  stopifnot(isSymmetric(graph))  
-  stopifnot(length(responses)==2)
-  if (any(diag(graph)!=0))
-  {
-    diag(graph) <- 0
-    warning("Diagonal set to 0")
-  }
-  method <- method[1]
-  if (! method %in% c("MH","CFTP","direct")) stop("method must be 'MH', 'CFTP', or 'direct'")
-
-  if (method %in% c("MH","CFTP"))
-  {
-    try <- 1
-
-    repeat{
-      Res <- IsingSamplerCpp(as.integer(n), graph, thresholds, beta, as.integer(nIter), as.integer(responses), as.logical(method == "CFTP"))                 
-
-      if (any(is.na(Res)) & method == "CFTP")
-      {
-        if (try < CFTPretry)
-        {
-          cat("\n Restarting CFTP chain, attempt: ",try,"\n")
-          try <- try + 1
-        } else 
-        {
-          warning(paste("NA's detected, this means that no CFTP sample was drawn after 100 couplings to the past and",CFTPretry,"resets of the chain. Use higher nIter value or method='MH' for inexact sampling."))
-          break
-        }
-      } else break
-    }
-  } else 
-  {
-    Res <- IsingDir(n, graph, thresholds, beta, responses)
-  }
-
-  return(Res)
-}
-
-## direct sampling function:
-IsingDir <- function(n, graph, thresholds, beta,responses = c(0L,1L))
-{
-  stopifnot(isSymmetric(graph))  
-  stopifnot(length(responses)==2)
-  if (any(diag(graph)!=0))
-  {
-    diag(graph) <- 0
-    warning("Diagonal set to 0")
-  }
-  N <- nrow(graph)
-  Allstates <- do.call(expand.grid,lapply(1:N,function(x)c(responses[1],responses[2])))
-  P <- exp(- beta * apply(Allstates,1,function(s)H(graph,s,thresholds)))  
-  return(Allstates[sample(1:nrow(Allstates),n,TRUE,prob=P),])
+# Wrapper function:
+IsingSampler <- function(n, graph, thresholds, beta = 1, nIter = 100, responses = c(0L, 1L), method = c("MH","CFTP","direct"), CFTPretry = 10, constrain)
+{
+  stopifnot(!missing(graph)|!missing(thresholds))
+  stopifnot(isSymmetric(graph))  
+  stopifnot(length(responses)==2)
+  if (any(diag(graph)!=0))
+  {
+    diag(graph) <- 0
+    warning("Diagonal set to 0")
+  }
+  method <- method[1]
+  if (! method %in% c("MH","CFTP","direct")) stop("method must be 'MH', 'CFTP', or 'direct'")
+
+  # Check constrains:
+  if (missing(constrain))
+  {
+    constrain <- matrix(NA,n,ncol(graph))
+  }
+
+  if (method %in% c("MH","CFTP"))
+  {
+    try <- 1
+
+    repeat{
+      Res <- IsingSamplerCpp(as.integer(n), graph, thresholds, beta, as.integer(nIter), as.integer(responses), as.logical(method == "CFTP"), constrain)                 
+
+      if (any(is.na(Res)) & method == "CFTP")
+      {
+        if (try < CFTPretry)
+        {
+          cat("\n Restarting CFTP chain, attempt: ",try,"\n")
+          try <- try + 1
+        } else 
+        {
+          warning(paste("NA's detected, this means that no CFTP sample was drawn after 100 couplings to the past and",CFTPretry,"resets of the chain. Use higher nIter value or method='MH' for inexact sampling."))
+          break
+        }
+      } else break
+    }
+  } else 
+  {
+    Res <- IsingDir(n, graph, thresholds, beta, responses)
+  }
+
+  return(Res)
+}
+
+## direct sampling function:
+IsingDir <- function(n, graph, thresholds, beta,responses = c(0L,1L))
+{
+  stopifnot(isSymmetric(graph))  
+  stopifnot(length(responses)==2)
+  if (any(diag(graph)!=0))
+  {
+    diag(graph) <- 0
+    warning("Diagonal set to 0")
+  }
+  N <- nrow(graph)
+  Allstates <- do.call(expand.grid,lapply(1:N,function(x)c(responses[1],responses[2])))
+  P <- exp(- beta * apply(Allstates,1,function(s)H(graph,s,thresholds)))  
+  return(Allstates[sample(1:nrow(Allstates),n,TRUE,prob=P),])
 }

R/LinTransform.R

@@ -1,16 +1,30 @@
-LinTransform <- function(graph, thresholds, from = c(0L, 1L), to = c(-1L, 1L), a, b)
-{
-  stopifnot(!missing(graph) & !missing(thresholds))
-
-  if (missing(a) & missing(b))
-  {  
-    a <- (to[1]-to[2])/(from[1]-from[2])    
-    b <- to[1] - a*from[1]
-  }
-
-  diag(graph) <- 0
-
-  return(list(
-    graph = graph/(a^2), 
-    thresholds = thresholds/a - (b*rowSums(graph))/a^2))
+binaryTransform <- function(x, from, to){
+  if (missing(from)){
+    from <- sort(unique(c(unlist(x))))
+  }
+  if (length(from) != 2){
+    stop("Binary data required")
+  }
+
+  x2 <- x
+  x2[x==from[1]] <- to[1]
+  x2[x==from[2]] <- to[2]
+  x2
+}
+
+LinTransform <- function(graph, thresholds, from = c(0L, 1L), to = c(-1L, 1L), a, b)
+{
+  stopifnot(!missing(graph) & !missing(thresholds))
+
+  if (missing(a) & missing(b))
+  {  
+    a <- (to[1]-to[2])/(from[1]-from[2])    
+    b <- to[1] - a*from[1]
+  }
+
+  diag(graph) <- 0
+
+  return(list(
+    graph = graph/(a^2), 
+    thresholds = thresholds/a - (b*rowSums(graph))/a^2))
 }