# R/Bayesian_Bricks.r In bbricks: Bayesian Methods and Graphical Model Structures for Statistical Modeling

#' @title Create objects of type '"BasicBayesian"'.
#' @description
#' A Basic Bayesian Object is with following conditional dependency structure:
#'      \deqn{theta|gamma \sim H(gamma)}
#'      \deqn{X|theta \sim F(theta)}
#' Where H(gamma) is usually called "the prior distribution", F(theta) is called "the observation distribution". Objects of type "LinearGaussianGaussian", "GaussianGaussian", "GaussianInvWishart", "GaussianNIW", "GaussianNIG", "CatDirichlet" and "CatDP" are all "BasicBayesian"s.
#'
#' @seealso \code{\link{LinearGaussianGaussian}} for Linear Gaussian and Gaussian conjugate structure, \code{\link{GaussianGaussian}} for Gaussian-Gaussian conjugate structure, \code{\link{GaussianInvWishart}} for Gaussian-Inverse-Wishart conjugate structure, \code{\link{GaussianNIW}} for Gaussian-NIW conjugate structure, \code{\link{GaussianNIG}} for Gaussian-NIG conjugate structure, \code{\link{CatDirichlet}} for Categorical-Dirichlet conjugate structure, \code{\link{CatDP}} for Categorical-DP conjugate structure ...
#' @param ENV The environment where you want to create the BasicBayesian object
#' @return An object of class "BasicBayesian".
#' @export
BasicBayesian <- function(ENV=parent.frame()){
object <- new.env(parent = ENV)

object$gamma <- list() object$H <- character(0)
object$F <- character(0) class(object) <- c("BasicBayesian","BayesianBrick") return(object) } #' @title Get sample sufficient statistics #' @description #' This is a generic function that will generate the sufficient statistics of a given Bayesian bricks object. i.e. \cr #' for the model structure: #' \deqn{theta|gamma \sim H(gamma)} #' \deqn{x|theta \sim F(theta)} #' get the sufficient statistics T(x). \cr #' For a given sample set x, each row of x is an observation, and a Bayesian bricks object obj. \code{sufficientStatistics()} return the sufficient statistics for different model structures: #' \subsection{class(obj)="LinearGaussianGaussian"}{ #' \deqn{x \sim Gaussian(A z + b, Sigma)} #' \deqn{z \sim Gaussian(m,S)} #' The sufficient statistics are: #' \itemize{ #' \item SA = \eqn{sum_{i=1:N} A_i^T Sigma^{-1} A_i} #' \item SAx = \eqn{sum_{i=1:N} A_i^T Sigma^{-1} (x_i-b_i)} #' } #' See \code{?sufficientStatistics.LinearGaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{mu \sim Gaussian(m,S)} #' Sigma is known. #' The sufficient statistics are: #' \itemize{ #' \item N: the effective number of samples. #' \item xsum: the row sums of the samples. #' } #' See \code{?sufficientStatistics.GaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianInvWishart"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' mu is known. \cr #' The sufficient statistics are: #' \itemize{ #' \item N: the effective number of samples. #' \item xsum: the sample scatter matrix centered on the mean vector. #' } #' See \code{?sufficientStatistics.GaussianInvWishart} for details. #' } #' \subsection{class(obj)="GaussianNIW"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' \deqn{mu \sim Gaussian(m,Sigma/k)} #' The sufficient statistics are: #' \itemize{ #' \item N: the effective number of samples. #' \item xsum: the row sums of the samples. #' \item S: the uncentered sample scatter matrix. #' } #' See \code{?sufficientStatistics.GaussianNIW} for details. #' } #' \subsection{class(obj)="GaussianNIG"}{ #' Where #' \deqn{x \sim Gaussian(X beta,sigma^2)} #' \deqn{sigma^2 \sim InvGamma(a,b)} #' \deqn{beta \sim Gaussian(m,sigma^2 V)} #' X is a row vector, or a design matrix where each row is an obervation. #' The sufficient statistics are: #' \itemize{ #' \item N: the effective number of samples. #' \item SXx: covariance of X and x #' \item SX: the uncentered sample scatter matrix. #' \item Sx: the variance of x #' } #' See \code{?sufficientStatistics.GaussianNIG} for details. #' } #' \subsection{class(obj)="CatDirichlet"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim Dirichlet(alpha)} #' The sufficient statistics of CatDirichlet object can either be x itself, or the counts of the unique labels in x.\cr #' See \code{?sufficientStatistics.CatDirichlet} for details. #' } #' \subsection{class(obj)="CatDP"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim DirichletProcess(alpha)} #' The sufficient statistics of CatDP object can either be x itself, or the counts of the unique labels in x.\cr #' See \code{?sufficientStatistics.CatDP} for details. #' } #' \subsection{class(obj)="DP"}{ #' Where #' \deqn{pi|alpha \sim DP(alpha,U)} #' \deqn{z|pi \sim Categorical(pi)} #' \deqn{theta_z|psi \sim H0(psi)} #' \deqn{x|theta_z,z \sim F(theta_z)} #' The sufficient statistics of "DP" object is the same sufficient statistics of the "BasicBayesian" inside the "DP". #' See \code{?sufficientStatistics.DP} for details. #' } #' \subsection{class(obj)="HDP"}{ #' Where #' \deqn{G|gamma \sim DP(gamma,U)} #' \deqn{pi_j|G,alpha \sim DP(alpha,G), j = 1:J} #' \deqn{z|pi_j \sim Categorical(pi_j)} #' \deqn{k|z,G \sim Categorical(G),\textrm{ if z is a sample from the base measure G}} #' \deqn{theta_k|psi \sim H0(psi)} #' The sufficient statistics of "HDP" object is the same sufficient statistics of the "BasicBayesian" inside the "HDP". #' See \code{?sufficientStatistics.HDP} for details. #' } #' \subsection{class(obj)="HDP2"}{ #' Where #' \deqn{G |eta \sim DP(eta,U)} #' \deqn{G_m|gamma,G \sim DP(gamma,G), m = 1:M} #' \deqn{pi_{mj}|G_m,alpha \sim DP(alpha,G_m), j = 1:J_m} #' \deqn{z|pi_{mj} \sim Categorical(pi_{mj})} #' \deqn{k|z,G_m \sim Categorical(G_m),\textrm{ if z is a sample from the base measure} G_m} #' \deqn{u|k,G \sim Categorical(G),\textrm{ if k is a sample from the base measure} G} #' \deqn{theta_u|psi \sim H0(psi)} #' \deqn{x|theta_u,u \sim F(theta_u)} #' The sufficient statistics of "HDP2" object is the same sufficient statistics of the "BasicBayesian" inside the "HDP2". #' See \code{?sufficientStatistics.HDP2} for details. #' } #' #' @seealso \code{\link{sufficientStatistics.LinearGaussianGaussian}} for Linear Gaussian and Gaussian conjugate structure, \code{\link{sufficientStatistics.GaussianGaussian}} for Gaussian-Gaussian conjugate structure, \code{\link{sufficientStatistics.GaussianInvWishart}} for Gaussian-Inverse-Wishart conjugate structure, \code{\link{sufficientStatistics.GaussianNIW}} for Gaussian-NIW conjugate structure, \code{\link{sufficientStatistics.GaussianNIG}} for Gaussian-NIG conjugate structure, \code{\link{sufficientStatistics.CatDirichlet}} for Categorical-Dirichlet conjugate structure, \code{\link{sufficientStatistics.CatDP}} for Categorical-DP conjugate structure ... #' @param obj a "BayesianBrick" object used to select a method. #' @param x a set of samples. #' @param ... further arguments passed to or from other methods. #' @return An object of corresponding sufficient statistics class, such as "ssGaussian" #' @export #' @examples #' x <- rGaussian(10,mu = 1,Sigma = 1) #' obj <- GaussianNIW() #an GaussianNIW object #' sufficientStatistics(obj=obj,x=x) sufficientStatistics <- function(obj,x,...) UseMethod("sufficientStatistics") #' @title Get weighted sample sufficient statistics #' @description #' This is a generic function that will generate the weighted sufficient statistics of a given "BayesianBrick" object. That is, for the model structure: #' \deqn{theta|gamma \sim H(gamma)} #' \deqn{x|theta \sim F(theta)} #' get the weighted sufficient statistics T(x). #' For a given sample set x, each row of x is an observation, the sample weights w, and a Bayesian bricks object obj. \code{sufficientStatistics_Weighted()} return the weighted sufficient statistics for different model structures: #' \subsection{class(obj)="LinearGaussianGaussian"}{ #' \deqn{x \sim Gaussian(A z + b, Sigma)} #' \deqn{z \sim Gaussian(m,S)} #' The sufficient statistics are: #' \itemize{ #' \item SA = \eqn{sum_{i=1:N} w_i A_i^T Sigma^{-1} A_i} #' \item SAx = \eqn{sum_{i=1:N} w_i A_i^T Sigma^{-1} (x_i-b_i)} #' } #' See \code{?sufficientStatistics.LinearGaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{mu \sim Gaussian(m,S)} #' Sigma is known. #' The sufficient statistics are: #' \itemize{ #' \item N: the effective number of samples. #' \item xsum: the row sums of the samples. #' } #' See \code{?sufficientStatistics_Weighted.GaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianInvWishart"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' mu is known. \cr #' The sufficient statistics are: #' \itemize{ #' \item N: the effective number of samples. #' \item xsum: the sample scatter matrix centered on the mean vector. #' } #' See \code{?sufficientStatistics_Weighted.GaussianInvWishart} for details. #' } #' \subsection{class(obj)="GaussianNIW"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' \deqn{mu \sim Gaussian(m,Sigma/k)} #' The sufficient statistics are: #' \itemize{ #' \item N: the effective number of samples. #' \item xsum: the row sums of the samples. #' \item S: the uncentered sample scatter matrix. #' } #' See \code{?sufficientStatistics_Weighted.GaussianNIW} for details. #' } #' \subsection{class(obj)="GaussianNIG"}{ #' Where #' \deqn{x \sim Gaussian(X beta,sigma^2)} #' \deqn{sigma^2 \sim InvGamma(a,b)} #' \deqn{beta \sim Gaussian(m,sigma^2 V)} #' X is a row vector, or a design matrix where each row is an obervation. #' The sufficient statistics are: #' \itemize{ #' \item N: the effective number of samples. #' \item SXx: covariance of X and x #' \item SX: the uncentered sample scatter matrix. #' \item Sx: the variance of x #' } #' See \code{?sufficientStatistics_Weighted.GaussianNIG} for details. #' } #' \subsection{class(obj)="CatDirichlet"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim Dirichlet(alpha)} #' The sufficient statistics of CatDirichlet object can either be x itself, or the counts of the unique labels in x.\cr #' See \code{?sufficientStatistics_Weighted.CatDirichlet} for details. #' } #' \subsection{class(obj)="CatDP"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim DirichletProcess(alpha)} #' The sufficient statistics of CatDP object can either be x itself, or the counts of the unique labels in x.\cr #' See \code{?sufficientStatistics_Weighted.CatDP} for details. #' } #' \subsection{class(obj)="DP"}{ #' Where #' \deqn{pi|alpha \sim DP(alpha,U)} #' \deqn{z|pi \sim Categorical(pi)} #' \deqn{theta_z|psi \sim H0(psi)} #' \deqn{x|theta_z,z \sim F(theta_z)} #' The sufficient statistics of "DP" object is the same sufficient statistics of the "BasicBayesian" inside the "DP". #' See \code{?sufficientStatistics_Weighted.DP} for details. #' } #' \subsection{class(obj)="HDP"}{ #' Where #' \deqn{G|gamma \sim DP(gamma,U)} #' \deqn{pi_j|G,alpha \sim DP(alpha,G), j = 1:J} #' \deqn{z|pi_j \sim Categorical(pi_j)} #' \deqn{k|z,G \sim Categorical(G),\textrm{ if z is a sample from the base measure G}} #' \deqn{theta_k|psi \sim H0(psi)} #' The sufficient statistics of "HDP" object is the same sufficient statistics of the "BasicBayesian" inside the "HDP". #' See \code{?sufficientStatistics_Weighted.HDP} for details. #' } #' \subsection{class(obj)="HDP2"}{ #' Where #' \deqn{G |eta \sim DP(eta,U)} #' \deqn{G_m|gamma,G \sim DP(gamma,G), m = 1:M} #' \deqn{pi_{mj}|G_m,alpha \sim DP(alpha,G_m), j = 1:J_m} #' \deqn{z|pi_{mj} \sim Categorical(pi_{mj})} #' \deqn{k|z,G_m \sim Categorical(G_m),\textrm{ if z is a sample from the base measure } G_m} #' \deqn{u|k,G \sim Categorical(G),\textrm{ if k is a sample from the base measure } G_m} #' \deqn{theta_u|psi \sim H0(psi)} #' \deqn{x|theta_u,u \sim F(theta_u)} #' The sufficient statistics of "HDP2" object is the same sufficient statistics of the "BasicBayesian" inside the "HDP2". #' See \code{?sufficientStatistics_Weighted.HDP2} for details. #' } #' @seealso \code{\link{sufficientStatistics_Weighted.LinearGaussianGaussian}} for Linear Gaussian and Gaussian conjugate structure, \code{\link{sufficientStatistics_Weighted.GaussianGaussian}} for Gaussian-Gaussian conjugate structure, \code{\link{sufficientStatistics_Weighted.GaussianInvWishart}} for Gaussian-Inverse-Wishart conjugate structure, \code{\link{sufficientStatistics_Weighted.GaussianNIW}} for Gaussian-NIW conjugate structure, \code{\link{sufficientStatistics_Weighted.GaussianNIG}} for Gaussian-NIG conjugate structure, \code{\link{sufficientStatistics_Weighted.CatDirichlet}} for Categorical-Dirichlet conjugate structure, \code{\link{sufficientStatistics_Weighted.CatDP}} for Categorical-DP conjugate structure ... #' @param obj a "BayesianBrick" object used to select a method. #' @param x a set of samples. #' @param w numeric, sample weights. #' @param ... further arguments passed to or from other methods. #' @return An object of corresponding sufficient statistics class, such as "ssGaussian" #' @export #' @examples #' \donttest{ #' x <- rGaussian(10,mu = 1,Sigma = 1) #' w <- runif(10) #' obj <- GaussianNIW() #an GaussianNIW object #' sufficientStatistics_Weighted(obj=obj,x=x,w=w) #' } sufficientStatistics_Weighted <- function(obj,x,w,...) UseMethod("sufficientStatistics_Weighted") #' @title update the prior distribution with sufficient statistics #' @description update the prior distribution with sufficient statistics #' @param obj A "BayesianBrick" object used to select a method. #' @param ... further arguments passed to or from other methods. #' @return None, or an error message if the update fails. posterior_bySufficientStatistics <- function(obj,...) UseMethod("posterior_bySufficientStatistics") #' @title update the prior distribution with sufficient statistics #' @description #' This is a generic function that will update the prior distribution of a "BayesianBrick" object by adding information of the observation's sufficient statistics. i.e. for the model structure: #' \deqn{theta|gamma \sim H(gamma)} #' \deqn{x|theta \sim F(theta)} #' update gamma to gamma_posterior by adding the information of x to gamma.\cr #' For a given sample set x or it's sufficient statistics ss, and a Bayesian bricks object obj, \code{posterior()} will update the posterior parameters in obj for different model structures: #' \subsection{class(obj)="LinearGaussianGaussian"}{ #' \deqn{x \sim Gaussian(A z + b, Sigma)} #' \deqn{z \sim Gaussian(m,S)} #' \code{posterior()} will update m and S in obj. #' See \code{?posterior.LinearGaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{mu \sim Gaussian(m,S)} #' Sigma is known. #' \code{posterior()} will update m and S in obj. #' See \code{?posterior.GaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianInvWishart"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' mu is known. #' \code{posterior()} will update v and S in obj. #' See \code{?posterior.GaussianInvWishart} for details. #' } #' \subsection{class(obj)="GaussianNIW"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' \deqn{mu \sim Gaussian(m,Sigma/k)} #' \code{posterior()} will update m, k, v and S in obj. #' See \code{?posterior.GaussianNIW} for details. #' } #' \subsection{class(obj)="GaussianNIG"}{ #' Where #' \deqn{x \sim Gaussian(X beta,sigma^2)} #' \deqn{sigma^2 \sim InvGamma(a,b)} #' \deqn{beta \sim Gaussian(m,sigma^2 V)} #' X is a row vector, or a design matrix where each row is an obervation. #' \code{posterior()} will update m, V, a and b in obj. #' See \code{?posterior.GaussianNIG} for details. #' } #' \subsection{class(obj)="CatDirichlet"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim Dirichlet(alpha)} #' \code{posterior()} will update alpha in obj. #' See \code{?posterior.CatDirichlet} for details. #' } #' \subsection{class(obj)="CatDP"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim DirichletProcess(alpha)} #' \code{posterior()} will update alpha in obj. #' See \code{?posterior.CatDP} for details. #' } #' \subsection{class(obj)="DP"}{ #' Where #' \deqn{pi|alpha \sim DP(alpha,U)} #' \deqn{z|pi \sim Categorical(pi)} #' \deqn{theta_z|psi \sim H0(psi)} #' \deqn{x|theta_z,z \sim F(theta_z)} #' \code{posterior()} will update alpha and psi in obj. #' See \code{?posterior.DP} for details. #' } #' \subsection{class(obj)="HDP"}{ #' Where #' \deqn{G|gamma \sim DP(gamma,U)} #' \deqn{pi_j|G,alpha \sim DP(alpha,G), j = 1:J} #' \deqn{z|pi_j \sim Categorical(pi_j)} #' \deqn{k|z,G \sim Categorical(G),\textrm{ if z is a sample from the base measure G}} #' \deqn{theta_k|psi \sim H0(psi)} #' \code{posterior()} will update gamma, alpha and psi in obj. #' See \code{?posterior.HDP} for details. #' } #' \subsection{class(obj)="HDP2"}{ #' Where #' \deqn{G |eta \sim DP(eta,U)} #' \deqn{G_m|gamma,G \sim DP(gamma,G), m = 1:M} #' \deqn{pi_{mj}|G_m,alpha \sim DP(alpha,G_m), j = 1:J_m} #' \deqn{z|pi_{mj} \sim Categorical(pi_{mj})} #' \deqn{k|z,G_m \sim Categorical(G_m),\textrm{ if z is a sample from the base measure } G_m} #' \deqn{u|k,G \sim Categorical(G),\textrm{ if k is a sample from the base measure} G} #' \deqn{theta_u|psi \sim H0(psi)} #' \deqn{x|theta_u,u \sim F(theta_u)} #' \code{posterior()} will update eta, gamma, alpha and psi in obj. #' See \code{?posterior.HDP2} for details. #' } #' #' @seealso \code{\link{posterior.LinearGaussianGaussian}} for Linear Gaussian and Gaussian conjugate structure, \code{\link{posterior.GaussianGaussian}} for Gaussian-Gaussian conjugate structure, \code{\link{posterior.GaussianInvWishart}} for Gaussian-Inverse-Wishart conjugate structure, \code{\link{posterior.GaussianNIW}} for Gaussian-NIW conjugate structure, \code{\link{posterior.GaussianNIG}} for Gaussian-NIG conjugate structure, \code{\link{posterior.CatDirichlet}} for Categorical-Dirichlet conjugate structure, \code{\link{posterior.CatDP}} for Categorical-DP conjugate structure ... #' @param obj A "BayesianBrick" object used to select a method. #' @param ... further arguments passed to or from other methods. #' @return None, or an error message if the update fails. #' @export posterior <- function(obj,...) UseMethod("posterior") #' @title update the prior distribution with sufficient statistics #' @description update the prior distribution with sufficient statistics #' @param obj A "BayesianBrick" object used to select a method. #' @param ... other parameters. #' @return None, or an error message if the update fails. posteriorDiscard_bySufficientStatistics <- function(obj,...) UseMethod("posteriorDiscard_bySufficientStatistics") #' @title update the prior distribution with sufficient statistics #' @description #' Contrary to posterior(), posteriorDiscard() a generic function that will update the prior distribution of a "BayesianBrick" object by removing the information provided by the observation's sufficient statistics. i.e. for the model structure: #' \deqn{theta|gamma \sim H(gamma)} #' \deqn{x|theta \sim F(theta)} #' update gamma to gamma_posterior by removing the information of x from gamma. #' For a given sample set x or it's sufficient statistics ss, and a Bayesian bricks object obj, \code{posteriorDiscard()} will update the posterior parameters in obj for different model structures: #' \subsection{class(obj)="LinearGaussianGaussian"}{ #' \deqn{x \sim Gaussian(A z + b, Sigma)} #' \deqn{z \sim Gaussian(m,S)} #' \code{posteriorDiscard()} will update m and S in obj. #' See \code{?posteriorDiscard.LinearGaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{mu \sim Gaussian(m,S)} #' Sigma is known. #' \code{posteriorDiscard()} will update m and S in obj. #' See \code{?posteriorDiscard.GaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianInvWishart"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' mu is known. #' \code{posteriorDiscard()} will update v and S in obj. #' See \code{?posteriorDiscard.GaussianInvWishart} for details. #' } #' \subsection{class(obj)="GaussianNIW"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' \deqn{mu \sim Gaussian(m,Sigma/k)} #' \code{posteriorDiscard()} will update m, k, v and S in obj. #' See \code{?posteriorDiscard.GaussianNIW} for details. #' } #' \subsection{class(obj)="GaussianNIG"}{ #' Where #' \deqn{x \sim Gaussian(X beta,sigma^2)} #' \deqn{sigma^2 \sim InvGamma(a,b)} #' \deqn{beta \sim Gaussian(m,sigma^2 V)} # X is a row vector, or a design matrix where each row is an obervation. #' \code{posteriorDiscard()} will update m, V, a and b in obj. #' See \code{?posteriorDiscard.GaussianNIG} for details. #' } #' \subsection{class(obj)="CatDirichlet"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim Dirichlet(alpha)} #' \code{posteriorDiscard()} will update alpha in obj. #' See \code{?posteriorDiscard.CatDirichlet} for details. #' } #' \subsection{class(obj)="CatDP"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim DirichletProcess(alpha)} #' \code{posteriorDiscard()} will update alpha in obj. #' See \code{?posteriorDiscard.CatDP} for details. #' } #' \subsection{class(obj)="DP"}{ #' Where #' \deqn{pi|alpha \sim DP(alpha,U)} #' \deqn{z|pi \sim Categorical(pi)} #' \deqn{theta_z|psi \sim H0(psi)} #' \deqn{x|theta_z,z \sim F(theta_z)} #' \code{posteriorDiscard()} will update alpha and psi in obj. #' See \code{?posteriorDiscard.DP} for details. #' } #' \subsection{class(obj)="HDP"}{ #' Where #' \deqn{G|gamma \sim DP(gamma,U)} #' \deqn{pi_j|G,alpha \sim DP(alpha,G), j = 1:J} #' \deqn{z|pi_j \sim Categorical(pi_j)} #' \deqn{k|z,G \sim Categorical(G),\textrm{ if z is a sample from the base measure } G} #' \deqn{theta_k|psi \sim H0(psi)} #' \code{posteriorDiscard()} will update gamma, alpha and psi in obj. #' See \code{?posteriorDiscard.HDP} for details. #' } #' \subsection{class(obj)="HDP2"}{ #' Where #' \deqn{G |eta \sim DP(eta,U)} #' \deqn{G_m|gamma,G \sim DP(gamma,G), m = 1:M} #' \deqn{pi_{mj}|G_m,alpha \sim DP(alpha,G_m), j = 1:J_m} #' \deqn{z|pi_{mj} \sim Categorical(pi_{mj})} #' \deqn{k|z,G_m \sim Categorical(G_m),\textrm{ if z is a sample from the base measure } G_m} #' \deqn{u|k,G \sim Categorical(G),\textrm{ if k is a sample from the base measure } G} #' \deqn{theta_u|psi \sim H0(psi)} #' \deqn{x|theta_u,u \sim F(theta_u)} #' \code{posteriorDiscard()} will update eta, gamma, alpha and psi in obj. #' See \code{?posteriorDiscard.HDP2} for details. #' } #' #' @seealso \code{\link{posteriorDiscard.LinearGaussianGaussian}} for Linear Gaussian and Gaussian conjugate structure, \code{\link{posteriorDiscard.GaussianGaussian}} for Gaussian-Gaussian conjugate structure, \code{\link{posteriorDiscard.GaussianInvWishart}} for Gaussian-Inverse-Wishart conjugate structure, \code{\link{posteriorDiscard.GaussianNIW}} for Gaussian-NIW conjugate structure, \code{\link{posteriorDiscard.GaussianNIG}} for Gaussian-NIG conjugate structure, \code{\link{posteriorDiscard.CatDirichlet}} for Categorical-Dirichlet conjugate structure, \code{\link{posteriorDiscard.CatDP}} for Categorical-DP conjugate structure ... #' @param obj A "BayesianBrick" object used to select a method. #' @param ... further arguments passed to or from other methods. #' @return None, or an error message if the update fails. #' @export posteriorDiscard <- function(obj,...) UseMethod("posteriorDiscard") #' @title Get the Maximum A Posteriori(MAP) estimate of a "BayesianBrick" object #' @description #' This is a generic function that will generate the MAP estimate of a given "BayesianBrick" object. For the model structure: \cr #' \deqn{theta|gamma \sim H(gamma)} #' \deqn{x|theta \sim F(theta)} #' MAP estimate of theta is theta_MAP = argmax_theta p(theta|gamma,x). #' For a given Bayesian bricks object obj, the MAP estimate will be: #' \subsection{class(obj)="LinearGaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(A z + b, Sigma)} #' \deqn{z \sim Gaussian(m,S)} #' \code{MAP()} will return the MAP estimate of z. #' See \code{?MAP.LinearGaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{mu \sim Gaussian(m,S)} #' Sigma is known. #' \code{MAP()} will return the MAP estimate of mu. #' See \code{?MAP.GaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianInvWishart"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' mu is known. #' \code{MAP()} will return the MAP estimate of Sigma. #' See \code{?MAP.GaussianInvWishart} for details. #' } #' \subsection{class(obj)="GaussianNIW"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' \deqn{mu \sim Gaussian(m,Sigma/k)} #' \code{MAP()} will return the MAP estimate of mu and Sigma. #' See \code{?MAP.GaussianNIW} for details. #' } #' \subsection{class(obj)="GaussianNIG"}{ #' Where #' \deqn{x \sim Gaussian(X beta,sigma^2)} #' \deqn{sigma^2 \sim InvGamma(a,b)} #' \deqn{beta \sim Gaussian(m,sigma^2 V)} #' X is a row vector, or a design matrix where each row is an obervation. #' \code{MAP()} will return the MAP estimate of beta and sigma^2. #' See \code{?MAP.GaussianNIG} for details. #' } #' \subsection{class(obj)="CatDirichlet"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim Dirichlet(alpha)} #' \code{MAP()} will return the MAP estimate of pi. #' See \code{?MAP.CatDirichlet} for details. #' } #' \subsection{class(obj)="CatDP"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim DirichletProcess(alpha)} #' \code{MAP()} will return the MAP estimate of pi. #' See \code{?MAP.CatDP} for details. #' } #' @seealso \code{\link{MAP.LinearGaussianGaussian}} for Linear Gaussian and Gaussian conjugate structure, \code{\link{MAP.GaussianGaussian}} for Gaussian-Gaussian conjugate structure, \code{\link{MAP.GaussianInvWishart}} for Gaussian-Inverse-Wishart conjugate structure, \code{\link{MAP.GaussianNIW}} for Gaussian-NIW conjugate structure, \code{\link{MAP.GaussianNIG}} for Gaussian-NIG conjugate structure, \code{\link{MAP.CatDirichlet}} for Categorical-Dirichlet conjugate structure, \code{\link{MAP.CatDP}} for Categorical-DP conjugate structure ... #' @param obj A "BayesianBrick" object used to select a method. #' @param ... further arguments passed to or from other methods. #' @return A list of the MAP estimates #' @export MAP <- function(obj,...) UseMethod("MAP") #' @title Get the Mean Posterior Estimate(MPE) of a "BayesianBrick" object #' @description #' This is a generic function that will generate the MPE estimate of a given "BayesianBrick" object. i.e. for the model structure: \cr #' \deqn{theta|gamma \sim H(gamma)} #' \deqn{x|theta \sim F(theta)} #' MPE estimate of theta is theta_MPE = E(theta|gamma,x), E() is the expectation function. #' For a given Bayesian bricks object obj, the MPE estimate will be: #' \subsection{class(obj)="LinearGaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(A z + b, Sigma)} #' \deqn{z \sim Gaussian(m,S)} #' \code{MPE()} will return the MPE estimate of z. #' See \code{?MPE.LinearGaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{mu \sim Gaussian(m,S)} #' Sigma is known. #' \code{MPE()} will return the MPE estimate of mu. #' See \code{?MPE.GaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianInvWishart"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' mu is known. #' \code{MPE()} will return the MPE estimate of Sigma. #' See \code{?MPE.GaussianInvWishart} for details. #' } #' \subsection{class(obj)="GaussianNIW"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' \deqn{mu \sim Gaussian(m,Sigma/k)} #' \code{MPE()} will return the MPE estimate of mu and Sigma. #' See \code{?MPE.GaussianNIW} for details. #' } #' \subsection{class(obj)="GaussianNIG"}{ #' Where #' \deqn{x \sim Gaussian(X beta,sigma^2)} #' \deqn{sigma^2 \sim InvGamma(a,b)} #' \deqn{beta \sim Gaussian(m,sigma^2 V)} #' X is a row vector, or a design matrix where each row is an obervation. #' \code{MPE()} will return the MPE estimate of beta and sigma^2. #' See \code{?MPE.GaussianNIG} for details. #' } #' \subsection{class(obj)="CatDirichlet"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim Dirichlet(alpha)} #' \code{MPE()} will return the MPE estimate of pi. #' See \code{?MPE.CatDirichlet} for details. #' } #' \subsection{class(obj)="CatDP"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim DirichletProcess(alpha)} #' \code{MPE()} will return the MPE estimate of pi. #' See \code{?MPE.CatDP} for details. #' } #' @seealso \code{\link{MPE.LinearGaussianGaussian}} for Linear Gaussian and Gaussian conjugate structure, \code{\link{MPE.GaussianGaussian}} for Gaussian-Gaussian conjugate structure, \code{\link{MPE.GaussianInvWishart}} for Gaussian-Inverse-Wishart conjugate structure, \code{\link{MPE.GaussianNIW}} for Gaussian-NIW conjugate structure, \code{\link{MPE.GaussianNIG}} for Gaussian-NIG conjugate structure, \code{\link{MPE.CatDirichlet}} for Categorical-Dirichlet conjugate structure, \code{\link{MPE.CatDP}} for Categorical-DP conjugate structure ... #' @param obj A "BayesianBrick" object used to select a method. #' @param ... further arguments passed to or from other methods. #' @return A list of MPE estimates #' @export MPE <- function(obj,...) UseMethod("MPE") #' @title Get the density from the posterior distribution. #' @description #' This is a generic function that will generate the the density value of the posterior distribution. i.e. for the model structure: \cr #' \deqn{theta|gamma \sim H(gamma)} #' \deqn{x|theta \sim F(theta)} #' get the probability density/mass from the distribution \eqn{theta \sim H(gamma)}. #' For a given Bayesian bricks object obj and an observation of theta, \code{dPosterior()} will calculate the density value for different model structures: #' \subsection{class(obj)="LinearGaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(A z + b, Sigma)} #' \deqn{z \sim Gaussian(m,S)} #' \code{dPosterior()} will return p(theta|m,S) #' See \code{?dPosterior.LinearGaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{mu \sim Gaussian(m,S)} #' Sigma is known. #' \code{dPosterior()} will return p(mu|m,S) #' See \code{?dPosterior.GaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianInvWishart"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' mu is known. #' \code{dPosterior()} will return p(Sigma|v,S) #' See \code{?dPosterior.GaussianInvWishart} for details. #' } #' \subsection{class(obj)="GaussianNIW"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' \deqn{mu \sim Gaussian(m,Sigma/k)} #' \code{dPosterior()} will return p(mu,Sigma|m,k,v,S) #' See \code{?dPosterior.GaussianNIW} for details. #' } #' \subsection{class(obj)="GaussianNIG"}{ #' Where #' \deqn{x \sim Gaussian(X beta,sigma^2)} #' \deqn{sigma^2 \sim InvGamma(a,b)} #' \deqn{beta \sim Gaussian(m,sigma^2 V)} #' X is a row vector, or a design matrix where each row is an obervation. #' \code{dPosterior()} will return p(beta,sigma^2|m,V,a,b) #' See \code{?dPosterior.GaussianNIG} for details. #' } #' \subsection{class(obj)="CatDirichlet"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim Dirichlet(alpha)} #' \code{dPosterior()} will return p(pi|alpha) #' See \code{?dPosterior.CatDirichlet} for details. #' } #' @seealso \code{\link{dPosterior.LinearGaussianGaussian}} for Linear Gaussian and Gaussian conjugate structure, \code{\link{dPosterior.GaussianGaussian}} for Gaussian-Gaussian conjugate structure, \code{\link{dPosterior.GaussianInvWishart}} for Gaussian-Inverse-Wishart conjugate structure, \code{\link{dPosterior.GaussianNIW}} for Gaussian-NIW conjugate structure, \code{\link{dPosterior.GaussianNIG}} for Gaussian-NIG conjugate structure, \code{\link{dPosterior.CatDirichlet}} for Categorical-Dirichlet conjugate structure ... #' @param obj A "BayesianBrick" object used to select a method. #' @param ... further arguments passed to or from other methods. #' @return numeric, the density value #' @export dPosterior <- function(obj,...) UseMethod("dPosterior") #' @title Generate random samples from the posterior distribution #' @description #' This is a generic function that will generate random samples from the posterior distribution. i.e. for the model structure: \cr #' \deqn{theta|gamma \sim H(gamma)} #' \deqn{x|theta \sim F(theta)} #' generate random sampels of theta from the distribution \eqn{theta \sim H(gamma)}. #' For a given Bayesian bricks object obj, \code{rPosterior()} will generate random samples for different model structures: #' \subsection{class(obj)="LinearGaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(A z + b, Sigma)} #' \deqn{z \sim Gaussian(m,S)} #' \code{rPosterior()} will generate random samples from Gaussian(m,S) #' See \code{?rPosterior.LinearGaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{mu \sim Gaussian(m,S)} #' Sigma is known. #' \code{rPosterior()} will generate random samples from Gaussian(m,S) #' See \code{?rPosterior.GaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianInvWishart"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' mu is known. #' \code{rPosterior()} will generate random samples from InvWishart(v,S) #' See \code{?rPosterior.GaussianInvWishart} for details. #' } #' \subsection{class(obj)="GaussianNIW"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' \deqn{mu \sim Gaussian(m,Sigma/k)} #' \code{rPosterior()} will generate random samples from NIW(m,k,v,S) #' See \code{?rPosterior.GaussianNIW} for details. #' } #' \subsection{class(obj)="GaussianNIG"}{ #' Where #' \deqn{x \sim Gaussian(X beta,sigma^2)} #' \deqn{sigma^2 \sim InvGamma(a,b)} #' \deqn{beta \sim Gaussian(m,sigma^2 V)} #' X is a row vector, or a design matrix where each row is an obervation. #' \code{rPosterior()} will generate random samples from NIG(m,V,a,b) #' See \code{?rPosterior.GaussianNIG} for details. #' } #' \subsection{class(obj)="CatDirichlet"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim Dirichlet(alpha)} #' \code{rPosterior()} will generate random samples from Dirichlet(alpha) #' See \code{?rPosterior.CatDirichlet} for details. #' } #' @seealso \code{\link{rPosterior.LinearGaussianGaussian}} for Linear Gaussian and Gaussian conjugate structure, \code{\link{rPosterior.GaussianGaussian}} for Gaussian-Gaussian conjugate structure, \code{\link{rPosterior.GaussianInvWishart}} for Gaussian-Inverse-Wishart conjugate structure, \code{\link{rPosterior.GaussianNIW}} for Gaussian-NIW conjugate structure, \code{\link{rPosterior.GaussianNIG}} for Gaussian-NIG conjugate structure, \code{\link{rPosterior.CatDirichlet}} for Categorical-Dirichlet conjugate structure ... #' @param obj A "BayesianBrick" object used to select a method. #' @param ... further arguments passed to or from other methods. #' @return numeric, the density value #' @export rPosterior <- function(obj,...) UseMethod("rPosterior") #' @title Get the marginal likelihood of a "BayesianBrick" object #' @description #' This is a generic function that will generate the marginal likelihood of a set of observations conditioned on a given "BayesianBrick" object. i.e. for the model structure: \cr #' \deqn{theta|gamma \sim H(gamma)} #' \deqn{x|theta \sim F(theta)} #' Marginal likelihood is p(x|gamma), p() is the probability density/mass function for continuous/discrete x. #' For a given Bayesian bricks object obj and a sample set x, \code{marginalLikelihood_bySufficientStatistics()} will calculate the marginal likelihood for different model structures: #' \subsection{class(obj)="LinearGaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(A z + b, Sigma)} #' \deqn{z \sim Gaussian(m,S)} #' \code{marginalLikelihood_bySufficientStatistics()} will return p(x|m,S,A,b,Sigma) #' See \code{?marginalLikelihood_bySufficientStatistics.LinearGaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{mu \sim Gaussian(m,S)} #' Sigma is known. #' \code{marginalLikelihood_bySufficientStatistics()} will return p(x|m,S,Sigma) #' See \code{?marginalLikelihood_bySufficientStatistics.GaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianInvWishart"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' mu is known. #' \code{marginalLikelihood_bySufficientStatistics()} will return p(x|mu,v,S) #' See \code{?marginalLikelihood_bySufficientStatistics.GaussianInvWishart} for details. #' } #' \subsection{class(obj)="GaussianNIW"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' \deqn{mu \sim Gaussian(m,Sigma/k)} #' \code{marginalLikelihood_bySufficientStatistics()} will return p(x|m,k,v,S) #' See \code{?marginalLikelihood_bySufficientStatistics.GaussianNIW} for details. #' } #' \subsection{class(obj)="GaussianNIG"}{ #' Where #' \deqn{x \sim Gaussian(X beta,sigma^2)} #' \deqn{sigma^2 \sim InvGamma(a,b)} #' \deqn{beta \sim Gaussian(m,sigma^2 V)} #' X is a row vector, or a design matrix where each row is an obervation. #' \code{marginalLikelihood_bySufficientStatistics()} will return p(x,X|m,V,a,b) #' See \code{?marginalLikelihood_bySufficientStatistics.GaussianNIG} for details. #' } #' \subsection{class(obj)="CatDirichlet"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim Dirichlet(alpha)} #' \code{marginalLikelihood_bySufficientStatistics()} will return p(x|alpha) #' See \code{?marginalLikelihood_bySufficientStatistics.CatDirichlet} for details. #' } #' \subsection{class(obj)="CatDP"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim DirichletProcess(alpha)} #' \code{marginalLikelihood_bySufficientStatistics()} will return p(x|alpha) #' See \code{?marginalLikelihood_bySufficientStatistics.CatDP} for details. #' } #' @seealso \code{\link{marginalLikelihood_bySufficientStatistics.LinearGaussianGaussian}} for Linear Gaussian and Gaussian conjugate structure, \code{\link{marginalLikelihood_bySufficientStatistics.GaussianGaussian}} for Gaussian-Gaussian conjugate structure, \code{\link{marginalLikelihood_bySufficientStatistics.GaussianInvWishart}} for Gaussian-Inverse-Wishart conjugate structure, \code{\link{marginalLikelihood_bySufficientStatistics.GaussianNIW}} for Gaussian-NIW conjugate structure, \code{\link{marginalLikelihood_bySufficientStatistics.GaussianNIG}} for Gaussian-NIG conjugate structure, \code{\link{marginalLikelihood_bySufficientStatistics.CatDirichlet}} for Categorical-Dirichlet conjugate structure, \code{\link{marginalLikelihood_bySufficientStatistics.CatDP}} for Categorical-DP conjugate structure ... #' @param obj A "BayesianBrick" object used to select a method. #' @param ... further arguments passed to or from other methods. #' @return numeric, the marginal likelihood #' @export marginalLikelihood_bySufficientStatistics <- function(obj,...) UseMethod("marginalLikelihood_bySufficientStatistics") #' @title Get the marginal likelihood of a "BayesianBrick" object #' @description #' This is a generic function that will generate the marginal likelihood of a set of observations conditioned on a given "BayesianBrick" object. i.e. for the model structure: \cr #' \deqn{theta|gamma \sim H(gamma)} #' \deqn{x|theta \sim F(theta)} #' Marginal likelihood is p(x|gamma), p() is the probability density/mass function for continuous/discrete x. #' For a given Bayesian bricks object obj and a sample set x, \code{marginalLikelihood()} will calculate the marginal likelihood for different model structures: #' \subsection{class(obj)="LinearGaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(A z + b, Sigma)} #' \deqn{z \sim Gaussian(m,S)} #' \code{marginalLikelihood()} will return p(x|m,S,A,b,Sigma) #' See \code{?marginalLikelihood.LinearGaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{mu \sim Gaussian(m,S)} #' Sigma is known. #' \code{marginalLikelihood()} will return p(x|m,S,Sigma) #' See \code{?marginalLikelihood.GaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianInvWishart"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' mu is known. #' \code{marginalLikelihood()} will return p(x|mu,v,S) #' See \code{?marginalLikelihood.GaussianInvWishart} for details. #' } #' \subsection{class(obj)="GaussianNIW"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' \deqn{mu \sim Gaussian(m,Sigma/k)} #' \code{marginalLikelihood()} will return p(x|m,k,v,S) #' See \code{?marginalLikelihood.GaussianNIW} for details. #' } #' \subsection{class(obj)="GaussianNIG"}{ #' Where #' \deqn{x \sim Gaussian(X beta,sigma^2)} #' \deqn{sigma^2 \sim InvGamma(a,b)} #' \deqn{beta \sim Gaussian(m,sigma^2 V)} #' X is a row vector, or a design matrix where each row is an obervation. #' \code{marginalLikelihood()} will return p(x,X|m,V,a,b) #' See \code{?marginalLikelihood.GaussianNIG} for details. #' } #' \subsection{class(obj)="CatDirichlet"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim Dirichlet(alpha)} #' \code{marginalLikelihood()} will return p(x|alpha) #' See \code{?marginalLikelihood.CatDirichlet} for details. #' } #' \subsection{class(obj)="CatDP"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim DirichletProcess(alpha)} #' \code{marginalLikelihood()} will return p(x|alpha) #' See \code{?marginalLikelihood.CatDP} for details. #' } #' @seealso \code{\link{marginalLikelihood.LinearGaussianGaussian}} for Linear Gaussian and Gaussian conjugate structure, \code{\link{marginalLikelihood.GaussianGaussian}} for Gaussian-Gaussian conjugate structure, \code{\link{marginalLikelihood.GaussianInvWishart}} for Gaussian-Inverse-Wishart conjugate structure, \code{\link{marginalLikelihood.GaussianNIW}} for Gaussian-NIW conjugate structure, \code{\link{marginalLikelihood.GaussianNIG}} for Gaussian-NIG conjugate structure, \code{\link{marginalLikelihood.CatDirichlet}} for Categorical-Dirichlet conjugate structure, \code{\link{marginalLikelihood.CatDP}} for Categorical-DP conjugate structure ... #' @param obj A "BayesianBrick" object used to select a method. #' @param ... further arguments passed to or from other methods. #' @return numeric, the marginal likelihood #' @export marginalLikelihood <- function(obj,...) UseMethod("marginalLikelihood") #' @title Get the density value of the posterior predictive distribution #' @description #' This is a generic function that will generate the the density value of the posterior predictive distribution. i.e. for the model structure: \cr #' \deqn{theta|gamma \sim H(gamma)} #' \deqn{x|theta \sim F(theta)} #' get the probability density/mass of the posterior predictive distribution of a new sample x_new: p(x_new|gamma). #' For a given Bayesian bricks object obj and a new sample x, \code{dPosteriorPredictive()} will calculate the marginal likelihood for different model structures: #' \subsection{class(obj)="LinearGaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(A z + b, Sigma)} #' \deqn{z \sim Gaussian(m,S)} #' \code{dPosteriorPredictive()} will return p(x|m,S,A,b,Sigma) #' See \code{?dPosteriorPredictive.LinearGaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{mu \sim Gaussian(m,S)} #' Sigma is known. #' \code{dPosteriorPredictive()} will return p(x|m,S,Sigma) #' See \code{?dPosteriorPredictive.GaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianInvWishart"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' mu is known. #' \code{dPosteriorPredictive()} will return p(x|mu,v,S) #' See \code{?dPosteriorPredictive.GaussianInvWishart} for details. #' } #' \subsection{class(obj)="GaussianNIW"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' \deqn{mu \sim Gaussian(m,Sigma/k)} #' \code{dPosteriorPredictive()} will return p(x|m,k,v,S) #' See \code{?dPosteriorPredictive.GaussianNIW} for details. #' } #' \subsection{class(obj)="GaussianNIG"}{ #' Where #' \deqn{x \sim Gaussian(X beta,sigma^2)} #' \deqn{sigma^2 \sim InvGamma(a,b)} #' \deqn{beta \sim Gaussian(m,sigma^2 V)} #' X is a row vector, or a design matrix where each row is an obervation. #' \code{dPosteriorPredictive()} will return p(x,X|m,V,a,b) #' See \code{?dPosteriorPredictive.GaussianNIG} for details. #' } #' \subsection{class(obj)="CatDirichlet"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim Dirichlet(alpha)} #' \code{dPosteriorPredictive()} will return p(x|alpha) #' See \code{?dPosteriorPredictive.CatDirichlet} for details. #' } #' \subsection{class(obj)="CatDP"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim DirichletProcess(alpha)} #' \code{dPosteriorPredictive()} will return p(x|alpha) #' See \code{?dPosteriorPredictive.CatDP} for details. #' } #' @seealso \code{\link{dPosteriorPredictive.LinearGaussianGaussian}} for Linear Gaussian and Gaussian conjugate structure, \code{\link{dPosteriorPredictive.GaussianGaussian}} for Gaussian-Gaussian conjugate structure, \code{\link{dPosteriorPredictive.GaussianInvWishart}} for Gaussian-Inverse-Wishart conjugate structure, \code{\link{dPosteriorPredictive.GaussianNIW}} for Gaussian-NIW conjugate structure, \code{\link{dPosteriorPredictive.GaussianNIG}} for Gaussian-NIG conjugate structure, \code{\link{dPosteriorPredictive.CatDirichlet}} for Categorical-Dirichlet conjugate structure, \code{\link{dPosteriorPredictive.CatDP}} for Categorical-DP conjugate structure ... #' @param obj A "BayesianBrick" object used to select a method. #' @param ... further arguments passed to or from other methods. #' @return numeric, the density value #' @export dPosteriorPredictive <- function(obj,...) UseMethod("dPosteriorPredictive") #' @title Generate random samples from the posterior predictive distribution #' @description #' This is a generic function that will generate random samples from the posterior predictive distribution. i.e. for the model structure: \cr #' \deqn{theta|gamma \sim H(gamma)} #' \deqn{x|theta \sim F(theta)} #' generate x_new from the posterior predictive distribution of x|gamma. #' For a given Bayesian bricks object obj, \code{rPosteriorPredictive()} will generate random samples from different model structures: #' \subsection{class(obj)="LinearGaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(A z + b, Sigma)} #' \deqn{z \sim Gaussian(m,S)} #' \code{rPosteriorPredictive()} will generate samples from the distribution of x|m,S,A,b,Sigma #' See \code{?rPosteriorPredictive.LinearGaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianGaussian"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{mu \sim Gaussian(m,S)} #' Sigma is known. #' \code{rPosteriorPredictive()} will generate samples from the distribution of x|m,S,Sigma #' See \code{?rPosteriorPredictive.GaussianGaussian} for details. #' } #' \subsection{class(obj)="GaussianInvWishart"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' mu is known. #' \code{rPosteriorPredictive()} will generate samples from the distribution of x|mu,v,S #' See \code{?rPosteriorPredictive.GaussianInvWishart} for details. #' } #' \subsection{class(obj)="GaussianNIW"}{ #' Where #' \deqn{x \sim Gaussian(mu,Sigma)} #' \deqn{Sigma \sim InvWishart(v,S)} #' \deqn{mu \sim Gaussian(m,Sigma/k)} #' \code{rPosteriorPredictive()} will generate samples from the distribution of x|m,k,v,S #' See \code{?rPosteriorPredictive.GaussianNIW} for details. #' } #' \subsection{class(obj)="GaussianNIG"}{ #' Where #' \deqn{x \sim Gaussian(X beta,sigma^2)} #' \deqn{sigma^2 \sim InvGamma(a,b)} #' \deqn{beta \sim Gaussian(m,sigma^2 V)} #' X is a row vector, or a design matrix where each row is an obervation. #' \code{rPosteriorPredictive()} will generate samples from the distribution of x,X|m,V,a,b #' See \code{?rPosteriorPredictive.GaussianNIG} for details. #' } #' \subsection{class(obj)="CatDirichlet"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim Dirichlet(alpha)} #' \code{rPosteriorPredictive()} will generate samples from the distribution of x|alpha #' See \code{?rPosteriorPredictive.CatDirichlet} for details. #' } #' \subsection{class(obj)="CatDP"}{ #' Where #' \deqn{x \sim Categorical(pi)} #' \deqn{pi \sim DirichletProcess(alpha)} #' \code{rPosteriorPredictive()} will generate samples from the distribution of x|alpha #' See \code{?rPosteriorPredictive.CatDP} for details. #' } #' @seealso \code{\link{rPosteriorPredictive.GaussianNIW}} for Gaussian-NIW conjugate structure, \code{\link{rPosteriorPredictive.GaussianNIG}} for Gaussian-NIG conjugate structure, \code{\link{rPosteriorPredictive.CatDirichlet}} for Categorical-Dirichlet conjugate structure, \code{\link{rPosteriorPredictive.CatDP}} for Categorical-DP conjugate structure ... #' @param obj A "BayesianBrick" object used to select a method. #' @param n integer, specify the number of samples to be generated. #' @param ... further arguments passed to or from other methods. #' @return a matrix or vector or list of random samples, depends on the type of 'obj'. #' @export rPosteriorPredictive <- function(obj,n,...) UseMethod("rPosteriorPredictive") #' Print the content of an BasicBayesian object #' @param x An object of the target type. #' @param ... other parameters passed to print. #' @return None. #' @export print.BasicBayesian <- function(x,...){ cat("Prior distribution: ",x$H,"\n")
cat("Observation distribution: ",x$F,"\n") cat("Parameters for prior distribution:\n") print(x$gamma)
}

