更新
通常,当人们问这样的问题时,有一个隐含的假设/要求——所有随机数都应该以相同的方式分布。这意味着,如果我从抽样数组中为索引为0的项绘制边际概率密度函数(pdf),我将得到与我为数组中最后一项绘制边际概率密度函数相同的分布。人们通常对随机数组进行抽样,然后将其传递给其他例程来做一些有趣的事情。如果项目0的边际PDF与上一个索引项目的边际PDF不同,则仅恢复数组将产生与使用此类随机值的代码截然不同的结果。
在这里,我使用我的采样程序,绘制了原始条件下([-50…50]和=300)项目0和最后一个项目(29)的随机数分布。看起来很像,不是吗?
好的,这是您的采样程序的图片,相同的原始条件([-50…50]和=300),相同的样本数
更新II
用户应该检查采样例程的返回值,如果(并且仅当)返回值为真,则接受并使用采样数组。这是接受/拒绝方法。如图所示,下面是用于柱状图样本的代码:
int[]hh=new int[100];//已分配柱状图
A
B
MathDotNet
using System;
using MathNet.Numerics.Distributions;
using MathNet.Numerics.Random;
class Program
{
static void SampleDirichlet(double alpha, double[] rn)
{
if (rn == null)
throw new ArgumentException("SampleDirichlet:: Results placeholder is null");
if (alpha <= 0.0)
throw new ArgumentException($"SampleDirichlet:: alpha {alpha} is non-positive");
int n = rn.Length;
if (n == 0)
throw new ArgumentException("SampleDirichlet:: Results placeholder is of zero size");
var gamma = new Gamma(alpha, 1.0);
double sum = 0.0;
for(int k = 0; k != n; ++k) {
double v = gamma.Sample();
sum += v;
rn[k] = v;
}
if (sum <= 0.0)
throw new ApplicationException($"SampleDirichlet:: sum {sum} is non-positive");
// normalize
sum = 1.0 / sum;
for(int k = 0; k != n; ++k) {
rn[k] *= sum;
}
}
static bool SampleBoundedDirichlet(double alpha, double sum, double lo, double hi, double[] rn)
{
if (rn == null)
throw new ArgumentException("SampleDirichlet:: Results placeholder is null");
if (alpha <= 0.0)
throw new ArgumentException($"SampleDirichlet:: alpha {alpha} is non-positive");
if (lo >= hi)
throw new ArgumentException($"SampleDirichlet:: low {lo} is larger than high {hi}");
int n = rn.Length;
if (n == 0)
throw new ArgumentException("SampleDirichlet:: Results placeholder is of zero size");
double mean = sum / (double)n;
if (mean < lo || mean > hi)
throw new ArgumentException($"SampleDirichlet:: mean value {mean} is not within [{lo}...{hi}] range");
SampleDirichlet(alpha, rn);
bool rc = true;
for(int k = 0; k != n; ++k) {
double v = lo + (mean - lo)*(double)n * rn[k];
if (v > hi)
rc = false;
rn[k] = v;
}
return rc;
}
static void Main(string[] args)
{
double[] rn = new double [30];
double lo = -50.0;
double hi = 50.0;
double alpha = 10.0;
double sum = 300.0;
for(int k = 0; k != 1_000; ++k) {
var q = SampleBoundedDirichlet(alpha, sum, lo, hi, rn);
Console.WriteLine($"Rng(BD), v = {q}");
double s = 0.0;
foreach(var r in rn) {
Console.WriteLine($"Rng(BD), r = {r}");
s += r;
}
Console.WriteLine($"Rng(BD), summa = {s}");
}
}
}
int[] hh = new int[100]; // histogram allocated
var s = 1.0; // step size
int k = 0; // good samples counter
for( ;; ) {
var q = SampleBoundedDirichlet(alpha, sum, lo, hi, rn);
if (q) // good sample, accept it
{
var v = rn[0]; // any index, 0 or 29 or ....
var i = (int)((v - lo) / s);
i = System.Math.Max(i, 0);
i = System.Math.Min(i, hh.Length-1);
hh[i] += 1;
++k;
if (k == 100000) // required number of good samples reached
break;
}
}
for(k = 0; k != hh.Length; ++k)
{
var x = lo + (double)k * s + 0.5*s;
var v = hh[k];
Console.WriteLine($"{x} {v}");
}