IsingSampler-package: Sampling methods and distribution functions for the Ising...

Description Details Author(s) Examples

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.

Details

Package: IsingSampler
Type: Package
Version: 1.0
Date: 2013-08-02
License: What license is it under?

Author(s)

Sacha Epskamp

Maintainer: Sacha Epskamp <[email protected]>

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
## This code compares the different sampling algorithms to the expected
## distribution of states in a tractible number of nodes.

## In the end are examples on how to obtain the distribution.

# Input:
N <- 5 # Number of nodes
nSample <- 5000 # Number of samples

# Ising parameters:
Graph <- matrix(sample(0:1,N^2,TRUE,prob = c(0.7, 0.3)),N,N) * rnorm(N^2)
Graph <- pmax(Graph,t(Graph)) / N
diag(Graph) <- 0
Thresh <- -(rnorm(N)^2)
Beta <- 1

# Response options (0,1 or -1,1):
Resp <- c(0L,1L)

# All posible states:
AllStates <- do.call(expand.grid,lapply(1:N,function(x)Resp))

# Simulate with metropolis:
MetData <- IsingSampler(nSample, Graph, Thresh, Beta, 1000/N, 
  responses = Resp, method = "MH")

# Simulate exact samples (CFTP):
ExData <- IsingSampler(nSample, Graph, Thresh, Beta, 100,
  responses = Resp, method = "CFTP")

# Direct simulation:
DirectData <- IsingSampler(nSample, Graph, Thresh, Beta, method = "direct")

# State distirbutions:
MetDist <- apply(AllStates,1,function(x)sum(colSums(t(MetData) == x)==N))
ExDist <- apply(AllStates,1,function(x)sum(colSums(t(ExData) == x)==N))
DirectDist <- apply(AllStates,1,function(x)sum(colSums(t(DirectData) == x)==N))
ExpDist <- exp(- Beta * apply(AllStates,1,function(s)IsingSampler:::H(Graph,s,Thresh)))  
ExpDist <- ExpDist/sum(ExpDist) * nSample

## Plot to compare distributions:
plot(MetDist, type="l", col="blue", pch=16, xlab="State", ylab="Freq", 
  ylim=c(0,max(MetDist,DirectDist,ExDist)))
points(DirectDist,type="l",col="red",pch=16)
points(ExpDist,type="l",col="green",pch=16)
points(ExDist,type="l",col="purple",pch=16)
legend("topright", col=c("blue","red","purple","green"), 
  legend=c("Metropolis","Direct","Exact","Expected"),lty=1,bty='n')

## Likelihoods:

# Sumscores:
IsingSumLikelihood(Graph, Thresh, Beta, Resp)

# All states:
IsingLikelihood(Graph, Thresh, Beta, Resp)

# Single state:
IsingStateProb(rep(Resp[1],N),Graph, Thresh, Beta, Resp)

Example output

Loading required package: Rcpp
  Sum           P
1   0 0.067954965
2   1 0.230508257
3   2 0.342812227
4   3 0.261425317
5   4 0.091251408
6   5 0.006047826
   Probability Var1 Var2 Var3 Var4 Var5
1  0.067954965    0    0    0    0    0
2  0.041274467    1    0    0    0    0
3  0.066918968    0    1    0    0    0
4  0.053638715    1    1    0    0    0
5  0.005752729    0    0    1    0    0
6  0.003494091    1    0    1    0    0
7  0.005665027    0    1    1    0    0
8  0.004540787    1    1    1    0    0
9  0.056126713    0    0    0    1    0
10 0.038751783    1    0    0    1    0
11 0.055271041    0    1    0    1    0
12 0.050360331    1    1    0    1    0
13 0.004751408    0    0    1    1    0
14 0.003280533    1    0    1    1    0
15 0.004678971    0    1    1    1    0
16 0.004263255    1    1    1    1    0
17 0.060435380    0    0    0    0    1
18 0.036789752    1    0    0    0    1
19 0.080873045    0    1    0    0    1
20 0.064969309    1    1    0    0    1
21 0.005116159    0    0    1    0    1
22 0.003114437    1    0    1    0    1
23 0.006846310    0    1    1    0    1
24 0.005499978    1    1    1    0    1
25 0.058461207    0    0    0    1    1
26 0.040454345    1    0    0    1    1
27 0.078231258    0    1    0    1    1
28 0.071440842    1    1    0    1    1
29 0.004949035    0    0    1    1    1
30 0.003424664    1    0    1    1    1
31 0.006622669    0    1    1    1    1
32 0.006047826    1    1    1    1    1
[1] 0.06795496

IsingSampler documentation built on May 29, 2017, 5:49 p.m.