如何從截斷的正態分佈中採樣特定方差的隨機數？

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)