Graviness Blog

»»¿ô¡¦¿ô³Ø¡¦²Ê³Ø¡¦ÅÅǾ¡¦»¨µ­¡¦¥¢¥Û¤Î½ç¤ÎÌ©Å٤ǵ­»ö¤¬¹½À®¤µ¤ì¤Þ¤¹¡£
<< February 2019 | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 >> ¥Ö¥í¥°¥é¥ó¥­¥ó¥°¡¦¤Ë¤Û¤ó¥Ö¥í¥°Â¼¤Ø
 
RECOMMEND
¥Ó¥Ã¥°¥Ð¥ó±§ÃèÏÀ (¾å)
¥Ó¥Ã¥°¥Ð¥ó±§ÃèÏÀ (¾å) ¡ÊJUGEM¥ì¥Ó¥å¡¼ »¡Ë
¥µ¥¤¥â¥ó¡¦¥·¥ó, ÀÄÌÚ ·°
RECENT COMMENT
  • Ë­¿Ã½¨µÈ¤ÈÁ½Ï¤Íø¿·º¸±ÒÌ礫¤é³Ø¤Ö¿ôÎó¤ÎÏÂ
    ¢Í ͥǵ (07/12)
  • ¡Ú狼²ò¤¤¤Æ¡ÛÁ²²½¼° a_(n+1) = f(n) * a_n ^ g(n) + h(n) ¤Î°ìÈ̹à
    ¢Í ͥǵ (02/18)
  • ¡Ú狼²ò¤¤¤Æ¡ÛÁ²²½¼° a_(n+1) = f(n) * a_n ^ g(n) + h(n) ¤Î°ìÈ̹à
    ¢Í S.S.+ (02/16)
  • Ë­¿Ã½¨µÈ¤ÈÁ½Ï¤Íø¿·º¸±ÒÌ礫¤é³Ø¤Ö¿ôÎó¤ÎÏÂ
    ¢Í ºä°æ¾¼ (03/19)
  • d/dx(x¢¬¢¬n): ¹â¤µ¤¬Äê¿ô¤Î¥Æ¥È¥ì¡¼¥·¥ç¥ó¤ÎÈùʬ - ¿ô³ØÅªµ¢Ç¼Ë¡¤òÍѤ¤¤ëÊýË¡
    ¢Í ­£ (09/30)
  • Á´¤Æ¤Î»°³Ñ·Á¤ÏÆóÅùÊÕ»°³Ñ·Á
    ¢Í ͥǵ (09/28)
  • Á´¤Æ¤Î»°³Ñ·Á¤ÏÆóÅùÊÕ»°³Ñ·Á
    ¢Í µµ¥ì¥¹ (09/28)
  • Á´¤Æ¤Î»°³Ñ·Á¤ÏÆóÅùÊÕ»°³Ñ·Á
    ¢Í ͥǵ (09/24)
  • Á´¤Æ¤Î»°³Ñ·Á¤ÏÆóÅùÊÕ»°³Ñ·Á
    ¢Í µµ¥ì¥¹ (09/23)
  • ¡Ṳ́²ò·è¡Û¿·¤·¤¤±é»»»Ò¤òÁϤë
    ¢Í $_ (09/10)
RECENT TRACKBACK
MOBILE
qrcode
PROFILE
 
ÅÀ¤Î½Å¤Ê¤ê¤ò¹Íθ¤·¤¿ºÇŬ¤ÊÅÀ¤Î¥µ¥¤¥º¤ÎƳ½Ð

¥Ç¡¼¥¿¤ò»¶ÉÛ¿ÞÅù¤Çɽ¼¨¤¹¤ë¤È¤­¡¢¥Ç¡¼¥¿¿ô¤¬Â¿¤¯¤Æ²èÌ̤¬ÅÀ¤ÇÅɤêÄÙ¤ì¤Æ¤·¤Þ¤¦·Ð¸³¤Ï¤Ê¤¤¤Ç¤·¤ç¤¦¤«¡£

