File wellmet/voronoi.py changed (mode: 100644) (index fbb8d26..f642ead) |
... |
... |
class ContactVoronoi: |
397 |
397 |
|
|
398 |
398 |
č co máme zadavat jako aweight u np.cov? |
č co máme zadavat jako aweight u np.cov? |
399 |
399 |
č pdf bodíku? Nebo jejich vahy? |
č pdf bodíku? Nebo jejich vahy? |
|
400 |
|
č Vahy. |
|
401 |
|
č Rozepíšu zde dětailněji, protože to není na první pohled zřejmný |
|
402 |
|
č a myslím si, že i kontrintuitivní |
|
403 |
|
č 1. Uvaha číslo jedna: když se snažíme odhadnout rozptyl vzorků, |
|
404 |
|
č výbraných z normálního rozdělení, tak jich hustotami nepřevažíme |
|
405 |
|
č 2. Uvaha číslo dva: máme sadu IS a chceme odhadnout rozptyl skutečného |
|
406 |
|
rozdělení. Chce se tu použit váhy, že ano? |
|
407 |
|
č 3. Uvaha číslo tří: Jak bychom počítali těžiště a momenty setrvačnosti |
|
408 |
|
rovinného obrazce s odlíšnou tloušťkou t? (označme "m" jako mass) |
|
409 |
|
|
|
410 |
|
sum[y_i * A_i * t_i] Int y*t dA |
|
411 |
|
y_c = ---------------------- = ---------- |
|
412 |
|
sum(A_i*t_i) Int t dA |
|
413 |
|
|
|
414 |
|
sum[y_i * m_i] Int y dm |
|
415 |
|
y_c = ---------------- = ---------- |
|
416 |
|
sum(m_i) Int dm |
|
417 |
|
|
|
418 |
|
I = sum[(y_i - y_c)**2 * A_i * t_i] = Int (y-y_c)**2 * t dA |
|
419 |
|
|
|
420 |
|
I = sum[(y_i - y_c)**2 * m_i] = Int (y-y_c)**2 dm |
|
421 |
|
|
|
422 |
|
f_i je jako tloušťka t |
|
423 |
|
1/h_i je jako integrovační ploška dA |
|
424 |
|
w = f_i / h_i vahy jsou jako hmotnosti |
|
425 |
|
Monte Carlo - integrovaní dělením na kusy stejné hmotnosti |
|
426 |
|
IS - kusy mají odlíšnou hmotnost |
|
427 |
|
|
|
428 |
|
|
|
429 |
|
|
|
430 |
|
|
400 |
431 |
|
|
401 |
|
č jak počítat skore? čím máme dělit pravděpodobnostní obsah? |
|
402 |
|
č 1. sum((w_i - w_mean)**2) / n To je rozptyl w |
|
403 |
|
č 2. sum((w_i - w_mean)**2) / n**2 To je rozptyl odhadu w_mean |
|
404 |
|
č 3. sqrt(sum((w_i - w_mean)**2)) / n Směrodatná odhylka odhadu w_mean |
|
|
432 |
|
č Score počítáme jako pravděpodobnostní obsah buňky dělený |
|
433 |
|
č směrodatnou odchylkou odhadu (3.) Není to ale jediná možnost |
|
434 |
|
č 1. sum((w_i - w_mean)**2) / n To je rozptyl w |
|
435 |
|
č 2. sum((w_i - w_mean)**2) / n**2 To je rozptyl odhadu w_mean |
|
436 |
|
č 3. sqrt(sum((w_i - w_mean)**2)) / n Směrodatná odchylka odhadu w_mean |
405 |
437 |
""" |
""" |
406 |
438 |
def __init__(self, sample_box, hull, model_space='G', \ |
def __init__(self, sample_box, hull, model_space='G', \ |
407 |
439 |
p_norm=2, ns=1000, p_base=0.5, auto_update=True, \ |
p_norm=2, ns=1000, p_base=0.5, auto_update=True, \ |
|
... |
... |
class ContactVoronoi: |
666 |
698 |
#č p_couple není finální pravděpodobnost |
#č p_couple není finální pravděpodobnost |
667 |
699 |
#č bodíky ještě musejí projít konvexní obálkou |
#č bodíky ještě musejí projít konvexní obálkou |
668 |
700 |
nodes.p_couple = np.sum(w[mask]) |
nodes.p_couple = np.sum(w[mask]) |
|
701 |
|
nodes.score = nodes.p_couple / np.sqrt(np.var(nodes.w)/len(nodes)) |
669 |
702 |
#č a doporučíme (dočasně) k integraci |
#č a doporučíme (dočasně) k integraci |
670 |
703 |
self._recommend_to_integration(idx) |
self._recommend_to_integration(idx) |
671 |
704 |
|
|
|
... |
... |
class ContactVoronoi: |
956 |
989 |
|
|
957 |
990 |
|
|
958 |
991 |
|
|
|
992 |
|
# for simplex: d = nvar+2 |
|
993 |
|
# for cell: d = base_r**2 |
|
994 |
|
def IS_like(f_plan, sampling_space='G', weights=None, nis=1000, d=1, design=None): |
|
995 |
|
""" |
|
996 |
|
takes sample and returns sampling plan with the same cov (in case d=1) |
|
997 |
|
covariance matrix we'll divide by d |
|
998 |
|
""" |
|
999 |
|
|
|
1000 |
|
plan = getattr(f_plan, sampling_space) |
|
1001 |
|
|
|
1002 |
|
|
|
1003 |
|
S_bc = np.cov(plan, rowvar=False, bias=True, aweights=weights) |
|
1004 |
|
barycenter = np.average(plan, axis=0, weights=weights) |
|
1005 |
|
|
|
1006 |
|
|
|
1007 |
|
#č matika |
|
1008 |
|
w, v = np.linalg.eig(S_bc) |
|
1009 |
|
|
|
1010 |
|
# use IS sampling density with center equal to the simplex's barycenter |
|
1011 |
|
# set the minimum distance as the standard deviation of IS densisty |
|
1012 |
|
#č u stats.norm zadáváme směrodatnou odchylku, to je asi správné |
|
1013 |
|
sigmas = np.sqrt(w/d) |
|
1014 |
|
mean = 0 |
|
1015 |
|
h_plan_N, pdf_N = get_norm_plan(nis, f_plan.nvar, mean, sigmas, design) |
|
1016 |
|
|
|
1017 |
|
#ёӵ здесь уже так легко не отделаемся. Трансформовать кароно. |
|
1018 |
|
h_plan_bc = (v @ h_plan_N.T).T |
|
1019 |
|
h_plan_sing = h_plan_bc + barycenter |
|
1020 |
|
|
|
1021 |
|
|
|
1022 |
|
|
|
1023 |
|
#č sice to má nazev h_plan, ale nese rozdělení a hustoty v f-ku |
|
1024 |
|
# |
|
1025 |
|
#č jelikož výtvaříme nový vzorek z "čístých" souřadnic |
|
1026 |
|
#č máme zaručeno, že h_plan vždy bude f_modelem (podle současného kódu CandyBox) |
|
1027 |
|
h_plan = f_plan.new_sample(h_plan_sing, sampling_space) |
|
1028 |
|
## desired: w = h_plan.pdf(sampling_space) / pdf_N # snad je to správně |
|
1029 |
|
w = np.divide(h_plan.pdf(sampling_space), pdf_N, out=pdf_N); del(pdf_N) |
|
1030 |
|
|
|
1031 |
|
#č vahy máme |
|
1032 |
|
#č zabalme do boxu |
|
1033 |
|
#č zbytek už nejsou naši starosti |
|
1034 |
|
return CandyBox(h_plan, w=w) |
|
1035 |
|
|
959 |
1036 |
|
|
960 |
1037 |
def sample_alike(self, nodes): |
def sample_alike(self, nodes): |
961 |
1038 |
""" |
""" |