Require real scalar return type for probability function

I’m getting started using Stan through PyStan and I’m encountering a problem that is similar to what was discussed here (Require real scalar return type for probability function). I don’t think I’m making the same mistake, however. My understanding from the answers to that issue is that the sampling notation is shorthand for “likelihood of the data given the parameters” and so the data is an implied argument for the likelihood functions. As you can see in my code below I am passing my data (‘OBS’ variable) to ‘foo_lpdf’ so I’m not sure where I’m going wrong. I’d appreciate any clarifications if I misunderstood the answers to the other issue.

EDIT: I think I see where I was I thinking about this incorrectly. I was taking the ‘Y ~ foo(theta);’ line as an assignment of the likelihood value to the Y variable. But the tilda is not an assignment operator. If I switch the sampling line to be ‘OBS ~ chi2(theta);’ then the program runs because it’s the likelihood of the of the data given the parameters, theta. Is this correct?

norm_code = """
    functions {
        real foo_lpdf(matrix OBS, real theta) {
            vector[num_elements(OBS[:,1])] tmp;
            real ln_like;
            for (i in 1:num_elements(OBS[:,1])) {
                tmp[i] = (OBS[i,1])^2/OBS[i,2];
            }
            ln_like = -0.5*sum(tmp);
            return ln_like;
        }
    }

    data {
        int N;
        matrix[N, 2] OBS;
    }

    parameters {
        real<lower=0> theta;
    }

    model {
        vector[N] Y;
        real alpha;
        real beta;
        alpha = 1.0;
        beta = 1.0;
        theta ~ gamma(alpha, beta);
        Y ~ foo(theta);
    }
"""

Yes.

y ~ foo(theta);

is just shorthand for

target += foo_lpdf(y | theta);

that also drops unnecessary constants.