²¼¿Þ¤Ï10000¸Ä¤Î¥Ç¡¼¥¿¤ò»¶ÉÛ¿Þ¤Çɽ¼¨¤·¤¿Îã¤Ç¤¹¡£²¼¿Þ¡Êº¸¡Ë¤ÏÃæ±û¤ÎÂÓ¤ÎÉôʬ¤¬ÅÀ¤ÇÅɤêÄÙ¤ì¤Æ¤¤¤Þ¤¹¡£ÉÁ²èÎΰè¤òÂ礭¤¯¤¹¤ë¤«¡¢²¼¿Þ¡Ê±¦¡Ë¤Î¤è¤¦¤ËÅÀ¤Î¥µ¥¤¥º¤ò¾®¤µ¤¯Ä´À᤹¤ë¤³¤È¤ÇÅɤêÄ٤줬¤Ê¤¯¤Ê¤ê¡¢Àº³Î¤Ê¥Ç¡¼¥¿Ê¬ÀÏ¡¢¤½¤·¤ÆÈþ¤·¤¯É½¼¨¤Ç¤­¤Þ¤¹¡£

¤·¤«¤·¡¢¸ÇÄêŪ¤Ê¥Ç¡¼¥¿¿ô¤ò°·¤¦¾ì¹ç¤ÏÌäÂê¤È¤Ê¤ê¤Þ¤»¤ó¤¬¡¢¥Ç¡¼¥¿¿ô¤¬ÊѤï¤ë¤è¤¦¤Ê¥·¥Á¥å¥¨¡¼¥·¥ç¥ó¤Î¾ì¹ç¤Ï¤½¤ÎÅÔÅÙÄ´À᤹¤ëɬÍפ¬¤¢¤êÌÌÅݤǤ¹¡£

Ëܵ­»ö¤Ç¤Ï¡¢¥Ç¡¼¥¿¿ô¤Ë±þ¤¸¤ÆºÇŬ¤ÊÅÀ¤Î¥µ¥¤¥º¤ò¼«Æ°·×»»¤¹¤ë¼°¤È¤½¤ÎƳ½Ð¤ò¼¨¤·¤Þ¤¹¡£·ëÏÀ¤ò¼¨¤¹¤È¡¢x¸Ä¤ÎÅÀ¤òÉÁ²è¤¹¤ë¤È¤­¤Î¸Ä¡¹¤ÎÅÀ¤Î¥µ¥¤¥ºs(x)¤Ï¼¡¼°¤ÇÆÀ¤é¤ì¤Þ¤¹¡£

s(x) = S1 * S¡ç / (S1 * x + S¡ç)

S1: ÅÀ¤ò1¸ÄÉÁ²è¤¹¤ë¤È¤­¤Î¡Êx=1¤Î¤È¤­¤Î¡Ë¡¢1¸Ä¤ÎÅÀ¤Ë¤è¤ëÉÁ²èÌÌÀÑ¡£
S¡ç: ÁíÉÁ²èÌÌÀѤξå¸Â¡£e.g. ²èÌ̤ÎÈ¾Ê¬ÄøÅÙ¤¬Ëä¤Þ¤ë¤Î¤Ç¤¢¤ì¤Ð¤½¤ÎÌÌÀÑ¡£

Î㤨¤Ð¡¢S1=36px¡¢S¡ç=3600px¤Î¾ì¹ç¡¢s(x) = 3600 / (x + 100) ¤È¤Ê¤ê¤Þ¤¹¡£100ÅÀÉÁ²è¤¹¤ë¤È¤­¸Ä¡¹¤ÎÅÀ¤Î¥µ¥¤¥º¤Ïs(100)=18px¤ÇÉÁ²è¤µ¤ì¡¢Æ±Íͤˤ·¤Æ1000ÅÀÉÁ²è¤¹¤ë¤È¤­s(1000)=3.3px¡¢10000ÅÀÉÁ²è¤¹¤ë¤È¤­s(10000)=0.4px¤È¤Ê¤ê¤Þ¤¹¡£

