scipy 给定随机变量列表和期望值,如何在python中生成概率分布

w41d8nur  于 2022-11-10  发布在  Python
关注(0)|答案(1)|浏览(169)

我有一个随机值x_1,x_2,... x_n的列表和一个期望值X。我想写一个函数,它取这两个值,并随机生成集合{1,2...n}上定义的许多离散概率分布中的一个,满足上面提到的约束。
将问题改写为生成一个长度为n的向量p,使得
0〈= p_i〈=1
| 有害物|= 1个
p.x = X
使用我找到的here信息,我破解了下面的解决方案,但它不起作用,因为概率不是b/w 0和1

def normal_random_distribution(data_array, data_len, X_array, expd_vl):
    e_mat = np.concatenate(([X_array], [np.ones(data_len)]), axis = 0)
    f_mat = np.array([expd_vl, 1])
    v_vec = np.linalg.lstsq(e_mat, f_mat,rcond=None)[0]
    e_null = sla.null_space(e_mat)
    lamda_lower = np.linalg.pinv(e_null) @ (-1 * v_vec)
    lamda_upper = np.linalg.pinv(e_null) @ ( 1 - v_vec)
    lamda = lamda_lower + random.random()*(lamda_upper - lamda_lower)
    sol = v_vec + e_null @ lamda
    return sol

如何生成p?sol有l1范数1,且|溶胶 * x| = X。我希望插值B/w lamda_lower和lamda_upper会使它b/w 0和1,但它没有。我想写函数,使它在每次调用时随机生成许多可能的解之一。我意识到一个明显的解是找到x_i〈X和x_j〉X,并插值b/w 2,但我希望函数(理论上)能够生成所有可能的分布。

0md85ypi

0md85ypi1#

在尝试创建一个递归解决方案几个小时后,我放弃了,自己编写了一个迭代的蛮力方法。

@numba.jit(nogil = True)  
def normal_random_distribution(data_array, X_array, expd_vl):
    # print("here")
    EPSILON = 0.001
    data_len = data_array.size
    data_array[:] = 1 / data_len
    current_E = np.sum(data_array * X_array)
    middle = np.argmax(X_array > expd_vl) - 1
    while abs(current_E - expd_vl) > EPSILON:
        # print(abs(current_E - expd_vl) / expd_vl)
        left_i = random.randint( 0, middle )
        rght_i = random.randint( middle + 1, data_len - 1 )
        if current_E < expd_vl:
            up_lim = min(data_array[left_i] , ( expd_vl - current_E ) / ( X_array[rght_i] - X_array[left_i] ) )
            mu = up_lim / 2
            sigma = up_lim / 6
            x = np.clip(np.random.normal( mu, sigma, 1 ), 0, up_lim)[0]
            data_array[left_i] -= x
            data_array[rght_i] += x
            current_E += x * ( X_array[rght_i] - X_array[left_i] )
        else:
            up_lim = min(data_array[rght_i] , ( current_E - expd_vl ) / ( X_array[rght_i] - X_array[left_i] ) ) 
            mu = up_lim / 2
            sigma = up_lim / 6
            x = np.clip(np.random.normal( mu, sigma, 1 ), 0, up_lim)[0]
            data_array[left_i] += x
            data_array[rght_i] -= x
            current_E -= x * ( X_array[rght_i] - X_array[left_i] )

相关问题