sympyで信頼区間

Python Sympy, quantiles of a distribution – StackOverflow

“distribution.confidence”メソッドが無くなっているなら流用すれば，と思ったけど，
sympy.mpmathが無くなって（それ自体は外部モジュールを使えばいいけど），
‘RandomSymbol’ objectのインプリメントもだいぶ変わっているっぽい．
各種インスタンスメソッドが無い．ちまちまと計算するしか無さそう．
（Sympyへの理解はまだまだ乏しいのでパッと方法が思い付かない）
　
　

計算的にやるなら，SciPyが楽で，　
　

from scipy import stats

print(stats.norm.interval(0.67, loc=0, scale=1))
print(stats.norm.interval(0.90, loc=0, scale=1))
print(stats.norm.interval(0.95, loc=0, scale=1))
print(stats.norm.interval(0.99, loc=0, scale=1))

(-0.9741138770593093, 0.97411387705930919)
(-1.6448536269514729, 1.6448536269514722)
(-1.959963984540054, 1.959963984540054)
(-2.5758293035489004, 2.5758293035489004)
　
　

ベイジアン信頼区間を求めるなら，
　

from scipy.stats import bayes_mvs, mvsdist

r = stats.norm.rvs(size=10000)
mean, var, std = mvsdist(r)

print(bayes_mvs(r, alpha=0.67))
print(bayes_mvs(r, alpha=0.90))
print(bayes_mvs(r, alpha=0.95))
print(bayes_mvs(r, alpha=0.99))

　

(Mean(statistic=-0.0046724190670562431, minmax=(-0.014453741852972777, 0.0051089037188602887)), Variance(statistic=1.0082673904053563, minmax=(0.99437744785478244, 1.0221573329559301)), Std_dev(statistic=1.0041251866203518, minmax=(0.99720874694945572, 1.0110416262912478)))

(Mean(statistic=-0.0046724190670562431, minmax=(-0.021188808618314344, 0.011843970484201851)), Variance(statistic=1.0082673904053563, minmax=(0.98481333342032484, 1.0317214473903877)), Std_dev(statistic=1.0041251866203518, minmax=(0.99244633556793849, 1.0158040376727651)))

(Mean(statistic=-0.0046724190670562431, minmax=(-0.024352911084510744, 0.015008072950398257)), Variance(statistic=1.0082673904053563, minmax=(0.98032015775875925, 1.0362146230519533)), Std_dev(statistic=1.0041251866203518, minmax=(0.99020897725772195, 1.0180413959829815)))

(Mean(statistic=-0.0046724190670562431, minmax=(-0.030536969868338346, 0.02119213173422586)), Variance(statistic=1.0082673904053563, minmax=(0.97153850080736881, 1.0449962800033437)), Std_dev(statistic=1.0041251866203518, minmax=(0.9858361873564212, 1.0224141858842823)))
　

mean.interval(0.67), mean.interval(0.90), mean.interval(0.95), mean.interval(0.99)

((-0.014453741852972777, 0.0051089037188602887),
(-0.021188808618314344, 0.011843970484201851),
(-0.024352911084510744, 0.015008072950398257),
(-0.030536969868338346, 0.02119213173422586))
　
　

任意の分布に従う乱数が欲しい場合は，

from scipy.stats import rv_continuous

class custom_distributions(rv_continuous):
    def _pdf(self, x):
        return np.exp(-x**2) / np.sqrt(np.pi) * 2


custom = custom_distributions(name='custom', a=0)
r = custom.rvs(size=10000)

mean, var, std = mvsdist(r)
print(bayes_mvs(r, alpha=0.67))
print(bayes_mvs(r, alpha=0.90))
print(bayes_mvs(r, alpha=0.95))
print(bayes_mvs(r, alpha=0.99))

(Mean(statistic=0.56167815956488576, minmax=(0.55753814777696609, 0.56581817135280543)), Variance(statistic=0.18062742671159079, minmax=(0.17813909414826815, 0.18311575927491344)), Std_dev(statistic=0.42500285494522361, minmax=(0.42207542453579339, 0.42793028535465383)))
(Mean(statistic=0.56167815956488576, minmax=(0.55468748469067197, 0.56866883443909955)), Variance(statistic=0.18062742671159079, minmax=(0.17642571791939227, 0.18482913550378932)), Std_dev(statistic=0.42500285494522361, minmax=(0.42005970133659659, 0.42994600855385062)))
(Mean(statistic=0.56167815956488576, minmax=(0.55334825667469234, 0.57000806245507918)), Variance(statistic=0.18062742671159079, minmax=(0.17562078188235009, 0.1856340715408315)), Std_dev(statistic=0.42500285494522361, minmax=(0.41911272412494244, 0.43089298576550478)))
(Mean(statistic=0.56167815956488576, minmax=(0.55073081148628722, 0.5726255076434843)), Variance(statistic=0.18062742671159079, minmax=(0.1740475800586197, 0.18720727336456189)), Std_dev(statistic=0.42500285494522361, minmax=(0.41726191088283709, 0.43274379900761012)))

mean.interval(0.67), mean.interval(0.90), mean.interval(0.95), mean.interval(0.99)

((0.55753814777696609, 0.56581817135280543),
(0.55468748469067197, 0.56866883443909955),
(0.55334825667469234, 0.57000806245507918),
(0.55073081148628722, 0.5726255076434843))

sns.boxplot(r)

sns.distplot(r)

　
　
　

class gaussian_gen(rv_continuous):
    def _pdf(self, x):
        return np.exp(-x**2 / 2.) / np.sqrt(2.0 * np.pi)


gaussian = gaussian_gen(name='gaussian')
r = gaussian.rvs(size=100)
sns.distplot(r)

class custom_distributions2(rv_continuous):
    def _pdf(self, x):
        return np.exp(-x**2) / np.sqrt(np.pi) * 2


custom2 = custom_distributions2(name='custom')
r = custom2.rvs(size=100)
sns.distplot(r)

質問のパラメータでやると結果が変．

f = exp(-x**2)/sqrt(pi)*2
plot(f)

g = exp(-x**2 / 2.) / sqrt(2.0 * pi)
plot(g)