Python¥×¥í¥°¥é¥à¤ÎÎã¤ò¼¨¤·¤Þ¤¹¡£¥Ç¡¼¥¿¿ôlen(data)¤Ë±þ¤¸¤ÆÅÀ¤ÎÉÁ²è¥µ¥¤¥ºs¤ò¼«Æ°·×»»¤·¤Æ¤¤¤Þ¤¹¡£

import numpy as np
import matplotlib
import matplotlib.pyplot as plt

data = np.random.normal(size=1000) # normal dist.

S_1, S_m = 36, 3600
s = S_1 * S_m / (S_1 * len(data) + S_m) # see

plt.title(f"count={len(data)}")
plt.ylim((-4.0, 4.0))
plt.scatter(x=np.arange(len(data)), y=data, s=s)
plt.savefig(f"sample.{len(data)}.png")
plt.close()

°Ê¹ß¤Ç·ëÏÀ¤Ë»ê¤ë¤Þ¤Ç¤ÎƳ½Ð¤ò¼¨¤·¤Þ¤¹¡£

JUGEM¥Æ¡¼¥Þ¡§¥³¥ó¥Ô¥å¡¼¥¿

³¤­¤òÆÉ¤à ¡ä¡ä
Smirnov-Grubbs¸¡Äê¤òÍѤ¤¤ë³°¤ìÃͽüµî¤ÎPython¼ÂÁõ

¥¹¥ß¥ë¥Î¥Õ¡Ý¥°¥é¥Ö¥¹¸¡Äê¤òÍѤ¤¤ë³°¤ìÃͽüµî¤ÎPython¡ÊNumPy/SciPy»ÈÍѡˤμÂÁõ¤Ç¤¹¡£

 

²¼µ­smirnov_grubbs´Ø¿ô¤Ë¥µ¥ó¥×¥ë¥Ç¡¼¥¿¤ÈÍ­°Õ¿å½à¦Á¡Êξ¦100¦Á¡ó¡Ë¤òÆþÎϤ¹¤ë¤È¡¢³°¤ìÃͤò½üµî¤·¤¿¥µ¥ó¥×¥ë¥Ç¡¼¥¿¤È¡¢Ãê½Ð¤·¤¿³°¤ìÃͥǡ¼¥¿¤ò½ÐÎϤ·¤Þ¤¹¡£

 

¥µ¥ó¥×¥ë¥Ç¡¼¥¿¤Î¤¦¤ÁÊ¿¶ÑÃͤ«¤é¤ÎÊк¹¤¬ºÇ¤âÂ礭¤¤¥Ç¡¼¥¿¤ò¡¢Í­°Õ¿å½à¡Êξ¦¡Ë¤Ç¥¹¥ß¥ë¥Î¥Õ¡Ý¥°¥é¥Ö¥¹¸¡Äꤷ¡¢µ¢Ìµ²¾Àâ¤ò´þµÑ¤·¤¿¾ì¹ç¤ËÅö³º¥Ç¡¼¥¿¤ò³°¤ìÃͤȤ·¤Þ¤¹¡£Åö³º¥Ç¡¼¥¿¤ò½ü¤¤¤¿¥µ¥ó¥×¥ë¥Ç¡¼¥¿¤ËÂФ·¡¢³°¤ìÃͤ¬¤Ê¤¯¤Ê¤ë¤Þ¤Ç¤³¤ÎÁàºî¤ò·«¤êÊÖ¤·¤Þ¤¹¡£

 

ref. ¥¹¥ß¥ë¥Î¥Õ¡Ý¥°¥é¥Ö¥¹¸¡Äê¡÷wikipedia

 

# coding: utf-8
#
# smirnov_grubbs(data, alpha)
#   Generate outlier-removed data with Smirnov-Grubbs.
#
#   Parameters:
#     data: array_like
#       Numeric data values.
#     alpha: float
#       Significance level. (two-sided)
#
#   Returns:
#     outlier-removed-data
#     outlier-data
#
#   e.g.
#     data = [1, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 1000]
#     alpha = 0.01
#     smirnov_grubbs(data, alpha)
#     ==> ([100, 101, 102, 103, 104, 105, 106, 107, 108, 109], [1000, 1])

import numpy as np
import scipy.stats as stats

