https://github.com/cran/BDgraph
Tip revision: 30ab69a5d52df4a5bb576d33e109b840362c0e7b authored by Reza Mohammadi on 14 November 2018, 17:30:12 UTC
version 2.53
version 2.53
Tip revision: 30ab69a
rgwish.R
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
# Copyright (C) 2012 - 2018 Reza Mohammadi |
# |
# This file is part of BDgraph package. |
# |
# BDgraph is free software: you can redistribute it and/or modify it under |
# the terms of the GNU General Public License as published by the Free |
# Software Foundation; see <https://cran.r-project.org/web/licenses/GPL-3>. |
# |
# Maintainer: Reza Mohammadi <a.mohammadi@uva.nl> |
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
# Sampling from G-Wishart distribution |
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
rgwish = function( n = 1, adj.g = NULL, b = 3, D = NULL, threshold = 1e-8 )
{
if( b <= 2 ) stop( "For G-Wishart distribution parameter 'b' must be more than 2" )
if( is.null( adj.g ) ) stop( "Adjacency matrix must be determined" )
if( is.matrix( adj.g ) ) G <- unclass( adj.g )
#if( class( adj.g ) == "graph" ) G <- unclass( adj.g )
if( class( adj.g ) == "sim" ) G <- adj.g $ G
if( class( adj.g ) == "bdgraph" ) G <- BDgraph::select( adj.g )
if( class( adj.g ) == "ssgraph" ) G <- BDgraph::select( adj.g )
if( !isSymmetric( G ) )
{
G[ lower.tri( G ) ] <- 0
G <- G + t( G )
}
if( ( sum( G == 0 ) + sum( G == 1 ) ) != ( nrow( G ) ^ 2 ) ) stop( " Element of matrix 'adj.g' must be 0 or 1" )
G <- as.matrix( G )
diag( G ) <- 0
p <- nrow( G )
if( p < 1 ) stop( "'p' must be more than or equal with 1" )
if( is.null( D ) ) D <- diag( p )
if( !isSymmetric( D ) ) stop( "Matrix 'D' must be positive definite matrix." )
if( nrow( D ) != p ) stop( "Dimension of matrix G and D must be the same." )
if( p == 1 )
return( rwish( n = n, p = p, b = b, D = D ) )
if( sum( G ) == ( p * ( p - 1 ) ) )
return( rwish( n = n, p = p, b = b, D = D ) )
Ti = chol( solve( D ) )
K = matrix( 0, p, p )
if( n > 1 )
{
samples = array( 0, c( p, p, n ) )
for( i in 1 : n )
{
result = .C( "rgwish_c", as.integer(G), as.double(Ti), K = as.double(K), as.integer(b), as.integer(p), as.double(threshold), PACKAGE = "BDgraph" )
samples[ , , i ] = matrix( result $ K, p, p )
}
}else{
result = .C( "rgwish_c", as.integer(G), as.double(Ti), K = as.double(K), as.integer(b), as.integer(p), as.double(threshold), PACKAGE = "BDgraph" )
samples = matrix( result $ K, p, p )
}
return( samples )
}