<jupyter_start><jupyter_text>fitting the FFDthere's lots of ways to fit a FFD, perhaps none of the good<jupyter_code>from scipy.odr import ODR, Model, Data, RealData
def f(B, x):
    if B[0] > 0:
        return np.inf
#     B[0] = -1 # fix slope? no
    return 10**(B[0]*x + B[1])
def quickfit(x,y,xe,ye):
    
    b00, b10 = -1., 30

    linear = Model(f)
    ok = np.where(np.isfinite(x) & np.isfinite(y) & np.isfinite(xe) & np.isfinite(ye))[0]
    mydata = RealData(x[ok], y[ok], sx=xe[ok], sy=ye[ok])
    myodr = ODR(mydata, linear, beta0=[b00, b10])
    myoutput = myodr.run()

    m, b = myoutput.beta[0], myoutput.beta[1]
    m_err, b_err = myoutput.sd_beta[0], myoutput.sd_beta[1]

    return m,b,m_err,b_err
    
plt.figure(figsize=(9,7))

plt.errorbar(ffd_x, ffd_y, xerr=np.sqrt(ffd_xerr**2+LUMIN_err**2), 
             yerr=ffd_yerr, 
             color='C0', alpha=0.25, linestyle='none', lw=1, marker='.')

plt.errorbar(ffd_xY, ffd_yY, xerr=np.sqrt(ffd_xerrY**2+LUMINY_err**2), 
             yerr=ffd_yerrY, 
             color='C1', alpha=0.25, linestyle='none', lw=1, marker='.')

m,b,m_err,b_err = quickfit(ffd_x, ffd_y, 
                           np.sqrt(ffd_xerr**2+LUMIN_err**2), ffd_yerr)
mY,bY,m_errY,b_errY = quickfit(ffd_xY, ffd_yY, 
                               np.sqrt(ffd_xerrY**2+LUMINY_err**2), ffd_yerrY)

plt.plot(ffd_x, 10**(m*ffd_x + b))
plt.plot(ffd_xY, 10**(mY*ffd_xY + bY))

plt.yscale('log')

plt.xlabel('log Energy (erg)')
plt.ylabel(r'Cumulative Flare Rate (day$^{-1}$)')
# plt.title('GJ 1243 (M4, P$_{rot}$=0.59d)', fontsize=14)

# plt.text(31.5, 5, 'Kepler',color='k')
plt.text(30.5, 0.7,'GJ 1243', color='C0')
plt.text(32.5, 0.7,'YZ CMi', color='C1')
plt.xlim(28.8,34)

plt.title('Fit: ODR')
# plt.savefig('YZCMI_vs_GJ1243.png', dpi=150, bbox_inches='tight', pad_inches=0.25)


print(m,b,m_err,b_err)
print(mY,bY,m_errY,b_errY)
def therng(thething, rng=[15,50,84]):
    pcent = np.percentile(thething, rng)
    q = np.diff(pcent)
    return pcent[1], q[0], q[1]
def FFD_mcmc(ffd_x, ffd_y, ffd_xerr, ffd_yerr, Nmin=10, Nsteps=2000):
    
    def _lnProb(theta, x, y, xerr, yerr):
        m, b, xmin = theta
        model = 10**(m * x + b)
        xmodel = (np.log10(y) - b)/m
        xok = x > xmin
        
        # I hate that Python lets variables creep into scope like this, but OK!
        if (sum(xok) >= Nmin) and (-5 < m < -0) and (0 < b < 100) and ( 0 < xmin < np.max(x)):
            
            # my 1st attempt at reduced chi-sq, trying to approximate X and Y errors
#             chisq = (np.sum(((y[xok] -  model[xok])/yerr[xok] )** 2) + 
#                      np.sum(((x[xok] - xmodel[xok])/xerr[xok] )** 2)) / np.float(np.sum(xok))

            # my 2nd attempt, better i think
            chisq = np.sum(((y[xok] -  model[xok])/yerr[xok])** 2 * 
                           ((x[xok] - xmodel[xok])/xerr[xok])** 2) / np.float(np.sum(xok))

            BIC = chisq + np.size(theta) * np.log(len(x)+(np.float(len(x) - np.sum(xok))))
            return -0.5 * BIC
#             return -0.5 * chisq

        return -np.inf

    # initialization
    xok = (ffd_x >= np.mean(ffd_x))
    fit, cov = np.polyfit(ffd_x[xok], np.log10(ffd_y[xok]), 1, 
                     w=1./((ffd_yerr[xok]**2 + (-1 * ffd_xerr[xok])**2)), cov=True)

    pos = np.array([fit[0], fit[1], np.mean(ffd_x)]) + 1e-4 * np.random.randn(32, 3)
    nwalkers, ndim = pos.shape
    
    sampler = emcee.EnsembleSampler(nwalkers, ndim, _lnProb, 
                                    args=(ffd_x, ffd_y, ffd_xerr, ffd_yerr))
    
    sampler.run_mcmc(pos, Nsteps, progress=True);
    samples = sampler.get_chain()
    flat_samples = sampler.get_chain(discard=int(Nsteps*0.1), thin=15, flat=True)
    
    return flat_samples
