human
I am working on the file "Foil Open Area/Open Area-pymc3.ipynb".
The first task is:
Looks to go all the way to all sides in X-Y.
Work to total open area calculations
Now we can model the total open area of the foil given the noise estimate per pixel and the pixels that are a part of the witness sample and the total area.
We model the observed background as Poisson with center at the real background:
$obsnbkg \sim Pois(nbkg)$
We model the observed witness sample, $obswit$, as Poisson with center of background per pixel times number of pixels in peak plus the number of real counts:
$obswit \sim Pois(nbkg/C + witc)$, $C = \frac{A_w}{A_t}$
This then leaves the number of counts in open areas of the system (excluding witness) as a Poisson with center of background per pixel times number of pixels in the system (less witness) plus the real number of counts.
$obsopen \sim Pois(nbkg/D + realc)$, $D=\frac{A_t - A_w}{A_t}$
Then then the open area is given by the ratio number of counts, $realc$, over an unknown area, $A_o$, as related to witness counts, $witc$, to the witness area, $A_w$, which is assumed perfect as as 0.6mm hole.
$\frac{A_o}{realc}=\frac{A_w}{witc} => A_o = \frac{A_w}{witc}realc $
Can you write Python code for it?

What is a the Python code you requested?
assistant
#Here is the Python code you requested.

Aw = np.pi*(0.2/2)**2 # mm**2
Af = 182.75 # mm**2  this is the area of the foil
W_F_ratio = Aw/Af

print(Aw, Af, W_F_ratio)

C = wit_pixels/n_pixels
D = (n_pixels-wit_pixels)/n_pixels
print('C', C, 'D', D)


with mc3.Model() as model:
    nbkg = mc3.Uniform('nbkg', 1, n_noise*5) # just 1 to some large number
    obsnbkg = mc3.Poisson('obsnbkg', nbkg, observed=n_noise)

    witc = mc3.Uniform('witc', 0, witness_counts*5) # just 0 to some large number
    obswit = mc3.Poisson('obswit', nbkg*C + witc, observed=witness_counts)

    realc = mc3.Uniform('realc', 0, total_counts*5) # just 0 to some large number
    obsopen = mc3.Poisson('obsopen', nbkg*D + realc, observed=total_counts-witness_counts)

#     @mc.deterministic(plot=True)
#     def open_area(realc=realc, witc=witc):
#         return realc*Aw/witc/Af
    
    
    open_area = mc3.Deterministic('open_area', realc*Aw/witc/Af)
    start = mc3.find_MAP()
    trace = mc3.sample(20000, start=start)
    

mc3.traceplot(trace);
mc3.summary(trace)


print("Foil 1 \n")

witc_mean = np.mean(trace['witc'])
witc_std = np.std(trace['witc'])

print("Found witness counts of {0} turn into {1} +/- {2}  ({3:.2f}%)\n".format(witness_counts, witc_mean, witc_std, witc_std/witc_mean*100))

realc_mean = np.mean(trace['realc'])
realc_std = np.std(trace['realc'])

print("Found non-witness counts of {0} turn into {1} +/- {2}  ({3:.2f}%)\n".format(total_counts-witness_counts, realc_mean, realc_std, realc_std/realc_mean*100))

nbkg_mean = np.mean(trace['nbkg'])
nbkg_std = np.std(trace['nbkg'])

print("Found noise counts of {0} turn into {1} +/- {2}  ({3:.2f}%)\n".format(0, nbkg_mean, nbkg_std, nbkg_std/nbkg_mean*100))

OA_median = np.median(trace['open_area'])
OA_mean = np.mean(trace['open_area'])
OA_std = np.std(trace['open_area'])
print("The open area fraction is {0:.6f} +/- {1:.6f}   ({2:.2f}%) at the 1 stddev level from 1 measurement\n".format(OA_mean, OA_std,OA_std/OA_mean*100 ))
print("Phil got {0} for 1 measurement\n".format(0.00139))
print("The ratio Brian/Phil is:...