#' a internal version of "is", only for this package
#' @param object an object to be test with
#' @param class2 the class you want to test
#' @return logical, return TRUE if 'class2' is in 'class(object)'.
#' @noRd
.is <- function(object,class2){
isTRUE(class2 %in% class(object))
}

#' log multivariate gamma function
#' @param a numeric, a numeric vector to be evaluated.
#' @param p numeric, the dimension.
#' @return numeric, the log value of the multivariate gamma function.
#' @noRd
.lmvgamma <- function(a,p){
sapply(a,function(ai){
p*(p-1)/4*log(pi)+ sum(lgamma(ai+(1-(1:p))/2))
},simplify = TRUE)
}

#' determinant of a matrix
#' @param m numeric or matrix. When m is a numeric vector, it's treated as a diagonal matrix.
#' @return logical, return TRUE if 'class2' is in 'class(object)'.
#' @noRd
.det <- function(m){
if(is.matrix(m)){
base::det(m)
}else if(is.vector(m)){
base::prod(m)
}
}

#' @title log sum exp
#' @description For each row l of a matrix x, calculate log(sum(exp(l))).
#' @param x matrix, the values in x a all logged. If x is a numeric vector, it will be converted to a matrix with 1 row.
#' @return numeric, the logsumexp of each row of x.
#' @export
#' @examples
#' \donttest{
#' ## Normalize the rows of x to make them sum up to 1
#' x <- matrix(runif(6,-1000,-20),nrow=2)
#' x <- x-logsumexp(x)
#' x <- exp(x)
#' rowSums(x)
#' }
logsumexp <- function(x){
if(is.vector(x)) x <- matrix(x,nrow =1L)
else if(!is.matrix(x)) x <- as.matrix(x)
x <- as.matrix(x)
a <- apply(x,1,max)
log(rowSums(exp(x-a)))+a
}

