Using a custom likelihood function in pystan

Modeling
1

HI everyone!

I am super new to Stan/ pystan and I try to fit a custom likelihood function to fit one parameter to behavior. From looking online I managed to build the following code. The model compiles and runs, but I think it doesn’t give me what I need because I get a different result using a simple grid search in python. I want to find the parameter k that maximizes the custom function I defined. This is my code:

    functions {
        real softmax_lpdf(real k, int N, vector a, vector b, vector c, vector d) {
            real ps;
            real pl;
            real sum_log_prob;
            vector [N] pred;
            for (n in 1:N) {
                ps = exp(a[n] / (1 + k * c[n]));
                pl = 1 - (exp(a[n] / (1 + k * c[n])));
                
                if (d[n] == 0)    
                   pred[n] = log(ps);
                else if (d[n] == 1) 
                    pred[n] = log(pl);
            }              
            sum_log_prob = sum(pred);
            return sum_log_prob;
        }                        
    }         
     data {
        int<lower=0> N; // number of trials
        vector[N] a;
        vector[N] b;
        vector[N] c;
        vector[N] d;
    }
    parameters {
        real<lower=0,upper=1> k; 
    }
    model {
        target += softmax_lpdf(k | N, a, b, c, d);  
    }   
           
    """

Thank you!

2

Do you know the truth? If you don’t, generate data from your model and fit it.

If you are concerned there are differences in your Stan and Python code, generate data from each, fit to each (so fit python data with python model, Stan data with Stan model), and then fit to each other (fit python data with Stan and Stan data with Python) to figure out if there are inconsistencies.

If you are writing the same code, it’s possible to estimate the parameters you want to estimate, and you’re using the same methods in both places, then this should all work.

3

Thanks! I am writing the same code in python and stan and I get different parameters and likelihoods when I sample my stan model, but I get the same likelihoods when I optimize it. So something there must be right. I just cannot understand what is wrong with the stan script…

4

Optimizing finds the mode, but sampling (nuts;hmc) finds the typical set.

See for example https://arxiv.org/abs/1701.02434

5

It makes sense, thanks! However, if I have a Gaussian posterior, both optimizing and sampling should give similar likelihoods, correct?

6

If it helps, print is a valid Stan function so you can print the arguments of the function in Stan and print the output value and see if they are what you also get from Python.

7

Thanks! I’ll try this.