flat_samplesY = FFD_mcmc(ffd_xY, ffd_yY, ffd_xerrY, ffd_yerrY, Nsteps=1000)
# flat_samplesY = FFD_mcmc(ffd_xY, ffd_yY, 
#                          np.sqrt(ffd_xerrY**2+(LUMINY_err)**2), ffd_yerrY, Nsteps=1000)
print(np.median(flat_samplesY, axis=0))
plt.figure(figsize=(9,7))
inds = np.random.randint(len(flat_samplesY), size=100)
for ind in inds:
    sample = flat_samplesY[ind]
    plt.plot([sample[2], max(ffd_xY)], 10**(sample[0]*np.array([sample[2], max(ffd_xY)]) + sample[1]),
             "k", alpha=0.1)
    
plt.errorbar(ffd_xY, ffd_yY, xerr=np.sqrt(ffd_xerrY**2+LUMINY_err**2), 
             yerr=ffd_yerrY, 
             color='C1', alpha=0.25, linestyle='none', lw=1, marker='.')
plt.plot(ffd_xY[ffd_xY >= np.median(flat_samplesY[:,2])], 
         10**(np.median(flat_samplesY[:,0]) * ffd_xY[ffd_xY >= np.median(flat_samplesY[:,2])] + np.median(flat_samplesY[:,1])), 
         c='r', lw=3)

plt.yscale('log')

plt.xlim(29.5,33.5)
plt.ylim(3e-2, 1e1)
plt.xlabel('log Energy (erg)')
plt.ylabel(r'Cumulative Flare Rate (day$^{-1}$)')

print(therng(flat_samplesY[:,0]))
print(therng(flat_samplesY[:,1]))
print(therng(flat_samplesY[:,2]))
# flat_samples = FFD_mcmc(ffd_x, ffd_y, np.sqrt(ffd_xerr**2+LUMIN_err**2), ffd_yerr, Nsteps=1000)
flat_samples = FFD_mcmc(ffd_x, ffd_y, ffd_xerr, ffd_yerr, Nsteps=1000)
print(np.median(flat_samples, axis=0))
plt.figure(figsize=(9,7))

inds = np.random.randint(len(flat_samples), size=100)
for ind in inds:
    sample = flat_samples[ind]
    plt.plot([sample[2], max(ffd_x)], 10**(sample[0]*np.array([sample[2], max(ffd_x)]) + sample[1]), "k", alpha=0.1)

    
plt.errorbar(ffd_x, ffd_y, xerr=np.sqrt(ffd_xerr**2+LUMIN_err**2), 
             yerr=ffd_yerr, 
             color='C0', alpha=0.25, linestyle='none', lw=1, marker='.')
plt.plot(ffd_x[ffd_x >= np.median(flat_samples[:,2])], 
         10**(np.median(flat_samples[:,0]) * ffd_x[ffd_x >= np.median(flat_samples[:,2])] + np.median(flat_samples[:,1])), 
         c='r', lw=3)

plt.yscale('log')
plt.xlim(29.5,33.5)
plt.ylim(3e-2, 1e1)

print(therng(flat_samples[:,0]))
print(therng(flat_samples[:,1]))
print(therng(flat_samples[:,2]))
plt.figure(figsize=(9,7))
    
plt.errorbar(ffd_x, ffd_y, xerr=ffd_xerr,#np.sqrt(ffd_xerr**2+LUMIN_err**2), 
             yerr=ffd_yerr, 
             color='C0', alpha=0.25, linestyle='none', lw=1, marker='.')
   
plt.errorbar(ffd_xY, ffd_yY, xerr=ffd_xerrY,#np.sqrt(ffd_xerrY**2+LUMINY_err**2), 
             yerr=ffd_yerrY, 
             color='C1', alpha=0.25, linestyle='none', lw=1, marker='.')
plt.plot(ffd_x[ffd_x >= np.median(flat_samples[:,2])], 
         10**(np.median(flat_samples[:,0]) * ffd_x[ffd_x >= np.median(flat_samples[:,2])] + np.median(flat_samples[:,1])), 
         c='r', lw=3)

plt.plot(ffd_xY[ffd_xY >= np.median(flat_samplesY[:,2])], 
         10**(np.median(flat_samplesY[:,0]) * ffd_xY[ffd_xY >= np.median(flat_samplesY[:,2])] + np.median(flat_samplesY[:,1])), 
         c='r', lw=3)