#' @title Inverse of a positive definite symmetric matrix
#' @description Use Cholesky decomposition to calculate the inverse of S,  where \eqn{S = A'A}, A is a upper diagonal matrix. inv(S) = inv(A)inv(A)'.
#' @param S a symmetric positive definitive matrix.
#' @param returnUpper logical, return inv(A) if returnUpper=TRUE,return inv(S) if returnUpper=FALSE, default FALSE.
#' @return A matrix, the inverse of "S".
#' @export
#' @examples
#' Sigma = matrix(c(2,1,1,2),2,2)      # some positive definite symmetric matrix
#' pdsInverse(Sigma)                   # get inv(Sigma)
#' pdsInverse(Sigma,returnUpper=TRUE)  # get inv(A), where Sigma=A'A, A is upper triangle
pdsInverse <- function(S,returnUpper=FALSE){
if(missing(S)) stop("'S' not specified!")
if(length(S)==1L){
if(returnUpper==FALSE)
return(1/S)
else
return(sqrt(S))
}else{
A <- chol(S)                              #Cholesky decomposition S=A'A
Ai <- backsolve(A,diag(nrow(S)))          #A_i = A^(-1)
if(returnUpper==FALSE)
return(crossprod(t(Ai)))
else
return(Ai)
}
}

#' @title determinant of a positive definite symmetric matrix
#' @description Use Cholesky decomposition to calculate the determinant of S, where \eqn{S = A'A}, A is a upper diagonal matrix. det(S) = det(A)*det(A).
#' @param S a symmetric positive definitive matrix.
#' @param LOG logical, return log(det(S)) if TRUE, return det(S) if FALSE, default FALSE.
#' @return A matrix, the determinant of "S".
#' @export
#' @examples
#' Sigma = matrix(c(1.3,1,1,2),2,2)      # some positive definite symmetric matrix
#' pdsDeterminant(Sigma)                 # get inv(Sigma)
#' pdsDeterminant(Sigma,LOG=TRUE)        # get inv(Sigma)
#' pdsInverse(Sigma,returnUpper=TRUE)  # get inv(A), where Sigma=A'A, A is upper triangle
pdsDeterminant <- function(S,LOG=FALSE){
if(missing(S)) stop("'S' not specified!")
if(length(S)==1L){
return(S)
}else{
A <- chol(S)                              #Cholesky decomposition S=A'A
if(LOG){
return(sum(log(diag(A))*2))
}else{
return(prod(diag(A))^2)
}
}
}

#' @description This is a list of additional questions you want devtools::release() to ask when releasing your package.
release_questions <- function(){
c("Have you run 'R CMD build PACKAGE_FOLDER' and 'R CMD check PACKAGE_VERSION.tar.gz' before the release?",
"Have you removed underscores '_' in \\text{...}? For exmple \\text{G_m} is not allowed!!",
"Have you set the correct version number?",
"Have you removed the irrelevant code blocks, such as 'if(FALSE){...}', from your R files?",
"Have you add all the files in ./R folder to DESCRIPTION?",
"Have you removed the unfinished examples from the vignette?",