天道酬勤,学无止境

如何正确钳制贝克曼分布(How to properly clamp beckmann distribution)

问题

我正在尝试实现 Microfacet BRDF 着色模型(类似于 Cook-Torrance 模型),但我在本文中定义的贝克曼分布方面遇到了一些问题:https://www.cs.cornell.edu/~srm/publications /EGSR07-btdf.pdf

在此处输入图片说明

其中 M 是微面法线,N 是宏观面法线,ab 是 [0, 1] 之间的“硬度”参数。

我的问题是这种分布通常会返回非常大的值,尤其是当 ab 非常小时。

例如,贝克曼分布用于计算根据以下方程生成微面法线 M 的概率:

在此处输入图片说明

概率必须在范围 [0,1] 之间,那么如果贝克曼分布给我的值大小为 1000000000+,那么如何使用上面的函数获得该范围内的值?

那么有没有合适的方法来限制分布? 还是我误解了它或概率函数? 如果值超过 1,我曾尝试将其简单地钳位为 1,但这并没有真正给我想要的结果。

回答1

我和你有同样的问题。

如果你读

http://blog.selfshadow.com/publications/s2012-shading-course/hoffman/s2012_pbs_physics_math_notes.pdf

http://blog.selfshadow.com/publications/s2012-shading-course/hoffman/s2012_pbs_physics_math_notebook.pdf

你会注意到这是完全正常的。 从链接引用:

“贝克曼 Αb 参数等于 RMS(均方根)微面斜率。因此,其有效范围是从 0(不包含 –0 对应于完美镜像或狄拉克三角洲并导致贝克曼公式中的除以 0 错误)并达到任意高的值。值 1 没有特殊意义——这只是意味着 RMS 斜率是 1/1 或 45°。(...)"

还有另一个引用:

