In wldyddl5510/SolMfd: Implementations of solution manifold algorithms

knitr::opts_chunk$set(echo = TRUE)

SolMfd: Algorithms for a solution manifold

Intended use of the package

Solution manifold is a mathematical concept that can be applied to various statistical applications. This package aims to provide the useful functions to utilize solution manifold algorithms. It is consists of three functions. First one is sampling points from the solution manifold (sol_mfd_points), second one is solving constraint likelihood estimation problem (constraint_likelihood), and the last one is calculating posterior density on solution manifold (post_density_solmfd). More statistical applications, such as integral estimation or solving density ridge problems, could be implemented in a later version.

Installation Instructions

Installing from github (currently available):

# install.packages("devtools")
devtools::install_github("wldyddl5510/SolMfd")

# # If you need vignettes ...
# install.packages(c("knitr", "formatR"))
# devtools::install_github("wldyddl5510/SolMfd", build_vignettes = T)

Installing from CRAN (not implemented yet):

# install.packages("SolMfd")

Examples

library(SolMfd)
set.seed(10) # for consistency
N = 20
# define a target function
phi = function(x) {return(pnorm(2, x[[1]], x[[2]]) - pnorm(-5, x[[1]], x[[2]]) - 0.5)}
d = 2
s = 1
# sampling points from sol_mfd
# using gaussian prior
res_point = sol_mfd_points(N, phi, d, s, gamma = 0.1, prior = "gaussian", mean = c(0, 3), sigma = matrix(c(1, 0.25, 0.25, 1), 2, 2))
res_point_tmp = res_point
# check close to 0
head(apply(res_point_tmp, 1, phi)) # are they on solution manifold? 
plot(res_point, xlab = "mean", ylab = "sigma") # how they are distributed

# constraint likelihood
# num of samples
n = 30
# data distribution
X = rnorm(n, mean = 1.5, sd = 3)
# negative log likelihood
nll = function(theta) {return(-sum(dnorm(X, theta[[1]], theta[[2]], log = TRUE)))}
# constraint
C = function(x) {return(pnorm(2, x[[1]], x[[2]]) - pnorm(-5, x[[1]], x[[2]]) - 0.5)}
theta = runif(2, 1, 3)
theta_updated = constraint_likelihood(nll, C, theta, 1)
# plot the convergences
const_val = apply(theta_updated, 1, C)
plot(x = seq(1, nrow(theta_updated)), const_val, xlab = "step", ylab = "constraint", type = 'o')
nll_val = apply(theta_updated, 1, nll)
plot(x = seq(1, nrow(theta_updated)), nll_val, xlab = "step", ylab = "nll", type = 'o')
plot(theta_updated, xlab = "mean", ylab = "sigma", type = 'o')
lines(theta_updated[seq(3, nrow(theta_updated), by = 2), ],  col = 'red')

# posterior density
k = get("dnorm", mode = 'function')
prob_density = function(x, theta) {return(dnorm(x, mean = theta[[1]], sd = theta[[2]]))}
n = 30
X = rnorm(n, 1.5, 3)
res_with_density = post_density_solmfd(X, prob_density, res_point, k)
res_with_density