Edit: OK - I can see why this doesn’t work. Not a Stan issue - a statistics confusion issue.

Dear all,

I would like to fit a power-law to an estimated power-spectrum, the residuals

of which are exponentially distributed when defined as spec_hat / spec_true. I

can’t work out why the code below doesn’t estimate the parameters correctly.

A reproducible example in R and cmdstanr is below. Thanks in advance for any

help.

```
## Reproducible example with cmdstanr
## Fitting a power-law to an estimated power spectrum
## I will simulate the power spectral estimates
## Stan code
stan_program <- "
data {
int<lower=0> N;
vector<lower=0>[N] freq;
vector[N] spec;
}
parameters {
real<lower=0> alpha;
real beta;
}
transformed parameters{
vector<lower = 0>[N] s_pow;
vector[N] residuals;
for (i in 1:N){
s_pow[i] = alpha * freq[i]^(-beta);
}
for (i in 1:N){
residuals[i] = (spec[i] / s_pow[i]);
}
}
model {
residuals ~ exponential(1);
// priors
beta ~ normal(0, 1);
alpha ~ normal(0.25, 0.1);
}
"
library(cmdstanr)
f <- write_stan_file(stan_program)
mod <- cmdstan_model(f)
alpha = 0.25
beta = 1
N <- 1e03
freq <- seq(1/N, 1/2, 1/N)
spec <- alpha * freq^(-beta)
# power spectra estimated by Fourier transform have proportional errors that are
# exponentially distributed
spec_hat <- spec * rexp(N/2, 1)
plot(spec~freq, log = "xy", type = "l", ylim = range(c(spec_hat, spec)), col = "red")
lines(freq, spec_hat)
hist(spec_hat / spec, freq = FALSE)
xax <- seq(0, 10, length.out = 100)
lines(xax, dexp(xax, 1), col = "red")
# I would like to fit the power-law model using the exponential PDF
stan_dat <- list(
N = N/2,
freq = freq,
spec = spec_hat
)
stan_mod <- mod$sample(
data = stan_dat,
chains = 4,
parallel_chains = 4,
refresh = 500
)
stan_mod$summary(variables = c("alpha", "beta"))
```