Inference on amplitude & noise in simple periodic pulse model Raw

pulse.r

library(rstan)
rstan_options(auto_write = TRUE)

simple_pulse_curve = function(time,amplitude,rate){
  period = 1/rate
  position = (time%%period)/period
  return(amplitude*dbeta(position,2,2))
}

n = 1e3
time = seq(0,10,length.out = n)
y = simple_pulse_curve(time,2,.5)+rnorm(n,0,.1)
plot(time,y,type='l')


mod = rstan::stan_model(file='pulse.stan')

post = rstan::sampling(
    object = mod
  , data = list(n=length(y),time=time,y=y)
  , chains = 4
  , cores = 4
)
post

pulse.Rproj

Version: 1.0

RestoreWorkspace: No
SaveWorkspace: No
AlwaysSaveHistory: No

EnableCodeIndexing: Yes
UseSpacesForTab: Yes
NumSpacesForTab: 2
Encoding: UTF-8

RnwWeave: Sweave
LaTeX: pdfLaTeX

AutoAppendNewline: Yes
StripTrailingWhitespace: Yes

pulse.stan

functions{
  vector pulse(vector time, real amplitude, real rate){
    vector[num_elements(time)] position ;
    vector[num_elements(time)] beta ;
    real period = 1/rate ;
    for(i in 1:num_elements(time)){
      position[i] = fmod(time[i],period)/period ;
      beta[i] = beta_lpdf(position[i]|2,2) ;
    }
    return amplitude*exp(beta) ;
  }
}
data{
  int n;
  vector[n] y;
  vector[n] time;
}
parameters{
  real<lower=0> amplitude;
  real<lower=0> noise;
}
model{
  amplitude ~ weibull(2,2) ;
  noise ~ normal(0,1) ;
  y ~ normal( pulse(time,amplitude,.5) , noise ) ;
}