```{r setup}
require(dirichletprocess)
```

# Mixing Distribution Functions

```{r}
#' Create a Gaussian Mixing Distribution with fixed variance.
#'
#'
#' @param priorParameters The prior parameters for the base measure.
#' @param sigma The fixed variance of the model.
#' @return A mixing distribution object.
#' @export
GaussianFixedVarianceMixtureCreate <- function(priorParameters=c(0,1),
                                               sigma){
  mdobj <- MixingDistribution("normalFixedVariance",
                              priorParameters,
                              "conjugate")
  mdobj$sigma <- sigma
  return(mdobj)
}

#' @export
#' @rdname Likelihood
Likelihood.normalFixedVariance <- function(mdObj, x, theta) {
  as.numeric(dnorm(x, theta[[1]], mdObj$sigma))
}

#' @export
#' @rdname PriorDraw
PriorDraw.normalFixedVariance <- function(mdObj, n = 1) {

  priorParameters <- mdObj$priorParameters

  mu <- rnorm(n, priorParameters[1], mdObj$sigma)
  theta <- list(array(mu, dim = c(1, 1, n)))
  return(theta)
}

#' @export
#' @rdname PosteriorParameters
PosteriorParameters.normalFixedVariance <- function(mdObj, x) {

  priorParameters <- mdObj$priorParameters

  n.x <- length(x)
  sigma <- mdObj$sigma
  mu0 <- priorParameters[1]
  sigma0 <- priorParameters[2]

  sigmaPosterior <- (1/sigma0^2 + n.x/sigma^2) ^ (-1)
  muPosterior <- sigmaPosterior * (mu0/sigma0^2 + sum(x)/sigma^2)
  posteriorParameters <- matrix(c(muPosterior, sigmaPosterior), ncol=2)
  return(posteriorParameters)
}

#' @export
#' @rdname PosteriorDraw
PosteriorDraw.normalFixedVariance <- function(mdObj, x, n = 1, ...) {

  posteriorParameters <- PosteriorParameters(mdObj, x)

  mu <- rnorm(n,
              posteriorParameters[1],
              posteriorParameters[2])
  theta <- list(array(mu, dim = c(1, 1, n)))
  return(theta)
}

#' @export
#' @rdname Predictive
Predictive.normalFixedVariance <- function(mdObj, x) {
  priorParameters <- mdObj$priorParameters
  sigma0 <- priorParameters[[2]]
  sigma <- mdObj$sigma

  predictiveArray <- numeric(length(x))

  for (i in seq_along(x)) {

    posteriorParameters <- PosteriorParameters(mdObj, x[i])

    predictiveArray[i] <- dnorm(x[i],
                                posteriorParameters[1],
                                sigma0^2 + sigma^2)
  }
  return(predictiveArray)
}

```

# Dirichlet Process Functions

```{r}
#' Create a Dirichlet Mixture of the Gaussian Distribution with fixed variance.
#'
#'
#' @param y Data.
#' @param sigma The fixed variance
#' @param g0Priors Base Distribution Priors.
#' @param alphaPriors Prior parameter distributions for the alpha concentration parameter.
#' @return Dirichlet process object
#'
#' @export
DirichletProcessGaussianFixedVariance <- function(y,
                                                  sigma,
                                                  g0Priors = c(0, 1),
                                                  alphaPriors = c(2, 4)) {

  mdobj <- GaussianFixedVarianceMixtureCreate(g0Priors, sigma)
  dpobj <- DirichletProcessCreate(y, mdobj, alphaPriors)
  dpobj <- Initialise(dpobj)
  return(dpobj)
}
```

# Fitting Functions (To Be Completed)

```{r}
UpdateSigma <- function(dpObj){
  
  print("Hitting Sigma Update")
  dpObj
}


Fit.gaussianFixedVariance <- function(dpObj, 
                                      its, 
                                      updatePrior = FALSE, 
                                      progressBar = interactive()){
  
    if (progressBar){
    pb <- txtProgressBar(min=0, max=its, width=50, char="-", style=3)
  }

  alphaChain <- numeric(its)
  likelihoodChain <- numeric(its)
  weightsChain <- vector("list", length = its)
  clusterParametersChain <- vector("list", length = its)
  priorParametersChain <- vector("list", length = its)
  labelsChain <- vector("list", length = its)

  for (i in seq_len(its)) {

    alphaChain[i] <- dpObj$alpha
    weightsChain[[i]] <- dpObj$pointsPerCluster / dpObj$n
    clusterParametersChain[[i]] <- dpObj$clusterParameters
    priorParametersChain[[i]] <- dpObj$mixingDistribution$priorParameters
    labelsChain[[i]] <- dpObj$clusterLabels


    likelihoodChain[i] <- sum(log(LikelihoodDP(dpObj)))

    dpObj <- ClusterComponentUpdate(dpObj)
    dpObj <- ClusterParameterUpdate(dpObj)
    dpObj <- UpdateAlpha(dpObj)
    dpObj <- UpdateSigma(dpObj)

    if (updatePrior) {
      dpObj$mixingDistribution <- PriorParametersUpdate(dpObj$mixingDistribution,
                                                        dpObj$clusterParameters)
    }
    if (progressBar){
      setTxtProgressBar(pb, i)
    }
  }

  dpObj$weights <- dpObj$pointsPerCluster / dpObj$n
  dpObj$alphaChain <- alphaChain
  dpObj$likelihoodChain <- likelihoodChain
  dpObj$weightsChain <- weightsChain
  dpObj$clusterParametersChain <- clusterParametersChain
  dpObj$priorParametersChain <- priorParametersChain
  dpObj$labelsChain <- labelsChain

  if (progressBar) {
    close(pb)
  }
  return(dpObj)
}
```


# Testing the Implementation

```{r}
testData <- c(rnorm(100, 5, 2),
              rnorm(100, -5, 2))
```

```{r}
dp <- DirichletProcessGaussianFixedVariance(testData, 2)

dp <- Fit.gaussianFixedVariance(dp, 10, progressBar = TRUE)
```
