This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
def khatri_rao(J, X): | |
# as needed to compute random effects model matrices | |
# references: | |
# - http://dx.doi.org/10.18637/jss.v067.i01 | |
# - https://stat.ethz.ch/R-manual/R-patched/library/Matrix/html/KhatriRao.html | |
Z = J[..., :, np.newaxis, :] * X[..., np.newaxis, :, :] | |
return Z.reshape((-1,) + Z.shape[2:]) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
set.seed(23) | |
# Create a positive-definite matrix | |
# a Wishart random variable + the identity matrix should fullfil this. Any positive-definite matrix is invertible. | |
m <- rWishart(.1, 5, diag(5))[,,1] + diag(5) | |
# set the bottom and right to 1 | |
# by this the top left is positive definite, but the complete matrix is singular, i.e. not invertible | |
m[4:5, ] <- 1 | |
m[, 4:5] <- 1 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
sig <- function(b, x) 1 / (1 + exp(-b * x)) | |
loss <- function(y, b, x) (y - (sig(b, x)) %*% x | |
bold <- b <- 1 | |
repeat | |
{ | |
bold <- b | |
b <- b - loss(y, b, x)[1] | |
if (abs(b - bold) < 0.000001) break | |
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# create a random matrix with non-negativ values (can be anything) | |
affin <- matrix(runif(100), 10, 10) | |
# drop self-loops | |
diag(affin) <- 0 | |
# create column stochastic transition matrix | |
trans <- sweep(affin, 2, colSums(affin), "/") | |
# create a matrix initial distribution where every column is one observation | |
p0 <- matrix(runif(10 * 10), nrow=10) | |
# column normalize the guys |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
library(microbenchmark) | |
bernoulli.loglik.derivative <- function(p, dat) | |
{ | |
-(sum(dat) - length(dat) * p) | |
} | |
optim <- function(dat) | |
{ | |
p.hat <- 1 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
## Example code for clustering on a three-component mixture model using the EM-algorithm. | |
### First we load some libraries and define some useful functions | |
library(mvtnorm) | |
library(MASS) | |
# Create a 'true' data set (an easy one) | |
.create.data <- function(n) | |
{ |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
using Gadfly | |
using Distributions | |
function df(x, y, b) | |
sum(- (y - x*b)' * x) | |
end | |
function gd() | |
rnorm = Normal() | |
x = rand(rnorm, 100) |