标签: numerical-methods

我一直在研究log-sum-exp问题.我有一个以对数形式存储的数字列表,我希望将它们相加并存储在对数中.

天真的算法是

def naive(listOfLogs):
    return math.log10(sum(10**x for x in listOfLogs))

许多网站包括: 在C中实现logsumexp？ 和 http://machineintelligence.tumblr.com/post/4998477107/ 推荐使用

def recommend(listOfLogs):
    maxLog = max(listOfLogs)
    return maxLog + math.log10(sum(10**(x-maxLog) for x in listOfLogs))

又名

def recommend(listOfLogs):
    maxLog = max(listOfLogs)
    return maxLog + naive((x-maxLog) for x in listOfLogs)

我不明白的是,如果推荐算法更好,为什么我们应该递归调用它？那能提供更多的好处吗？

def recursive(listOfLogs):
    maxLog = max(listOfLogs)
    return maxLog + recursive((x-maxLog) for x in listOfLogs)

虽然我问是否还有其他技巧可以让这个计算在数值上更稳定？

我试图在不使用 math.h 的情况下找到 e ^{x 。}当x大于或小于 ~\xc2\xb120时，我的代码给出错误的答案。我尝试将所有双精度类型更改为长双精度类型，但它在输入上产生了一些垃圾。

我的代码是：

#include <stdio.h>\n\ndouble fabs1(double x) {\n    if(x >= 0){\n        return x;\n    } else {\n        return x*(-1);\n    }\n}\n\ndouble powerex(double x) {\n    double a = 1.0, e = a;\n    for (int n = 1; fabs1(a) > 0.001; ++n) {\n        a = a * x / n;\n        e += a;\n    }\n    return e;\n}\n\nint main(){\n    freopen("input.txt", "r", stdin);\n    freopen("output.txt", "w", stdout);\n    int n;\n    scanf("%d", &n);\n    for(int i = 0; i<n; i++) {\n        double …

我希望生成一个伪随机数/置换,它在一个范围内"占据"一个完整的周期或整个周期.通常,"线性同余发生器"(LCG)可用于生成此类序列,使用如下公式:

X = (a*Xs+c) Mod R

其中Xs是种子,X是结果,a和c是相对的素数常数,R是最大值(范围).

(通过完整周期/完整周期,我的意思是可以选择常数,使得任何X在一些随机/置换序列中仅出现一次,并且将在0到R-1或1到R的范围内).

LCG几乎满足了我的所有需求.我对LCG的问题是奇数/偶数结果的非随机性,即:对于种子Xn,结果X将交替奇数/偶数.

问题:

1. 有没有人知道如何创造类似的东西,不会交替奇/偶？

2. 我相信可以建立'复合LCG',但我没有细节.有人可以举一个这个CLCG的例子吗？

3. 是否有替代公式可能符合上述细节和下面的约束？

约束:

4. 我想要一些基于种子的简单公式.即:为了获得下一个数字,我提供种子并获得置换序列中的下一个"随机数".具体来说,我不能使用预先计算的数组.(见下一点)
5. 序列绝对必须是'全周期/完整周期'
6. 范围R可能是几百万甚至32比特/ 40亿.

7. 计算不应该溢出并且有效/快速,即:没有大的指数或几十个乘法/除法.

8. 序列不一定非常随机或安全 - 我不需要加密随机性(但如果可行则可以使用它),只需要"好"随机性或明显随机性,没有奇数/偶数序列.

任何想法都赞赏 - 提前感谢.

更新:理想情况下,Range变量可能不是2的精确幂,但在任何一种情况下都应该有效.

R可以使用样条线库中的splinefun()生成样条函数.但是,我需要在其一阶和二阶导数上评估此函数.有没有办法做到这一点？

例如

library(splines)
x <- 1:10
y <- sin(pi/x) #just an example
f_of_x <- splinefun(x,y)

如何评估f'(x)为x的向量？

有没有办法在PHP中获取浮点数的二进制表示？像Java的Double.doubleToRawLongBits()之类的东西.

给定一个正浮点数,我想获得最大可表示的浮点数,该数字小于该数.在Java中,我可以这样做:

double x = Double.longBitsToDouble(Double.doubleToRawLongBits(d) - 1);

但我在PHP中没有看到类似的东西.

我正在尝试将时间戳(仅秒的小数部分)从纳秒(10 ^ -9秒的单位)重新缩放到NTP时间戳的下半部分(单位为2 ^ -32秒).实际上这意味着乘以4.2949673.但我需要在没有浮点数学的情况下完成它,并且不使用大于32位的整数(事实上,我实际上是为8位微控制器编写这个,所以即使32位数学也很昂贵,特别是划分).

我已经提出了一些工作得相当好的算法,但我没有真正的数值方法基础,所以我很欣赏任何有关如何改进它们的建议,或任何其他更准确的算法. /或更快.

算法1

uint32_t intts = (ns >> 16) * 281474 + (ns << 16) / 15259 + ns / 67078;