def smirnov_grubbs(data, alpha):
    x, o = list(data), []
    while True:
        n = len(x)
        t = stats.t.isf(q=(alpha / n) / 2, df=n - 2)
        tau = (n - 1) * t / np.sqrt(n * (n - 2) + n * t * t)
        i_min, i_max = np.argmin(x), np.argmax(x)
        myu, std = np.mean(x), np.std(x, ddof=1)
        i_far = i_max if np.abs(x[i_max] - myu) > np.abs(x[i_min] - myu) else i_min
        tau_far = np.abs((x[i_far] - myu) / std)
        if tau_far < tau: break
        o.append(x.pop(i_far))
    return (np.array(x), np.array(o))

 

JUGEM¥Æ¡¼¥Þ¡§¥³¥ó¥Ô¥å¡¼¥¿

³¤­¤òÆÉ¤à ¡ä¡ä
Brunner-Munzel¸¡Äê¤ÎPython¼ÂÁõ

Brunner-Munzel¸¡Äê¤ÎPython¡ÊNumPy/SciPy»ÈÍѡˤμÂÁõ¤Ç¤¹¡£

 

²¼µ­bmtest´Ø¿ô¤Ë2·²¤Î¥Ç¡¼¥¿¤òÆþÎϤ¹¤ë¤È¡¢¸¡ÄêÅý·×ÎÌW¡¢pÃÍ¡Êξ¦¡Ë¡¢¼«Í³ÅÙ¡¢¤ª¤è¤Óµ¢Ìµ²¾Àâ¤Ç¼¨¤µ¤ì¤ë³ÎΨ¤Î¿äÄêÃͤò½ÐÎϤ·¤Þ¤¹¡£

ÊÒ¦¸¡Äê¤Ë¤¹¤ë¾ì¹ç¡Ê¤É¤Á¤é¤Ë¤º¤ì¤Æ¤¤¤ë¤«¤¬½ÅÍפʾì¹ç¡Ë¡¢½ÐÎϤµ¤ì¤¿pÃͤò0.5Çܤ·¤Æ¡¢W¤ÎÃͤ¬ÀµÉé¤ÇȽÃǤ·¤Þ¤¹¡£

Ãí°Õ: ¾¯¤·¤â¤«¤Ö¤ê¤â¤·¤Ê¤¤2·²¤òÆþÎϤ¹¤ë¤È¡¢Ê¬»¶¤¬0¤È¤Ê¤Ã¤Æ¡¢¸¡ÄêÅý·×ÎÌW¤ä¼«Í³Å٤η׻»Ãæ¤Ë Runtimewarning; divide by zero ¤¬½Ð¤Þ¤¹¡£

 

ref. Brunner-Munzel¸¡Äê¡÷±ü¼ À²É§¤µ¤ó

 

# coding: utf-8
#
# bmtest(x1, x2)
#   Calculate Brunner-Munzel-test scores.
#
#   Parameters:
#     x1, x2: array_like
#       Numeric data values from sample 1, 2.
#
#   Returns:
#     w:
#       Calculated test statistic.
#     p_value:
#       Two-tailed p-value of test.
#     dof:
#       Degree of freedom.
#     p:
#       "P(x1 < x2) + 0.5 P(x1 = x2)" estimates.
#
#   References:
#     * https://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.html

import numpy as np
import scipy.stats as stats

def bmtest(x1, x2):
    n1, n2 = len(x1), len(x2)
    R = stats.rankdata(list(x1) + list(x2))
    R1, R2 = R[:n1], R[n1:]
    r1_mean, r2_mean = np.mean(R1), np.mean(R2)
    Ri1, Ri2 = stats.rankdata(x1), stats.rankdata(x2)
    var1 = np.var([r - ri for r, ri in zip(R1, Ri1)], ddof=1)
    var2 = np.var([r - ri for r, ri in zip(R2, Ri2)], ddof=1)
    w = ((n1 * n2) * (r2_mean - r1_mean)) / ((n1 + n2) * np.sqrt(n1 * var1 + n2 * var2))
    dof = (n1 * var1 + n2 * var2) ** 2 / ((n1 * var1) ** 2 / (n1 - 1) + (n2 * var2) ** 2 / (n2 - 1))
    c = stats.t.cdf(abs(w), dof) if not np.isinf(w) else 0.0
    p_value = min(c, 1.0 - c) * 2.0
    p = (r2_mean - r1_mean) / (n1 + n2) + 0.5
    return (w, p_value, dof, p)

 

