<jupyter_start><jupyter_text>We shall create an oscillating signal with an exponentially decaying component, which will then have some Gaussian noise applied to it. This can be represented using the following equation$$ y\left(t\right) = \alpha\sin\left(\omega t + \phi\right)e^{-2\beta t} $$We shall choose some initial values of $\alpha=10$, $\beta=0.1$; with a angular frequency of $\omega = \frac{20}{2\pi}$.The phase shift $\phi$ is defined as being some random number $-\pi<\phi\leq\pi$ that is applied to the sinusoidal waveform.<jupyter_code>N_smooth = 360
N = 360
C = 1      # Defines length of data in integers of 2*pi
alpha = 10
beta = 0.5
phi = 2*np.pi*np.random.random(1)-np.pi
#phi = 0.5*np.pi
#phi = -2.59929972
print(phi)

t = np.linspace(0, np.pi*C, N)
t_smooth = np.linspace(0, np.pi*C, N_smooth)
frequency = 3
omega = frequency*2*np.pi
print(omega)


def sine(t,alpha,omega,phi):
    return alpha*np.sin(omega*t+phi)
def expn(t,beta):
    return np.exp(-2 * beta * t)

decay = sine(t,alpha,omega,phi) * expn(t,beta)

f_sg = 1
f_true = np.random.randn(N) * f_sg + decay
f_obs_sg = np.ones(N) * 1
f_obs = f_true + np.random.randn(N) * f_obs_sg<jupyter_output>[-2.51836314]
18.84955592153876<jupyter_text>With the noise added to our data, let's plot it to see its effect and see what we will be aiming to make inferences from.<jupyter_code>plt.figure(1, figsize=(8,5))
plt.plot(t, decay, '--', label=('True'))
plt.plot(t, f_obs, alpha=0.6, label=('Observed'))
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Now let's introduce and implement our Stan code.<jupyter_code>code = '''
functions {
    real sine(real omega, real alpha, real t, real phi){
        return alpha * sin(omega*t+phi);
    }
    real expn(real t, real beta){
        return exp(-2 * beta * t);
    }
}
data {
    int N;
    real t[N];
    real f[N];
    real f_unc[N];
}
parameters {
    real<lower = 0> omega;
    real<lower = 0> beta;
    real phi;
    real ln_height;

}
transformed parameters {
    real height;
    real freq;
    height = exp(ln_height);
    freq = omega/(2 * pi());
}
model {

    real mod[N];
    for (i in 1:N){
        mod[i] = sine(omega, height, t[i], phi) * expn(t[i], beta);
    }
    mod ~ normal(f,f_unc);
    freq ~ normal(3, 0.1);
    height ~ normal(10.0, 0.3);
    beta ~ normal(0.5, 0.1);
    phi ~ normal(0, 1.5);
}
'''
import pystan
sm = pystan.StanModel(model_code=code)
stan_data = {'N': len(t),
             't': t, 
             'f': f_obs,
             'f_unc': f_obs_sg}
# start = {'omega': 19}
# start = {'phi': 0}
nchains = 4
fit = sm.sampling(data=stan_data, iter=1000, chains=nchains, init=[start for n in range(nchains)])
print(fit)
#%matplotlib tk
fit.plot()
plt.show()
%matplotlib inline
import corner
rc("font", family="serif", size=12)
data = np.vstack([fit['freq'], np.exp(fit['ln_height']), fit['beta'], fit['phi']]).T
corner.corner(data, labels=[r'$freqency$', r'$alpha$', r'$beta$', r'$phi$'], truths=[frequency, alpha, beta, phi])
plt.show()<jupyter_output><empty_output><jupyter_text>Comments on model failureAs we can see from the corner plot, the model has converged very well for the values of alpha, beta and omega, however the distribution in phi is clearly bimodal. The model always seems to converge correctly if the initial function has a phase shift very close to intervals of $\frac{\pi}{4}$ (ie looks like $\pm\sin$ or $\pm\cos$).At intermediate values of $\phi$ however, the model seems completely inconsistent with regards to whether it converges correctly for the phase. Indeed, running the fit multiple times for the current phase of $\phi = -2.59929972$ seems to usually cause the bimodal peak, however occasionally it works first time.The plot below shows the issue with this, since the mean value of the array is no longer representative of the most likely parameter values. This results in the model being in antiphase to the data.<jupyter_code>plt.plot(t, decay, label='Data')
plt.plot(t_smooth, sine(t_smooth, np.exp(fit['ln_height'].mean()), fit['omega'].mean(),
                 np.mean(fit['phi'])) * expn(t_smooth, fit['beta'].mean()), label='Model')
plt.legend()<jupyter_output><empty_output><jupyter_text>It is clear to me the reason why it is flicking between two values, since they produce the same value in the sine function, in the same way that $\sin\left(\frac{3\pi}{4}\right) = \sin\left(-\frac{\pi}{4}\right)$. This is very roughly demonstrated in the plot below, taking the max and min estimated values for phi, and showing they represent the same function.I tried to fix this issue by adding some limits in the stan code ```real  phi``` however this made the model worse when I ran it, with omega additionally becoming bimodal. (I've uncommented that line from the stan code, but feel free to re-run with it to see what I mean).<jupyter_code>plt.plot(t, sine(t,10,3.18,min(fit['phi'])))
plt.plot(t, sine(t,10,3.18,max(fit['phi'])))
plt.scatter(np.arange(0,len(fit['phi']),1),fit['phi'])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14505624629958555, 'fraction_numerical': 0.02466942964278666, 'mean_word_length': 4.2193614830072095, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '410020', 'license': 'MIT', 'n_tokens_mistral': 1862, 'n_tokens_neox': 1723, 'n_words': 574, 'path': 'Dan_notebooks/development/OscillationsHBM.ipynb', 'repo_name': 'daw538/y4project', 'size': '5067'}