# Function bodies error with custom likelihood

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

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

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