JUGEM¥Æ¡¼¥Þ¡§¥³¥ó¥Ô¥å¡¼¥¿

³¤­¤òÆÉ¤à ¡ä¡ä
d/dx(x¢¬¢¬n): ¹â¤µ¤¬Äê¿ô¤Î¥Æ¥È¥ì¡¼¥·¥ç¥ó¤ÎÈùʬ

²áµî¤Îµ­»ö¡Öd/dx(x↑↑n): ¹â¤µ¤¬Äê¿ô¤Î¥Æ¥È¥ì¡¼¥·¥ç¥ó¤ÎÈùʬ - ¿ô³ØÅªµ¢Ç¼Ë¡¤òÍѤ¤¤ëÊýË¡¡×¤Ë¤Ä¤¤¤Æ¡¢¹â¤µ¤¬Äê¿ô¤Î¥Æ¥È¥ì¡¼¥·¥ç¥ó¤ÎƳ´Ø¿ô¤ò¿ô³ØÅªµ¢Ç¼Ë¡¤òÍѤ¤¤º¤Ëµá¤á¤Þ¤¹¡£É½µ­¤È¤·¤Æ¡¢nx := x↑↑n ¤òÍѤ¤¤Þ¤¹¡£

 

---

³¤­¤òÆÉ¤à ¡ä¡ä
¡Ú狼²ò¤¤¤Æ¡ÛÁ²²½¼° a_(n+1) = f(n) * a_n ^ g(n) + h(n) ¤Î°ìÈ̹à

²áµî¤ÎÆó¤Ä¤Îµ­»ö¤ÎÁ²²½¼°¤òÊñ´Þ¤¹¤ë a_(n+1) = f(n) * a_n ^ g(n) + h(n) ¤Î·Á¼°¤Î°ìÈ̹à¤òµá¤á¤¿¤¤¤Ç¤¹¡£

* Á²²½¼° a_(n+1) = f(n) * a_n + g(n) ¤Î°ìÈ̹à

* Á²²½¼° a_(n+1) = f(n) * a_n ^ g(n) ¤Î°ìÈ̹à

 

Ä̶лþ¤Ê¤É¤Ç¿ô¥ö·î¹Í¤¨¤Æ¤Þ¤¹¤¬¡¢¤³¤ì¤¬¤É¤¦¤ä¤Ã¤Æ¤â²ò¤±¤Ê¤¤¡£¡£¡£¶²¤é¤¯½éÅù´Ø¿ô¤ÎÈϰϤǤϲò¤±¤Ê¤¤¤Î¤Ç¤Ï¤Ê¤¤¤«¤È»×¤¤¤Þ¤¹¡£

 

´ØÏ¢¤·¤Æ¡¢ÎÙÀÜ»°¹à´ÖÁ²²½¼° a_(n+2) = f(n) * a_(n+1) + g(n) * a_n ¤Ï¡¢Î¾ÊÕ¤ò a_(n+1) ¤Ç³ä¤ë¤È¡¢a_(n+2)/a_(n+1) = f(n) + g(n) * a_n/a_(n+1) ¤È¤Ê¤ê¡¢b_n = a_(n+1)/a_n ¤È¤ª¤±¤Ð¡¢b_(n+1) = g(n) b_n ^ (-1) + f(n) ¤È¤Ê¤ë¤¿¤á¡¢Ëܵ­»ö¤ÎÁ²²½¼°¤¬²ò¤±¤ì¤ÐÎÙÀÜ»°¹à´ÖÁ²²½¼°¤â²ò¤¯¤³¤È¤¬¤Ç¤­¤Þ¤¹¡£

 

JUGEM¥Æ¡¼¥Þ¡§³ØÌ䡦³Ø¹»