Function bodies error with custom likelihood

Modeling
1

Hi - I’m trying to fit a model with the following log-likelihood function:

log[L(b|x)] = nlog(\frac{b+1}{x_{max}^{(b+1)} + x_{min}^{b+1)}} + \sum_{k=1}^{K}c_klogx_k

as described in this paper :https://www.int-res.com/articles/suppl/m636p019_supp1.pdf

My attempt at this in Stan is

functions{
  real paretocounts_lpdf(vector xk, real b, real xmin, real xmax, vector ck){
    return(log((b + 1)/(xmax^(b+1) - xmin^(b+1))*exp(b*xk^ck)));
    }
}
    
data {
	int<lower=0> N;
	real<lower=0> xmin;
	real<lower=0> xmax;
	vector[N] xk;
	vector[N] ck;
}

parameters {
	real b;
}

model {
	b ~ normal(-2, 1);
	target += paretocounts_lpdf(xk|b, xmin, xmax, ck);
}

with data that look something like this:

xk = c(0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 
         0.2, 0.24, 0.25, 0.26, 0.29, 0.33, 0.37, 0.38, 0.41, 0.43, 0.55, 
         0.61, 0.67, 0.68, 0.72, 0.73, 3.2, 4.05)
  
ck = c(3L, 2L, 3L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 3L, 1L, 1L, 1L, 3L, 
           1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L)
N = length(ck)
xmin = min(xk)
xmax = max(xk)

dat = list(N = N,
           xmin = xmin,
           xmax = xmax,
           xk = xk,
           ck = ck)

This returns the following error:

Error in stanc(file = file, model_code = model_code, model_name = model_name, :
0
Semantic error in ‘string’, line 2, column 2 to line 4, column 5:
Function bodies must contain a return statement of correct type in every branch.

I’m quite new to direct coding in Stan (though have used brms a lot). What is the error here and how can I fix it? Thank you.

1 Like
2

It looks like your function’s return statement should return a vector, but the function is declared to return a real. Perhaps you intended to return the sum?

1 Like
3

Thank you! This solved the problem and revealed and additional mistake in the equation I wrote. I needed to sum the vectors first in a for loop. The code below now works.

stan
functions{
  real paretocounts_lpdf(vector xk, real b, real xmin, real xmax, vector ck){
    vector[num_elements(xk)] temp; //array that stores densities
        real sumckxk;
        for(i in 1:num_elements(xk)){  //loop over data and store each density
            temp[i] = ck[i]*log(xk[i]); 
        }
    sumckxk = sum(temp);
    return(sum(ck)*log((b + 1)/(xmax^(b+1) - xmin^(b+1))) + b*sumckxk);
    }
}
    
data {
	int<lower=0> N;
	real<lower=0> xmin;
	real<lower=0> xmax;
	vector[N] xk;
	vector[N] ck;
}

parameters {
	real b;
}

model {
	b ~ normal(-2, 5);
	target += paretocounts_lpdf(xk|b, xmin, xmax, ck);
}
1 Like