如何从截断的正态分布中采样特定方差的随机数？

Question

我试图从给定数字的特定方差和界限的截尾正态分布抽样数字，例如，我需要平均值为0和单位方差的数字，但它们必须在一定范围内，例如[-2,2]

我无法弄清楚如何截断分布，同时保持方差。

import math 
import numpy as np 
import scipy.stats as stats 


truncation = 2 
lower, upper = -truncation, truncation 
mu, sigma = 0, 1 
num_samples = 1000 
if truncation: 
    n = stats.truncnorm((lower - mu)/sigma, (upper - mu)/sigma, loc=mu, scale=sigma) 
    samples = n.rvs(num_samples) 
    std_trunc = np.std(samples) 

    n = stats.norm(loc=mu, scale=sigma) 
    samples = n.rvs(num_samples) 
    std_simple = np.std(samples) 

print(std_trunc, std_simple, sep='\n') 

# outputs 
# 0.859167285015 # I need number close to 1 here 
# 1.01735583631 # like here, but here it's not truncated 

Answer 1

维基百科页面给出表达式为observed mean and variance，我们可以用它来反转，找出什么样的价值观，我们应该传递给truncnorm给我们我们想要的结果。

我们不会利用基于标准正常工作的任何简化，部分原因是一般原因，部分原因是我还没有吃过早餐，所以我不想做任何算术。可能你可以用一个简单的计算来代替整个最小化。

import numpy as np 
import scipy.stats as stats 
import scipy.optimize 

def truncated_mean_std(mu, sigma, lower, upper): 
    # N.B. lower/upper are the actual values, not Z-scaled 
    alpha = (lower - mu)/sigma 
    beta = (upper - mu)/sigma 
    d_pdf = (stats.norm.pdf(alpha) - stats.norm.pdf(beta)) 
    wd_pdf = (alpha * stats.norm.pdf(alpha) - beta * stats.norm.pdf(beta)) 
    d_cdf = stats.norm.cdf(beta) - stats.norm.cdf(alpha) 
    mu_trunc = mu + sigma * (d_pdf/d_cdf) 
    var_trunc = sigma**2 * (1 + wd_pdf/d_cdf - (d_pdf/d_cdf)**2) 
    std_trunc = var_trunc**0.5 
    return mu_trunc, std_trunc 

def trunc_samples(mu, sigma, lower, upper, num_samples=1000): 
    n = stats.truncnorm((lower - mu)/sigma, (upper - mu)/sigma, loc=mu, scale=sigma) 
    samples = n.rvs(num_samples) 
    return samples 

def corrector(mu, sigma, lower, upper): 
    target = np.array([mu, sigma]) 
    result = scipy.optimize.minimize(
     lambda x: ((target - truncated_mean_std(x[0], x[1], lower, upper))**2).sum(), 
     x0=[mu, sigma]) 
    return result.x 

这给了我：

In [79]: s = trunc_samples(mu=0, sigma=1, lower=-2, upper=2, num_samples=10**7) 

In [80]: s.mean(), s.std() 
Out[80]: (-9.8821067931585576e-05, 0.87951241887015619) 

In [81]: mu_to_use, sigma_to_use = corrector(0, 1, -2, 2) 

In [82]: mu_to_use, sigma_to_use 
Out[82]: (-7.4553057719882245e-09, 1.3778928137492246) 

In [83]: s = trunc_samples(mu=mu_to_use, sigma=sigma_to_use, lower=-2, upper=2, num_samples=10**7) 

In [84]: s.mean(), s.std() 
Out[84]: (0.0004091647648333381, 0.99991490259048865) 

In [85]: s.min(), s.max() 
Out[85]: (-1.9999995310631815, 1.9999997070340947)