Nothing
#' Internal function to sample latent clusters (for observations)
#'
#' This function samples the latent \emph{z} parameter within the Gibbs sampler.
#' Calls on the \code{rmultinom1} function written in C++. Not for calling
#' directly by users.
#'
#' @param nobs numeric. The total number of rows in the dataset.
#' @param nmaxclust numeric. A single number indicating the maximum number of
#' clusters to test.
#' @param dat A data frame containing only columns of the discretized data
#' streams for all observations.
#' @param ltheta numeric. A vector of log-transformed estimates for parameter
#' theta.
#' @param lphi A list containing log-transformed estimates for each data stream
#' of the phi parameter.
#' @param ndata.types numeric. The number of data streams being analyzed.
#'
#' @return A vector with estimates for \emph{z} for each observation within
#' \code{dat}.
#'
#'@importFrom stats "runif"
#'
sample.z.mixmod=function(nobs, nmaxclust, dat, ltheta, lphi, ndata.types){
lprob=matrix(ltheta,nobs,nmaxclust,byrow=T)
for (i in 1:nmaxclust){
for (j in 1:ndata.types){
#account for NA
tmp=lphi[[j]][i,dat[,j]]
cond=is.na(dat[,j])
tmp[cond]=0
#finish calculation
lprob[,i]=lprob[,i]+tmp
}
}
max1=apply(lprob,1,max)
lprob=lprob-max1
tmp=exp(lprob)
prob=tmp/rowSums(tmp)
z=rmultinom1(prob=prob, randu=runif(nobs))
z+1
}
#-----------------------------------
#' Internal function to sample parameter for truncated stick-breaking prior
#'
#' This function samples the latent \emph{v} parameter within the Gibbs sampler.
#' Not for calling directly by users.
#'
#' @param z A vector of latent cluster estimates provided by
#' \code{\link{sample.z.mixmod}}.
#' @param gamma1 numeric. Hyperparameter for the truncated stick-breaking prior.
#' @param nmaxclust numeric. A single number indicating the maximum number of
#' clusters to test.
#'
#' @return A list with estimates for \emph{v} and \emph{theta} for each of the possible states.
#'
#'
#' @importFrom stats "rbeta"
#'
sample.v.mixmod=function(z, gamma1, nmaxclust){
tmp=table(z)
tmp1=rep(0,nmaxclust)
tmp1[as.numeric(names(tmp))]=tmp
theta=v=rep(NA,nmaxclust)
aux=1
for (i in 1:(nmaxclust-1)){
seq1=(i+1):nmaxclust
v[i]=rbeta(1,tmp1[i]+1,sum(tmp1[seq1])+gamma1)
theta[i]=v[i]*aux
aux=aux*(1-v[i])
}
v[nmaxclust]=1
theta[nmaxclust]=aux
list(theta=theta,v=v)
}
#-----------------------------------
#' Internal function to sample bin estimates for each movement variable
#'
#' Estimates values of \emph{phi} matrix for use in characterizing distributions
#' of the movement variables. Not for calling directly by users.
#'
#' @param alpha numeric. A hyperparameter for the Dirichlet distribution.
#' @param nmaxclust numeric. A single number indicating the maximum number of
#' clusters to test.
#' @param nbins numeric. A vector of the number of bins used to discretize each
#' data stream. These must be in the same order as the columns within
#' \code{dat}.
#' @param ndata.types numeric. The number of data streams being analyzed.
#' @param nmat A list based on \code{\link{SummarizeDat}} C++ function to help
#' with multinomial draws.
#'
#' @return A list of proportion estimates that characterize distributions
#' (bins) for each data stream and possible behavioral state.
#'
sample.phi.mixmod=function(alpha, nmaxclust, nbins, ndata.types, nmat){
phi=list()
for (j in 1:ndata.types){
tmp2=matrix(NA,nmaxclust,nbins[j])
for (i in 1:nmaxclust){
tmp2[i,]=MCMCpack::rdirichlet(1,nmat[[j]][i,]+alpha)
}
phi[[j]]=tmp2
}
phi
}
#------------------------------------
#' Internal function to calculate the log-likelihood for iteration of mixture model
#'
#' Calculates the log-likelihood of the mixture model based on estimates for
#' \emph{theta} and \emph{phi}.
#'
#' @param phi A list of proportion estimates that characterize distributions
#' (bins) for each data stream and possible behavioral state.
#' @param theta numeric. A vector of values that sum to one.
#' @param ndata.types numeric. The number of data streams being analyzed.
#' @param dat A data frame containing only columns of the discretized data
#' streams for all observations.
#' @param nobs numeric. The total number of rows in the dataset.
#' @param nmaxclust numeric. A single number indicating the maximum number of
#' clusters to test.
#'
#' @return A numeric value of the log-likelihood based upon the current values
#' for \emph{phi} and \emph{theta}.
#'
#'
get.llk.mixmod=function(phi, theta, ndata.types, dat, nobs, nmaxclust){
prob=matrix(theta,nobs,nmaxclust,byrow=T)
for (i in 1:nmaxclust){
for (j in 1:ndata.types){
#account for NA
tmp=phi[[j]][i,dat[,j]]
cond=is.na(dat[,j])
tmp[cond]=1
#finish calculation
prob[,i]=prob[,i]*tmp
}
}
sum(log(rowSums(prob)))
}
#------------------------------------
#' Internal function to sample the gamma hyperparameter
#'
#' @param v numeric. A vector of proportions for each of the possible clusters.
#' @param ngroup numeric. The total number of possible clusters.
#' @param gamma.possib numeric. A vector of possible values that gamma can take
#' ranging between 0.1 and 1.
#'
#' @return A single numeric value for gamma that falls within
#' \code{gamma.possib} for calculation of the log-likelihood.
#'
#' @importFrom stats "rmultinom"
#'
sample.gamma.mixmod=function(v, ngroup, gamma.possib){
#calculate the log probability associated with each possible value of gamma
cond=v>0.9999999
v[cond]=0.9999999
ngamma=length(gamma.possib)
soma=sum(log(1-v[-ngroup]))
k=(ngroup-1)*(lgamma(1+gamma.possib)-lgamma(gamma.possib))
res=k+(gamma.possib-1)*soma
#exponentiate and normalize probabilities
res=res-max(res)
res1=exp(res)
res2=res1/sum(res1)
#sample from a categorical distribution
tmp=rmultinom(1,size=1,prob=res2)
ind=which(tmp==1)
gamma.possib[ind]
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.