plt.yscale('log')
plt.xlim(29.5,33.5)
plt.ylim(3e-2, 1e1)
def nrpt(x3, y3, m, b):
    ''' find the nearest point on a line '''
    # https://stackoverflow.com/a/47198877
    x1 = 0.
    x2 = x3
    y1 = x1*m + b
    y2 = x2*m + b
    
    dx, dy = x2-x1, y2-y1
    det = dx*dx + dy*dy
    a = (dy*(y3-y1) + dx*(x3-x1)) / det
    return x1+a*dx, y1+a*dy

xok = (ffd_x >= 32)
fit = np.polyfit(ffd_x[xok], ffd_ylog[xok], 1)

print(fit)
# fit = [-1, -1]
plt.errorbar(ffd_x[xok], ffd_ylog[xok], xerr=ffd_xerr[xok], yerr=ffd_yerrlog[xok])
# plt.scatter(ffd_x, ffd_ylog, c='C1')
plt.plot(ffd_x[xok], np.polyval(fit, ffd_x[xok]), c='C2')

xnew, ynew = nrpt(ffd_x, ffd_ylog, *fit)
plt.scatter(xnew[xok], ynew[xok], c='C3')

normalchisq = np.sum( ((ffd_ylog[xok] - np.polyval(fit, ffd_x[xok]))/ffd_yerrlog[xok])**2 ) / np.float(np.sum(xok))

dist = np.sqrt((ffd_x[xok]-xnew[xok])**2 + (ffd_ylog[xok]-ynew[xok])**2)
errs = np.sqrt((ffd_xerr[xok])**2 + (ffd_yerrlog[xok])**2)
newchisq = np.sum( (dist/errs)**2 ) / np.float(np.sum(xok))

print(normalchisq, newchisq)
def _lnProb(theta, x, y, xerr, yerr):
    Nmin = 10
    m, b, xmin = theta
    xok = x > xmin

    # implement some strict "priors" here
    if ((sum(xok) >= Nmin) and (-2 < m < -0.5) and (20 < b < 40) and 
        (np.min(x) < xmin < np.max(x))):

        # orthogonal distance of points to the line
        xnew, ynew = nrpt(x, y, m, b)

        # the York2004 version of the chisq
#             chisq = np.sum((y[xok] - model[xok])**2 / 
#                            (yerr[xok]**2 + (m * xerr[xok])**2)) / np.float(np.sum(xok))

#             chisq = np.sum(((x[xok]-xnew[xok])/xerr[xok])**2 *
#                            ((y[xok]-ynew[xok])/yerr[xok])**2 ) / np.float(np.sum(xok))
        # orthog distance chisq 
        dist = np.sqrt((x[xok]-xnew[xok])**2 + (y[xok]-ynew[xok])**2)
        errs = np.sqrt((xerr[xok])**2 + (yerr[xok])**2)
        chisq = np.sum( (dist/errs)**2 ) / np.float(np.sum(xok))

        BIC = chisq + np.size(theta) * np.log(len(x)+(np.float(len(x) - np.sum(xok))))
        return -0.5 * BIC

    return -np.inf
def FFD_odr(ffd_x, ffd_y, ffd_xerr, ffd_yerr, Nsteps=2000):
    
    # initialization
    xok = (ffd_x >= np.mean(ffd_x))
    fit = np.polyfit(ffd_x[xok], ffd_y[xok], 1)

    pos = np.array([fit[0], fit[1], np.min(ffd_x)]) 
    pos = pos + 1e-4 * np.random.randn(2**5, pos.size)
    nwalkers, ndim = pos.shape
    
    sampler = emcee.EnsembleSampler(nwalkers, ndim, _lnProb, 
                                    args=(ffd_x, ffd_y, ffd_xerr, ffd_yerr))
    
    sampler.run_mcmc(pos, Nsteps, progress=True);
    samples = sampler.get_chain()
    flat_samples = sampler.get_chain(discard=int(Nsteps*0.1), thin=15, flat=True)
    
    return flat_samples
# Our errors seem too small (chisq << 0 near the best model)

# so this actually seems to work oddly well by down-weighting the errors. Huh.
fudge = 7.
flat_samplesY2 = FFD_odr(ffd_xY, ffd_yYlog, ffd_xerrY/fudge, ffd_yerrYlog/fudge,  Nsteps=1000)
plt.figure(figsize=(9,7))
inds = np.random.randint(len(flat_samplesY2), size=100)
for ind in inds:
    sample = flat_samplesY2[ind]
    plt.plot([sample[2], max(ffd_xY)], (sample[0]*np.array([sample[2], max(ffd_xY)]) + sample[1]),
             "k", alpha=0.1)
    
