Skip to content

Instantly share code, notes, and snippets.

@dirmeier
dirmeier / khatri_rao.py
Last active May 19, 2020 22:48
Khatri-Rao product for any kind of matrices.
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:])
@dirmeier
dirmeier / mi.R
Last active April 23, 2018 10:04
Matrix inverses
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
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
}
@dirmeier
dirmeier / diffusion.R
Last active October 24, 2017 18:34
Analytical solution to network propagation using the Markov random walk with restart
# 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
@dirmeier
dirmeier / bernoulli_mle.R
Last active October 16, 2017 20:24
Maximum likelihood estimation of the success probability on a Bernoulli experiment in R.
library(microbenchmark)
bernoulli.loglik.derivative <- function(p, dat)
{
-(sum(dat) - length(dat) * p)
}
optim <- function(dat)
{
p.hat <- 1
@dirmeier
dirmeier / gmm_em.R
Last active June 19, 2020 19:03
Fitting of Gaussian mixture models using the EM in R.
## 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)
{
@dirmeier
dirmeier / gradient_descent.jl
Last active December 24, 2021 14:45
Gradient descent example in Julia.
using Gadfly
using Distributions
function df(x, y, b)
sum(- (y - x*b)' * x)
end
function gd()
rnorm = Normal()
x = rand(rnorm, 100)