“微面方向的统计分布是通过微面正态分布函数 D(m) 定义的。与 F() 不同,D() 的值不限于介于 0 和 1 之间——尽管值必须是非负的,它们可以任意大(表示法线指向特定方向的微平面非常集中。(...)”

你应该在谷歌上搜索 Self Shadow 的 Physically Based Shading 课程,里面有很多有用的材料(每年有一篇博文:2010、2011、2012 和 2013)

受限制的 HTML

  • 允许的HTML标签:<a href hreflang> <em> <strong> <cite> <blockquote cite> <code> <ul type> <ol start type> <li> <dl> <dt> <dd> <h2 id> <h3 id> <h4 id> <h5 id> <h6 id>
  • 自动断行和分段。
  • 网页和电子邮件地址自动转换为链接。

相关推荐
  • How to properly clamp beckmann distribution
    I am trying to implement a Microfacet BRDF shading model (similar to the Cook-Torrance model) and I am having some trouble with the Beckmann Distribution defined in this paper: https://www.cs.cornell.edu/~srm/publications/EGSR07-btdf.pdf Where M is a microfacet normal, N is the macrofacet normal and ab is a "hardness" parameter between [0, 1]. My issue is that this distribution often returns obscenely large values, especially when ab is very small. For instance, the Beckmann distribution is used to calculate the probability of generating a microfacet normal M per this equation : A probability
  • 【基于物理的渲染(PBR)白皮书】(五)几何函数相关总结
    本文由@浅墨_毛星云 出品,首发于知乎专栏,转载请注明出处 文章链接: https://zhuanlan.zhihu.com/p/81708753 在基于物理的渲染中,几何函数(Geometry Function)是保证Microfacet BRDF理论上能量守恒,逻辑上自洽的重要一环。其描述了微平面自阴影的属性,表示具有半矢量法线的微平面(microfacet)中,同时被入射方向和反射方向可见(没有被遮挡的)的比例,即未被遮挡的m = h微表面的百分比。 在Microfacet Specular BRDF的D,G,F三项中,如果说法线分布函数是最核心的一项,那么几何函数则是核心的辅助项,而且是三项中最复杂的一项。 历史上主流的几何函数建模,按提出或归纳的时间进行排序,可以总结为: Smith [1967]V-cavity(Cook-Torrance)[1982]Schlick-Smith [1994]Neumann [1999]Kelemen [2001]Implicit [2010] 其中,Smith遮蔽函数(Smith masking function)是现在业界所采用的主流遮蔽函数,Eric Heitz在2014年[Heitz 2014]将其拓展为Smith联合遮蔽阴影函数(Smith Joint Masking-Shadowing Function),该函数具有四种形式:
  • 钳制真实(固定/浮点)值的最快方法?(Fastest way to clamp a real (fixed/floating point) value?)
    问题 有没有比使用if语句或三元运算符更有效的方法来钳制实数? 我想为双打和32位修订点实现(16.16)都这样做。 我并不是在要求可以处理这两种情况的代码。 它们将在单独的功能中处理。 显然,我可以做类似的事情: double clampedA; double a = calculate(); clampedA = a > MY_MAX ? MY_MAX : a; clampedA = a < MY_MIN ? MY_MIN : a; 或者 double a = calculate(); double clampedA = a; if(clampedA > MY_MAX) clampedA = MY_MAX; else if(clampedA < MY_MIN) clampedA = MY_MIN; 定点版本将使用功能/宏进行比较。 这是在代码的性能关键部分中完成的,因此我正在寻找一种尽可能高效的方法(我怀疑这会涉及位操作) 编辑:它必须是标准/便携式C,特定于平台的功能在这里没有任何意义。 另外, MY_MIN和MY_MAX与我要限制的值具有相同的类型(在上面的示例中为两倍)。 回答1 对于16.16表示形式,简单的三元级不可能在速度方面得到改善。 对于双打,因为您需要标准/便携式C,所以任何形式的位纠结都将以糟糕的结果结束。 即使有可能进行小提琴演奏(我对此表示怀疑)
  • 对数深度缓冲区OpenGL(Logarithmic Depth Buffer OpenGL)
    问题 我已经设法在OpenGL中实现对数深度缓冲区,主要是来自Outerra的文章(您可以在此处,此处和此处阅读)。 但是,我遇到了一些问题,并且不确定这些问题是否是使用对数深度缓冲区所固有的,或者是否有我无法想到的解决方法。 刚开始时,这就是我如何计算顶点着色器内的对数深度的方法: gl_Position = MVP * vec4(inPosition, 1.0); gl_Position.z = log2(max(ZNEAR, 1.0 + gl_Position.w)) * FCOEF - 1.0; flogz = 1.0 + gl_Position.w; 这就是我修复片段着色器中的深度值的方法: gl_FragDepth = log2(flogz) * HALF_FCOEF; 其中ZNEAR = 0.0001 , ZFAR = 1000000.0 , FCOEF = 2.0 / log2(ZFAR + 1.0) ,并且HALF_FCOEF = 0.5 * FCOEF 。 在我的情况下, C为1.0,以简化代码并减少计算量。 首先,我对获得的精度水平感到非常满意。 使用普通的深度缓冲(znear = 0.1,zfar = 1000.0),我可以在视距边缘进行很多z角战斗。 现在,将我的MUCH进一步设为znear:zfar,我将第二个地平面放置在第一个地平面以下0
  • 在Swift中将数字“限制”在两个值之间的标准方法(Standard way to “clamp” a number between two values in Swift)
    问题 鉴于: let a = 4.2 let b = -1.3 let c = 6.4 我想知道将这些值钳位到给定范围(例如0...5的最简单,最快捷的方法: a -> 4.2 b -> 0 c -> 5 我知道我可以执行以下操作: let clamped = min(max(a, 0), 5) 或类似的东西: let clamped = (a < 0) ? 0 : ((a > 5) ? 5 : a) 但是我想知道在Swift中是否还有其他方法可以做到这一点-特别是我想知道(并在SO上进行文档记录,因为在Swift中似乎没有关于夹紧数字的问题)中是否包含任何内容。 Swift标准库专门用于此目的。 可能没有,如果是的话,那也是我很乐意接受的答案。 回答1 斯威夫特4/5 Comparable/Strideable扩展类似于标准Swift库中的ClosedRange.clamped(to:_) -> ClosedRange 。 extension Comparable { func clamped(to limits: ClosedRange<Self>) -> Self { return min(max(self, limits.lowerBound), limits.upperBound) } } // Swift < 5.1 extension Strideable where
  • PBR Specular D的几何学含义
    PBR Specular D的几何学含义 jplee 2019.9.03 PBR( Physically Based Rendering )里成为大趋势的Cook-Torrance的Specular BRDF [1]定义如下: 这里面D函数是对于微表面(microfacet)的倾斜度的分布(Distribution)函数。本文档将对于本函数在几何学上的意义进行说明。 一个表面(surface)是由几个微表面构成,这是PBR理论的根本。而且构成表面的各微表面都有着不同的法线(normal)。另外微表面上发生了镜面反射(mirror reflection, 正反射)。 很多文档里面将这样的微表面的法线用m表示,也会用h表示。应该是microfacet normal和half vector的略写。那么图1里面表面和微表面的关系就显而易见了。 图1:由微表面构成的表面[2]. 微表面的法线用m来表示很容易理解,但是用h(half vector)来表示就有点费解了。Burley和Cook-Torrance [1, 3]一致记述了用半向量(half vector)来表现微表面的法线。 这里最好是用半向量的思路来理解这个含义。上面我们已经假设微表面能够镜面反射,现在我们来反过来看一下这个问题。 表面内,能够给我们展示镜面反射的、带有倾斜度的微表面究竟有多少呢? 这就是我们要通过Specular
  • 将颜色值从浮点数0..1转换为字节0..255(Converting color value from float 0..1 to byte 0..255)
    问题 将颜色值从浮点数转换为字节的正确方法是什么? 起初我以为b=f*255.0应该这样做,但现在我在想,在这种情况下,只有精确的1.0会转换为255 ,而0.9999已经是254 ,这可能不是我想要的... 似乎b=f*256.0会更好,除了在精确1.0的情况下不希望产生256的情况。 最后,我使用了这个: #define F2B(f) ((f) >= 1.0 ? 255 : (int)((f)*256.0)) 回答1 1.0是唯一可能出错的情况,因此请单独处理该情况: b = floor(f >= 1.0 ? 255 : f * 256.0) 同样,可能有必要强迫f实际上为0 <= f <= 1,以避免由于舍入误差(例如f = 1.0000001)而导致的错误行为。 f2 = max(0.0, min(1.0, f)) b = floor(f2 == 1.0 ? 255 : f2 * 256.0) 替代安全解决方案: b = (f >= 1.0 ? 255 : (f <= 0.0 ? 0 : (int)floor(f * 256.0))) 或者 b = max(0, min(255, (int)floor(f * 256.0))) 回答2 我总是做round(f * 255.0) 。 不需要测试(特殊情况为1)和/或其他答案。 对于您的目的而¨
  • Directx11进阶教程PBR(2)之BRDF
    PBR单位 电磁辐射的测量是辐射度学。有很多的辐射量单位来描述物体表面接受的光照,但是图形渲染经常用RGB三色谱来作为单位, 也就是float3. BRDF基本属性 普通认为, 光照shading是在局部进行的, 在上面PBR物理现象的图17可以看出来。这种情况下,我们用 V 代表作为从物体表面观察点到人眼睛的单位方向向量,用 L 代表物体表面观察点的入射光的反方向单位向量。物体对光的交互关系用双向反射分布函数(BRDF)来描述,也就是有关于:“一个物体表面点在给定的入射光照方向L和给定的观察方向V,会表现如何的函数。如下图所示: 当然,入射光照方向L和观察方向V通常都是认定为表面之上的, 也就是dot(n, L)和dot(n, V)都是非负数的。 BRDF现象可以表达为两种方式,下面两种都是正确有效的 (1)多束不同方向的光在一个表面点的反射方向都为V。 (2)一束特定方向的入射光反射出多个方向的光。 (这里在IBL的specular重要性采样中计算可以用上面的理论来解释) BRDF计算采用的是光谱单位,但是上面我们提过了图形渲染通常采用的RGB三色谱单位,所以BRDF里面代入的各种光照值也是RGB三色谱值。 BRDF在光照方程式中使用,方程式: 这个光照方程式告诉我们一个点的表现等于沿着该点各种方向的入射光照的微积分计算. BRDF函数在光照计算中,遵循两大定律: (1
  • 如何转换RGB-> YUV-> RGB(双向)(How to convert RGB -> YUV -> RGB (both ways))
    问题 我想要一对转换算法,一个从RGB到YUV,另一个从YUV到RGB,它们彼此相反。 也就是说,往返转换应保持该值不变。 (如果愿意,请用Y'UV,YUV,YCbCr,YPbPr替换YUV。) 这样的事情存在吗? 如果是这样,那是什么? 发布的解决方案(如何在C / C ++中执行RGB-> YUV转换?,http://www.fourcc.org/fccyvrgb.php,http://en.wikipedia.org/wiki/YUV)只是反演(当省略钳位到[0,255]时,两个3x3矩阵是相反的)。 但是,忽略钳位会产生负亮度之类的东西,这会对YUV空间中的图像处理造成极大的破坏。 保持钳位会使转换成为非线性,这使得定义逆转换变得棘手。 回答1 是的,存在可逆转换。 equasys GmbH发布了从RGB到YUV,YCbCr和YPbPr的可逆转换,并解释了每种情况适合的情况,这种钳制的真正意义以及参考文献链接。 (就像一个很好的答案)。 对于我自己的应用程序(jpg图像,而不是模拟电压),YCbCr是合适的,因此我为这两个转换编写了代码。 的确,对于许多图像,反复来回的值相差不到256分之一。 前后图像在视觉上无法区分。 PIL的色彩空间转换YCbCr-> RGB因提及equasys的网页而倍受赞誉。 其他可能无疑会提高equasys的准确性和简洁性的答案: https:/
  • 如何使用mysql用户定义函数生成高斯分布(how to generate a gaussian distribution using mysql user-defined function)
    问题 我喜欢使用MySQL进行定量分析和统计。 我想使用以下形式的MySQL用户定义函数:sample_gaussian(mean,stdev),该函数返回从具有用户输入参数均值和标准差的高斯分布采样的单个随机值。 MySQL已经有一个函数rand()返回一个随机数,所以我只需要知道一些伪代码来约束/转换该值,以便它进入正确的分布。 有什么建议么? 顺便说一句,这是我的第一个stackoverflow问题,因此,如果这个问题在本网站上吸引了太多用户,请原谅我。 回答1 为了回答我自己的问题,这是一个MySQL用户定义函数,该函数返回具有给定均值和标准差的,从高斯分布中采样的单个随机值。 DROP FUNCTION IF EXISTS gauss; DELIMITER // CREATE FUNCTION gauss(mean float, stdev float) RETURNS float BEGIN set @x=rand(), @y=rand(); set @gaus = ((sqrt(-2*log(@x))*cos(2*pi()*@y))*stdev)+mean; return @gaus; END // DELIMITER ; 为了验证这实际上是在返回高斯分布,您可以生成一系列的高斯分布,然后绘制直方图: create temporary table temp (id
  • 如何正确触发对链接的 sql server 的插入?(How to properly trigger an insert to a linked sql server?)
    问题 我在两个不同的 sql 服务器中有同一个表(一个是 SqlServer 2000,另一个是 2008)。 我正在使用 sql 服务器管理工​​作室。 我希望每次在 SqlServer 2000 表 (Table_1) 中的表上进行插入时都会发生一个触发器,并且该记录也将插入到 SqlServer 2008 表(也是 Table_1)中的同一个表中。 sql server 2008 被定义为链接服务器,并且可以使用 sql server management studio 从 2000 db 连接运行查询并在 2008 db 上执行插入。 触发器定义如下: ALTER TRIGGER [dbo].[trgrTable] ON [dbo].[Table_1] AFTER INSERT AS BEGIN INSERT INTO [TLVSQL].[AVI_DEV].[dbo].[Table_1](ID, Value) SELECT INSERTED.ID AS ID, INSERTED.Value AS Value FROM INSERTED END [TLVSQL].[AVI_DEV]是 2008 数据库名称。 但是每次我在 2000 表上执行插入时,我都会收到一条消息,指出由于“sqloledb 无法启动分布式事务链接服务器......”,插入失败。
  • 如何从整数范围生成正态分布的随机数?(How to generate normally distributed random from an integer range?)
    问题 给定整数范围的开始和结尾,如何计算该范围之间的正态分布随机整数? 我意识到正态分布变成-+无穷大。 我猜尾巴可能会被截断,因此当在范围之外计算随机数时,请重新计算。 这会提高整数在该范围内的可能性,但是只要可以容忍此效果(<5%),就可以了。 public class Gaussian { private static bool uselast = true; private static double next_gaussian = 0.0; private static Random random = new Random(); public static double BoxMuller() { if (uselast) { uselast = false; return next_gaussian; } else { double v1, v2, s; do { v1 = 2.0 * random.NextDouble() - 1.0; v2 = 2.0 * random.NextDouble() - 1.0; s = v1 * v1 + v2 * v2; } while (s >= 1.0 || s == 0); s = System.Math.Sqrt((-2.0 * System.Math.Log(s)) / s); next_gaussian = v2 * s
  • 如何生成R中给定的均值,SD,偏斜和峰度的分布?(How to generate distributions given, mean, SD, skew and kurtosis in R?)
    问题 是否可以在R中生成均值,SD,偏斜和峰度已知的分布? 到目前为止,最好的方法似乎是创建随机数并相应地对其进行转换。 如果有适合生成特定发行版的软件包,而我可能还没有找到它。 谢谢 回答1 SuppDists软件包中包含Johnson发行版。 约翰逊将为您提供与力矩或分位数相匹配的分布。 其他评论是正确的,4个时刻不分配。 但是约翰逊肯定会尝试。 这是一个将Johnson拟合到一些样本数据的示例: require(SuppDists) ## make a weird dist with Kurtosis and Skew a <- rnorm( 5000, 0, 2 ) b <- rnorm( 1000, -2, 4 ) c <- rnorm( 3000, 4, 4 ) babyGotKurtosis <- c( a, b, c ) hist( babyGotKurtosis , freq=FALSE) ## Fit a Johnson distribution to the data ## TODO: Insert Johnson joke here parms<-JohnsonFit(babyGotKurtosis, moment="find") ## Print out the parameters sJohnson(parms) ## add the Johnson
  • C ++ 11中的随机数生成:如何生成,它如何工作? [关闭](Random number generation in C++11: how to generate, how does it work? [closed])
    问题 在这里很难说出要问什么。 这个问题是模棱两可,含糊,不完整,过于宽泛或夸张的,因此不能以目前的形式合理地回答。 如需帮助澄清此问题以便可以重新打开,请访问帮助中心。 9年前关闭。 最近,我遇到了一种在C ++ 11中生成随机数的新方法,但是无法消化我所读到的论文(这是什么引擎,数学术语,如“分布” ,即“产生的所有整数同等可能”)。 所以任何人都可以解释一下 这些是什么? 他们是什么意思? 怎么产生的? 他们如何工作? ETC 您可以在一个有关随机数生成的常见问题解答中将其全部调用。 回答1 这个问题太宽泛,无法给出完整的答案,但让我挑几个有趣的观点: 为什么“同样可能” 假设您有一个简单的随机数生成器,该生成器以相等的概率生成数字0、1,...,10(将其视为经典的rand() )。 现在,您需要0、1、2范围内的随机数,每个概率均等。 您的下意识反应是采用rand() % 3 。 但是,等等,余数0和1比余数2更频繁地出现,所以这是不正确的! 这就是为什么我们需要适当的分布,该分布需要一个统一的随机整数源并将其转变为我们想要的分布,例如示例中的Uniform[0,2] 。 最好把它留给一个好的图书馆! 引擎 因此,所有随机性的核心是一个好的伪随机数生成器,该生成器生成一个在一定间隔内均匀分布的数字序列,理想情况下,其周期很长。 rand()的标准实现通常不是最好的
  • EPB功能安全笔记(5):EPB系统软件接口定义
    本文要点在上文中基于VDA 305对EPB系统架构进行了阐述,并对系统架构中的各个软件模块的功能做了说明。其中,绿色部分属于Brake Assy,包括电机控制软件(PBC, parking brake controller)和制动卡钳;蓝色部分属于ESC Assy,又称Brake Host,提供EPB ECU和供电电源、CAN 通讯接口、硬线接口等外围设备和EPB宿主软件。Brake Assy和ESC Assy一起构成完整的EPB系统。EPB系统架构,绿色:Brake Assy; 蓝色:ESC Assy定义EPB系统架构是为了接下来分析技术安全要求,而技术安全要求最终需要落实到具体的交互信号上的。从这个角度上看,上图所示只是一个基于EPB功能需求初步定义的框架,框架中各个软件模块间的接口还需要明确定义。所以本文将对这部分内容进行补充说明。1.软件接口定义从EPB系统框架图中可以看出,确定各个软件模块间的接口,本质上就是确定Brake Assy的PBC软件和ESC Assy的Host软件间的接口。下面将从这个角度出发定义PBC软件和ESC Host软件的接口。1.1.激活请求接口(Actuation Request Interface)激活请求接口如下图所示,其中包含了SSM模块与PBC之间的双向通讯信号。PBC与Host之间的激活请求接口接口描述InterfaceTask and
  • EPB功能安全笔记(5):EPB系统软件接口定义
    本文要点在上文中基于VDA 305对EPB系统架构进行了阐述,并对系统架构中的各个软件模块的功能做了说明。其中,绿色部分属于Brake Assy,包括电机控制软件(PBC, parking brake controller)和制动卡钳;蓝色部分属于ESC Assy,又称Brake Host,提供EPB ECU和供电电源、CAN 通讯接口、硬线接口等外围设备和EPB宿主软件。Brake Assy和ESC Assy一起构成完整的EPB系统。EPB系统架构,绿色:Brake Assy; 蓝色:ESC Assy定义EPB系统架构是为了接下来分析技术安全要求,而技术安全要求最终需要落实到具体的交互信号上的。从这个角度上看,上图所示只是一个基于EPB功能需求初步定义的框架,框架中各个软件模块间的接口还需要明确定义。所以本文将对这部分内容进行补充说明。1.软件接口定义从EPB系统框架图中可以看出,确定各个软件模块间的接口,本质上就是确定Brake Assy的PBC软件和ESC Assy的Host软件间的接口。下面将从这个角度出发定义PBC软件和ESC Host软件的接口。1.1.激活请求接口(Actuation Request Interface)激活请求接口如下图所示,其中包含了SSM模块与PBC之间的双向通讯信号。PBC与Host之间的激活请求接口接口描述InterfaceTask and
  • 未烘烤的莫吉贝克(Unbaking mojibake)
    问题 当您错误地解码了字符时,如何识别原始字符串的可能候选者? Ä×èÈÄÄî▒è¤ô_üiâAâjâüâpâXüj_10òb.png 我知道这个图像文件名应该是一些日语字符。 但是由于对urllib引用/取消引用,编码和解码iso8859-1,utf8进行了各种猜测,所以我一直无法取消和获取原始文件名。 腐败是可逆的吗? 回答1 您可以使用chardet(通过pip安装): import chardet your_str = "Ä×èÈÄÄî▒è¤ô_üiâAâjâüâpâXüj_10òb" detected_encoding = chardet.detect(your_str)["encoding"] try: correct_str = your_str.decode(detected_encoding) except UnicodeDecodeError: print("Could not estimate encoding") 结果:时间测试视点(动画通过)_10秒(不知道这是否正确) 对于Python 3(编码为utf8的源文件): import chardet import codecs falsely_decoded_str = "Ä×èÈÄÄî¦è¤ô_üiâAâjâüâpâXüj_10òb" try: encoded_str = falsely_decoded
  • 在Matlab中的Wigner-Ville分布上应用窗口函数(To apply window function on Wigner-Ville Distribution in Matlab)
    问题 我们在这里思考如何创建重叠64的Hamming-64窗口。 h = hamming(64); h2 = hamming(38); h = conv(h, h2); 现在,我们正在考虑如何将此窗口函数应用于“时频工具箱”中Auger等人的Wigner-Ville分布函数的结果变量。 函数tfrwv.m没有用于窗口函数的任何参数。 所以我们有这些变量 [B,T,F] = tfrwv(data, 1:length(data), length(data)); 这是有关问题的一个答案,但并不完全相同。 有人说将window函数应用于结果 逐点相乘 h的尺寸是101x1的两倍,而T和F的尺寸是5001x1的两倍。 因此,如果逐点相乘,似乎需要对窗口矢量进行外推。 这里还有一个解释 在第二个代码块的大约一半处,我将窗口函数应用于缓冲的信号。 这实际上是窗口函数与时间序列数据的每个缓冲块的向量乘法。 我只是使用一个偷偷摸摸的对角矩阵技巧来有效地做到这一点。 如何将窗口函数应用于变量B,T和F ? 回答1 查看有关Wigner分布的本文。 从第8页到第11页。我认为tfr(indices,icol) = x(ti+tau,1) .* conj(x(ti-tau,xcol)); 在代码中实现公式(23)。 总和和指数部分等于tfr= fft(tfr); 。 或者
  • 在C / C ++中遵循正态分布生成随机数(Generate random numbers following a normal distribution in C/C++)
    问题 如何在C或C ++中按照正态分布轻松生成随机数? 我不想使用Boost。 我知道Knuth详细讨论了这个问题,但是我现在没有他的书。 回答1 有许多方法可以从常规RNG生成高斯分布的数字。 通常使用Box-Muller变换。 它会正确产生具有正态分布的值。 数学很简单。 您生成两个(均匀)随机数,然后对它们应用公式,就得到两个正态分布的随机数。 返回一个,并将另一个保存为下一个随机数请求。 回答2 C ++ 11 C ++ 11提供了std :: normal_distribution,这就是我今天要去的方式。 C或更旧的C ++ 以下是按复杂度升序排列的一些解决方案: 将0到1之间的12个均匀随机数相加并减去6。这将与正态变量的均值和标准差相匹配。 一个明显的缺点是范围限制为±6,这与真实的正态分布不同。 Box-Muller变换。 这已在上面列出,并且实现起来相对简单。 但是,如果您需要非常精确的样本,请注意,将Box-Muller变换与某些均匀生成器结合使用会遇到一种称为Neave Effect 1的异常现象。 为了获得最佳精度,我建议绘制制服并应用逆累积正态分布以得出正态分布变量。 这是一个很好的逆累积正态分布算法。 1. HR Neave,“关于将Box-Muller变换与乘性同余伪随机数生成器一起使用”,《应用统计》,1973年第22、92-97页 回答3
  • 用Python计算累积分布函数(CDF)(Calculate the Cumulative Distribution Function (CDF) in Python)
    问题 如何在python中计算累积分布函数(CDF)? 我想从我拥有的点数组(离散分布)中进行计算,而不是从scipy具有的连续分布中进行计算。 回答1 (我对问题的解释可能是错误的。如果问题是如何从离散的PDF转换为离散的CDF,则将np.cumsum除以合适的常数将在样本np.cumsum 。不np.cumsum ,则将数组的np.cumsum乘以点之间的距离即可。) 如果您有一个离散的样本数组,并且想知道样本的CDF,则可以对数组进行排序。 如果查看排序结果,您将意识到最小值代表0%,最大值代表100%。 如果您想知道分布的50%处的值,只需查看位于已排序数组中间的array元素即可。 让我们用一个简单的例子仔细看一下: import matplotlib.pyplot as plt import numpy as np # create some randomly ddistributed data: data = np.random.randn(10000) # sort the data: data_sorted = np.sort(data) # calculate the proportional values of samples p = 1. * np.arange(len(data)) / (len(data) - 1) # plot the sorted