plt.errorbar(ffd_xY, ffd_yYlog, xerr=ffd_xerrY,#np.sqrt(ffd_xerrY**2+LUMINY_err**2), 
             yerr=ffd_yerrYlog, 
             color='C1', alpha=0.25, linestyle='none', lw=1, marker='.')

plt.plot(ffd_xY[ffd_xY >= np.median(flat_samplesY2[:,2])], 
         (np.median(flat_samplesY2[:,0]) * ffd_xY[ffd_xY >= np.median(flat_samplesY2[:,2])] + 
          np.median(flat_samplesY2[:,1])), 
         c='r', lw=3)

plt.xlim(29.5,33.5)
plt.ylim(-2, 1)
plt.xlabel('log Energy (erg)')
plt.ylabel(r'log Cumulative Flare Rate (day$^{-1}$)')

print(therng(flat_samplesY2[:,0]))
print(therng(flat_samplesY2[:,1]))
print(therng(flat_samplesY2[:,2]))
flat_samples2 = FFD_odr(ffd_x[:-30], ffd_ylog[:-30], ffd_xerr[:-30]/fudge, ffd_yerrlog[:-30]/fudge, Nsteps=1000)
plt.figure(figsize=(9,7))
inds = np.random.randint(len(flat_samples2), size=100)
for ind in inds:
    sample = flat_samples2[ind]
    plt.plot([sample[2], max(ffd_x)], (sample[0]*np.array([sample[2], max(ffd_x)]) + sample[1]),
             "k", alpha=0.1)
    
plt.errorbar(ffd_x, ffd_ylog, xerr=ffd_xerr, 
             yerr=ffd_yerrlog, 
             color='C0', alpha=0.25, linestyle='none', lw=1, marker='.')

plt.plot(ffd_x[ffd_x >= np.median(flat_samples2[:,2])], 
         (np.median(flat_samples2[:,0]) * ffd_x[ffd_x >= np.median(flat_samples2[:,2])] + 
          np.median(flat_samples2[:,1])), 
         c='r', lw=3)

plt.xlim(29.5,33.5)
plt.ylim(-2, 1)
plt.xlabel('log Energy (erg)')
plt.ylabel(r'log Cumulative Flare Rate (day$^{-1}$)')

print(therng(flat_samples2[:,0]))
print(therng(flat_samples2[:,1]))
print(therng(flat_samples2[:,2]))
Ndraws = 50
plt.figure(figsize=(9,7))
inds = np.random.randint(len(flat_samples2), size=Ndraws)
for ind in inds:
    sample = flat_samples2[ind]
    plt.plot([sample[2], max(ffd_x)], (sample[0]*np.array([sample[2], max(ffd_x)]) + sample[1]),
             "k", alpha=0.2)
    
plt.errorbar(ffd_x, ffd_ylog, xerr=ffd_xerr, yerr=ffd_yerrlog, 
             color='DodgerBlue', alpha=0.25, linestyle='none', lw=1, marker='.')

plt.plot(ffd_x[ffd_x >= np.median(flat_samples2[:,2])], 
         (np.median(flat_samples2[:,0]) * ffd_x[ffd_x >= np.median(flat_samples2[:,2])] + 
          np.median(flat_samples2[:,1])), 
         c='cyan', lw=3)

inds = np.random.randint(len(flat_samplesY2), size=Ndraws)
for ind in inds:
    sample = flat_samplesY2[ind]
    plt.plot([sample[2], max(ffd_xY)], (sample[0]*np.array([sample[2], max(ffd_xY)]) + sample[1]),
             "k", alpha=0.2)
    
plt.errorbar(ffd_xY, ffd_yYlog, xerr=ffd_xerrY, yerr=ffd_yerrYlog, 
             color='Firebrick', alpha=0.25, linestyle='none', lw=1, marker='.')

plt.plot(ffd_xY[ffd_xY >= np.median(flat_samplesY2[:,2])], 
         (np.median(flat_samplesY2[:,0]) * ffd_xY[ffd_xY >= np.median(flat_samplesY2[:,2])] + 
          np.median(flat_samplesY2[:,1])), 
         c='r', lw=3)

plt.xlim(30.1,33.1)
plt.ylim(-1.7, .9)
plt.xlabel('log Energy (erg)')
plt.ylabel(r'log Cumulative Flare Rate Per Day')

plt.savefig('YZCMI_vs_GJ1243_mcmc.png', dpi=150, bbox_inches='tight', pad_inches=0.25)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20019912690510838, 'fraction_numerical': 0.04174006280156238, 'mean_word_length': 3.3951531470885223, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 29, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '498181', 'license': 'MIT', 'n_tokens_mistral': 6247, 'n_tokens_neox': 5929, 'n_words': 890, 'path': 'yzcmi_vs_gj1243.ipynb', 'repo_name': 'jradavenport/YZCMi_TESS', 'size': '13057'}