File wellmet/IS_stat.py changed (mode: 100644) (index 789008d..5e55e0e) |
... |
... |
def IS_like(f_plan, sampling_space='G', weights=None, nis=1000, d=1, design=None |
366 |
366 |
|
|
367 |
367 |
|
|
368 |
368 |
|
|
369 |
|
|
|
370 |
369 |
# for simplex: d = nvar+2 |
# for simplex: d = nvar+2 |
371 |
370 |
# for cell: it's not so easy |
# for cell: it's not so easy |
372 |
371 |
def sample_alike(plan, weights=None, nis=1000, d=1, design=None): |
def sample_alike(plan, weights=None, nis=1000, d=1, design=None): |
373 |
372 |
""" |
""" |
374 |
373 |
takes sample and returns sampling plan with the same cov (in case d=1) |
takes sample and returns sampling plan with the same cov (in case d=1) |
375 |
374 |
covariance matrix we'll divide by d |
covariance matrix we'll divide by d |
376 |
|
""" |
|
377 |
375 |
|
|
|
376 |
|
Returns: sampling_plan, h_pdfs |
|
377 |
|
""" |
378 |
378 |
nsim, ndim = plan.shape |
nsim, ndim = plan.shape |
379 |
379 |
|
|
|
380 |
|
#ё зрада отменяется! |
|
381 |
|
#č Tahle funkce se pokusí aspoň něco nejak navzorkovat, |
|
382 |
|
#č ale musí být aspoň-alespoň dva vzorky. |
|
383 |
|
#č nechvám zodpovědnost za to na volajícím kódě! |
|
384 |
|
assert nsim > 1 |
|
385 |
|
|
|
386 |
|
S_bc = np.cov(plan, rowvar=False, bias=True, aweights=weights) |
|
387 |
|
barycenter = np.average(plan, axis=0, weights=weights) |
|
388 |
|
|
|
389 |
|
|
|
390 |
|
|
380 |
391 |
if nsim < ndim + 1: |
if nsim < ndim + 1: |
381 |
|
return False, None, None |
|
|
392 |
|
S2 = np.diag(S_bc) |
|
393 |
|
if not np.all(S2 > 0): |
|
394 |
|
std = np.sqrt(np.sum(S2) / ndim / d) |
|
395 |
|
else: |
|
396 |
|
std = np.sqrt(S2 / d) |
|
397 |
|
return get_norm_plan(nis, ndim, mean=barycenter, \ |
|
398 |
|
std=std, design=design) |
382 |
399 |
|
|
383 |
|
S_bc = np.cov(plan, rowvar=False, bias=True, aweights=weights) |
|
|
400 |
|
|
384 |
401 |
#č matika |
#č matika |
385 |
402 |
w, v = np.linalg.eig(S_bc) |
w, v = np.linalg.eig(S_bc) |
386 |
403 |
|
|
387 |
404 |
if not np.all(w > 0): |
if not np.all(w > 0): |
388 |
|
return False, None, None |
|
389 |
|
|
|
390 |
|
barycenter = np.average(plan, axis=0, weights=weights) |
|
391 |
|
|
|
|
405 |
|
S2 = np.diag(S_bc) |
|
406 |
|
if not np.all(S2 > 0): |
|
407 |
|
std = np.sqrt(np.sum(w) / ndim / d) |
|
408 |
|
else: |
|
409 |
|
std = np.sqrt(S2 / d) |
|
410 |
|
return get_norm_plan(nis, ndim, mean=barycenter, \ |
|
411 |
|
std=std, design=design) |
392 |
412 |
|
|
393 |
413 |
# use IS sampling density with center equal to the simplex's barycenter |
# use IS sampling density with center equal to the simplex's barycenter |
394 |
414 |
# set the minimum distance as the standard deviation of IS densisty |
# set the minimum distance as the standard deviation of IS densisty |
|
... |
... |
def sample_alike(plan, weights=None, nis=1000, d=1, design=None): |
401 |
421 |
h_plan_bc = (v @ h_plan_N.T).T |
h_plan_bc = (v @ h_plan_N.T).T |
402 |
422 |
h_plan_sing = h_plan_bc + barycenter |
h_plan_sing = h_plan_bc + barycenter |
403 |
423 |
|
|
404 |
|
return True, h_plan_sing, pdf_N |
|
405 |
|
|
|
406 |
|
|
|
|
424 |
|
return h_plan_sing, pdf_N |
407 |
425 |
|
|
408 |
426 |
|
|
409 |
427 |
|
|