Commit 578d37bf authored by hang yu's avatar hang yu
Browse files

MCMC functions for tf fitting

parent 12153408
......@@ -570,6 +570,30 @@ def update_par_dict_from_data_backup(freq, G_data, G_err, par_dict0,
########################################
### MCMC functions ###
########################################
def log_likelihood(par, freq, GG, coh, n_avg,
n_zr, n_pr, n_zc, n_pc):
GG_model = sysID.par_to_TF_vect(freq, par,
n_zr, n_pr, n_zc, n_pc)
SNRsq = n_avg * coh / (1.-coh)
log_prob = -0.5 * np.sum( SNRsq * np.abs(GG_model - GG)**2./np.abs(GG)**2.)
return log_prob
def log_probability(par, freq, GG, coh, n_avg,
n_zr, n_pr, n_zc, n_pc,
log_prior_func):
lp = log_prior_func(par,
n_zr, n_pr, n_zc, n_pc)
if not np.isfinite(lp):
return -np.inf
log_prob = lp + log_likelihood(par, freq, GG, coh, n_avg,
n_zr, n_pr, n_zc, n_pc)
return log_prob
########################################
......@@ -707,3 +731,34 @@ def plot_prob_contour(sigma, theta, idx,
ax.plot(x1[0, :], x2[0, :], label=label, **ax_kwargs)
ax.plot(x1[1, :], x2[1, :], **ax_kwargs)
return ax
def myDecimateFunc(xx0, fs0, fs, detrend=True):
yy = xx0.copy()
if detrend:
# so far just include a linear detrend
tt0=np.arange(0., len(xx0)/fs0, 1./fs0)
poly_coef = np.polyfit(tt0, xx0, 1)
yy -= (poly_coef[0]*tt0+poly_coef[1])
down_factor = int(fs0/fs)
fir_aa = sig.firwin(20*down_factor+2, 0.8/down_factor, \
window='blackmanharris')
yy=sig.decimate(yy, down_factor, \
ftype=sig.dlti(fir_aa[1:-1], 1), zero_phase=True)
if detrend:
# add back the linear trend
tt = np.arange(0., len(yy)/fs, 1./fs)
yy += (poly_coef[0]*tt + poly_coef[1])
return yy
def rms_from_psd(psd, freq):
nPt = len(psd)
rms = np.zeros(nPt)
for i in range(nPt):
rms[i] = np.trapz(psd[i:], freq[i:])
rms = np.sqrt(rms)
return rms
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment