<jupyter_start><jupyter_text>The gamma function in stan is parameterised in a way I'm not used to (wiki says the stan way is common in Bayesian stats).  Let's investigate.${\rm Gamma}(\alpha, \beta) = \frac{\beta^{\alpha}}{\Gamma(\alpha)} x^{\alpha - 1} e^{-\beta x}$.This parameterisation gives a mean value $\frac{\alpha }{\beta}$ and a variance of $\frac{\alpha}{\beta^{2}}$.  $\alpha$ is called the shape and $\beta$ is called the rate.$\alpha$ is the same parameter as the shape $k$ in my familiar notation.  The rate $\beta$ is then the recipricol of the scale, i.e., $\theta = \frac{1}{\beta}$.When fitting a power spectrum the data are the model multiplied by the noise distribution.  If we calculate the data divided by the model we are left with the noise distirbution.<jupyter_code>code = '''
functions {
    real harvey(real f, real a, real b, real c){
        return 0.9*a^2/b/(1.0 + (f/b)^c);
    }
    real gaussian(real f, real numax, real width, real height){
        return height * exp(-0.5 * (f - numax)^2 / width^2);
    }
    real apod(real f, real nyq){
        real x = 3.14 / 2.0 * f / nyq;
        return (sin(x) / x)^2;
    }
}
data {
    int N;
    int M;
    int dof;
    vector[N] f;
    real p[N, M];
    vector[M] numax_est;
    vector[M] white_est;

}
parameters {
    real<lower = 0> numax[M];
    real<lower = 0> white[M];
    real ae;
    real be;
    real ce;
    real de;
    real ak;
    real bk;
    real ck;
    real dk;
    real wk;
    real we;
    real hk;
    real he; 
}
transformed parameters {
    vector[M] a;
    vector[M] b;
    vector[M] c;
    vector[M] d;
    vector[M] w;
    vector[M] h;
    for (j in 1:M){
        a[j] = 10^(ak + numax[j] * ae);
        b[j] = 10^(bk + numax[j] * be);
        c[j] = 10^(ck + numax[j] * ce);
        d[j] = 10^(dk + numax[j] * de);
        w[j] = 10^(wk + numax[j] * we);
        h[j] = 10^(hk + numax[j] * he);
    }
}
model {
    real beta[N, M];
    for (j in 1:M){
        for (i in 1:N){
            beta[i, j] = dof / (apod(f[i], 288.8)
                        * (harvey(f[i], a[j], b[j], 4.0)
                        + harvey(f[i], c[j], d[j], 4.0)
                        + gaussian(f[i], 10^numax[j], w[j], h[j]))
                        + 10^white[j]);
        }
        p[1:N, j] ~ gamma(dof, beta[1:N, j]);
    }
    ak ~ normal(3.45, 0.01); // log10
    ae ~ normal(-0.51, 0.01);
    bk ~ normal(-0.46, 0.01); //log10
    be ~ normal(0.89, 0.01);
    ck ~ normal(3.59, 0.01); // log10
    ce ~ normal(-0.61, 0.01);
    dk ~ normal(0.04, 0.01); //log10
    de ~ normal(0.95, 0.01);
    wk ~ normal(-0.83, 0.005); // log10
    we ~ normal(1.0, 0.01);
    hk ~ normal(6.93, 0.008); // log10
    he ~ normal(-2.18, 0.005);
    numax ~ normal(numax_est, 0.01);
    white ~ normal(white_est, 0.03);
}
'''
sm = pystan.StanModel(model_code=code, model_name='backfit')
dat = {'N': len(freq),
       'M': len(numaxs),
       'dof': 100,
      'f': freq,
      'p': data,
      'numax_est': np.log10(numaxs),
      'white_est': np.log10(whites)}
# Note dof is not actually dof but is in fact the number of bins over which we have smoothed.
fit = sm.sampling(data=dat, iter=3000, chains=1)
fit.plot()
print(fit)
diff = 10**fit['numax'].mean(axis=0) - numaxs
print(diff)
n = 7
fig, ax = plt.subplots(n,n, figsize=[14, 14])
ak = fit['ak'].mean()
ae = fit['ae'].mean()
bk = fit['bk'].mean()
be = fit['be'].mean()
ck = fit['ck'].mean()
ce = fit['ce'].mean()
dk = fit['dk'].mean()
de = fit['de'].mean()
wk = fit['wk'].mean()
we = fit['we'].mean()
hk = fit['hk'].mean()
he = fit['he'].mean()
numax = fit['numax'].mean(axis=0)
white = fit['white'].mean(axis=0)
for i in range(0, n*n):
    fit0 = model(freq, ak, ae, bk, be, ck, ce, dk, de, numax[i], wk, we, hk, he, white[i])
    x = i // n
    y = i % n
    ax[x, y].plot(freq, data[:,i], 'r-')
    ax[x, y].plot(freq, fit0, 'k-')
    ax[x, y].set_xscale('log')
    ax[x, y].set_yscale('log')
    ax[x, y].set_ylim([1e2, 1e6])
    ax[x, y].set_xlim([3 ,288])
    ax[x, y].axvline(10**numax[i], c='g', linestyle='--')
    ax[x, y].axvline(numaxs[i], c='b')
    ax[x, y].set_yticks([])
    ax[x, y].set_xticks([])
print(np.median(10**fit['numax'], axis=0))
print(np.std(10**fit['numax'], axis=0))<jupyter_output>[ 31.4149798   31.45894589  37.57075429  42.47402715  56.88204336
  28.84883866  27.70131284  37.5220494   58.79817734 138.45279579
  76.60970644  31.64296256  25.49177255 101.74826563  33.90250513
  35.72466656  30.90824476  25.43971381  29.18495357  45.87369101
  67.77544892  95.06666128   8.73306362  32.02105109 160.22211534
  46.11934755  18.67576145  65.04358203   5.07789577  20.58497811
  32.38746369  47.84290893  73.80351875  31.52608199  83.70505532
  27.44265419  32.92950357  17.98541613 106.57214918  93.63507675
   7.46651818  82.36621635  90.87221975  31.94371335  62.25468309
  31.70562452  72.66386222  60.39747266  32.01999436]
[0.15331653 0.15120198 0.17683773 0.18938059 0.26963975 0.14371295
 0.13202383 0.15475066 0.30860386 1.65173268 0.43789048 0.1655967
 0.12227564 0.75080011 0.16020358 0.17074256 0.16728922 0.12413313
 0.15791135 0.19888931 0.36218673 0.6340508  0.06630429 0.16632847
 2.32285701 0.21461794 0.12131078 0.37737593 0.04977025 0.11123779
 0.16168469 0.2627[...]
{'fraction_non_alphanumeric': 0.18488194973343489, 'fraction_numerical': 0.18145468392993144, 'mean_word_length': 2.6127922971114166, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '420919', 'license': 'MIT', 'n_tokens_mistral': 2820, 'n_tokens_neox': 2234, 'n_words': 548, 'path': 'ClusterFit_old_data.ipynb', 'repo_name': 'grd349/ClusterUp', 'size': '5252'}