选择前两个常数略微下冲,而不是超过正确的数字,并且根据经验确定最终因子67078以纠正此问题.产生+/- 4 NTP单位的正确值,结果为+/- 1 ns - 可接受,但残差随之变化ns.我想我可以添加另一个术语.

算法2

uint32_t ns2 = (2 * ns) + 1;
uint32_t intts = (ns2 << 1)
  + (ns2 >> 3) + (ns2 >> 6) + (ns2 >> 8) + (ns2 >> 9) + (ns2 >> 10)
  + (ns2 >> 16) + (ns2 >> 18) + (ns2 …

我有两个长度为n和m的序列.每个都是形式(x,y)的一系列点,并表示图像中的曲线.我需要找到不同(或类似)这些序列给出的事实

1. 一个序列可能比另一个序列长(即,一个序列可以是另一个序列的一半或四分之一,但如果它们跟踪大致相同的曲线,则它们是相同的)

2. 这些序列可能是相反的方向(即序列1从左到右,而序列2从右到左)

   我研究了Levenshtein之类的一些差异估计以及蛋白质折叠的结构相似性匹配中的编辑距离,但它们似乎都没有.我可以编写自己的暴力方法,但我想知道是否有更好的方法.

谢谢.

<math.h>log(1+x)为doubles的计算提供了更准确的方法.

是否有类似的精确计算方式log(1-x)？

我问的原因是因为我试图在对数空间中做一些工作以获得更高的精度(我主要是乘法和求和数非常接近于零).我发现编写一个log( exp(log_of_a) + exp(log_of_b) ) = log( a + b )使用的函数很容易log1p.我试图为差异做一个类似的功能:

log( exp(log_of_a) - exp(log_of_b) ) = log( a - b ) a > b当然,在哪里.

从本质上讲,只要没有log_a或log_b == -inf,函数会简单的返回:

return log( 1 - exp(log_b-log_a) ) + log_a;

在我的log_add功能中,我最终得到了一个log( 1 + ... ),所以我使用了log1p.但我在这里log( 1 - ... ).为了以防万一,我甚至用谷歌搜索log1m,但没有运气......

当参数x在范围内时[-inf, 1),我可以简单地使用log1p(-x)(给出我的断言a > …

问题的简要说明:我使用Newton Raphson算法在多项式中进行根发现,但在某些情况下不起作用.为什么？

我从"c ++中的数字配方"中得到了一个N​​ewton Raphson混合算法,如果New-Raph没有正确收敛(具有低导数值或者如果收敛速度不快),则将其一分为二.

我用几个多项式检查了算法并且它有效.现在我在我的软件内部进行测试,并且我总是遇到一个特定多项式的错误.我的问题是,我不知道为什么这个多项式只是没有得到结果,而其他人做的很多.因为我想改进任何多项式的算法,需要知道哪一个是没有收敛的原因,所以我可以正确对待它.

下面我将发布我可以提供的有关算法和我遇到错误的多项式的所有信息.

多项式:

f(t)= t^4 + 0,557257315256597*t^3 - 3,68254086033178*t^2 +
+ 0,139389107255627*t + 1,75823776590795

它的第一个衍生物:

 f'(t)= 4*t^3 + 1.671771945769790*t^2 - 7.365081720663563*t + 0.139389107255627

情节:

根(由Matlab提供):

  -2.133112008595826          1.371976341295347          0.883715461977390 
  -0.679837109933505

算法:

double rtsafe(double* coeffs, int degree, double x1, double x2,double xacc,double xacc2)
    {
    int j;
    double df,dx,dxold,f,fh,fl;
    double temp,xh,xl,rts;
    double* dcoeffs=dvector(0,degree);
    for(int i=0;i<=degree;i++)
        dcoeffs[i]=0.0;
    PolyDeriv(coeffs,dcoeffs,degree);
    evalPoly(x1,coeffs,degree,&fl);
    evalPoly(x2,coeffs,degree,&fh);
    evalPoly(x2,dcoeffs,degree-1,&df);
if ((fl > 0.0 && fh > 0.0) || (fl < 0.0 && fh < 0.0))
    nrerror("Root must …

TL; DR：问题是关于乘法精度

我必须将矩阵相乘A（100x8000），B（8000x27）和C（27x1）。

由于矩阵B和C是常数并且A是可变的，因此我更喜欢将其计算为：ABC = np.dot(A, np.dot(B, C))。但是我不知道，它可能是数值恶化（在以下方面准确性）比np.dot(np.dot(a, B), C)。

可能重要的是：矩阵A并B包含8000个分别（分别）100个和27个相关特征的样本。

乘法在数值上（在精度上）是否最优？如果是，我该如何确定？

特殊情况

可以假定A和B矩阵都是非负的。此外：

C = np.linalg.solve(cov(B, k), X)

其中X是27x1矩阵，其中包含27个（可能是相关的）未知分布的随机变量cov = lambda X, k: np.dot(X.T, X) + k * np.eye(X.shape[1])，并且k是将表达式最小化的非负常数：

sum((X[i, 0] - np.dot(np.dot(B[:, [i]].T, drop(B, i)),
                      np.linalg.solve(cov(drop(B, i), k),
                                      np.delete(X, …