Using exponential_lcdf to fit exponentially distributed residuals - resolved

Modeling
1

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