Using custom likelihood in matlabStan

Modeling
1

Hi all,

First of all, I want to congratulate you all for this awesome work. Since I couldn’t find a solution in this forum I wanted to ask you in a new thread. I’ve searched in the manual and tried for hours but could find a satisfying solution.

So far I have implemented in matlab my own likelihood function for a different framework. This gives the logposterior and its gradient for a given theta in the following way:
[LogPost,GradLogPost]=likelihood(theta);

The framework I used before doesn’t have a NUTS implementation, why I try to run my model on Stan in Matlab. My aim is to only use the sampler of Stan, rather than the modeling with it.

Is there a way to import this into STAN in the following way:

#include likelihood.m
parameters {
  real<lower=0,upper=1000> theta[4];
}
model {
  target += likelihood(theta);
}

This solution wouldn’t use the gradient calculation of the matlab-skript, but would be a great start. If there is a solution to even use the gradient calculation, I would even be more happy.
My major problem is the lack of C++ skills and I hope that you can provide a solution.

Edit: The Data would be already a part of likelihood.m, why I wouldn’t need to implement either. I’m only interested in the properties of theta.

Thank you and best regards,
ly92

2

No, that is not possible. However, the MATLAB language is not that different from the Stan language, so you should be able to define a user-defined function in the functions block of the Stan program that does the same thing as your likelihood function in MATLAB.

3

Thank you for the quick response. I tried that, but get at the end a likelihood of 0 all the time. Is there a way I can check if the loglikelihood is calculated correctly for a specific parameter? Something like:

theta_test=[0.1, - 1.5, - 0.7];
lp_test=Stan_likelihood(theta_test)

-2.4

4

In R, you can do rstan::expose_stan_functions("filename.stan") and then call your user-defined Stan functions, but not with other interfaces.

1 Like
5

This case study shows how to define user defined likelihoods and test them in R. You may found this R for Matlab users page useful.

1 Like
6

Thank you for your answer.
I tried defining my own lpdf. I also tried the following source: https://link.springer.com/content/pdf/10.3758%2Fs13428-016-0746-9.pdf

Unfortunately the sampled parameters are uniformly distributed. Can you please tell me, what I do wrong in my code?

functions {
	real 	newlprob_lpdf(vector x, real delta, real beta, real t0, real kTl, real sigma) {
        vector[num_elements(x)] xtemp;
        real t[num_elements(x)];
        real sumlprob;
        vector[num_elements(x)] lprob;
        for (i in 1:num_elements(x)){
            t[i]=(i-1)*0.2;
            if (t0<t[i])
                xtemp[i] = kTl * (exp(-beta*(t[i]-t0)) - exp(-delta*(t[i]-t0))) / (delta-beta);
            else
                xtemp[i] = 0;
        }
        for (i in 1:num_elements(x)){
            lprob[i]=log(2*pi()*sigma^2)+((x[i]-xtemp[i])/sigma)^2;            
            }
		sumlprob=- 1/2 * sum(lprob);
        
		return sumlprob;
	}
}
data {
	int<lower=1> Length;
	vector[Length] meas;
}
parameters {
	real<lower=0.01,upper=100> sigma;
	real<lower=0.00001,upper=100000> kTl;
	real<lower=0.00001,upper=100000> beta;
	real<lower=0.00001,upper=100000> delta;
	real<lower=0.01,upper=10> t0;
}
model{
	meas~newlprob(delta,beta,t0,kTl,sigma);
}
7

If you can install R (e.g. RStudio is easy to install), you could call your newlprob_lpdf from R and check that it returns correct values given different data and parameter values. If you don’t want to install R, then you could add print statements to your model block to print out what your function is returning.

btw, it is not recommend to use <lower=0.00001,upper=100000> type constraints unless you have logical constraints for these. Now these look like arbitrary chosen values. The underlying transformation to unconstrained space makes the initialization and HMC to work badly for these. If these need to be positive, remove the upper constraint and add priors.

8

Thank you, after all the struggle I finally installed RStudio with RStan(the Getting Started walkthrough didn’t work for me, but I managed to run it in a virtual environment), and could find the problem in my function with the help of rstan::expose_stan_functions("filename.stan") .

1 Like