 # Using exponential_lcdf to fit exponentially distributed residuals - resolved

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"))

``````
1 Like