1

我有一个算法：

  • 接受两个输入 (x, y)，其中 x 和 y 是两个独立泊松分布中的“lambda”变量 ( http://en.wikipedia.org/wiki/Poisson_distribution#Definition )
  • 计算每个 x 和 y 的泊松向量，
  • 计算两个泊松向量的矩阵乘积
  • 返回 [sum(upper_quadrant), sum(lower_quadrant), sum(diagonal)]

所以：

import math, random

def poisson(m, n):
    p=math.exp(-m)
    r=[p]
    for i in range(1, n):
        p*=m/float(i)
        r.append(p)
    return r

def simulate(mx, my, n):
    r=[0.0 for i in range(3)]
    px, py = (poisson(mx, n), 
              poisson(my, n))
    for i in range(n):
        for j in range(n):
            if i > j:
                k=0
            elif i < j:
                k=1
            else:
                k=2
            r[k]+=px[i]*py[j]
    return r

现在我需要求解 x 和 y，给定一组特定的输出。我已经破解了以下求解器函数：

def solve(p, n, generations, decay, tolerance):
    def rms_error(X, Y):
        return (sum([(x-y)**2 
                     for x, y in zip(X, Y)])/float(len(X)))**0.5
    def calc_error(x, y, n, target):
        guess=simulate(x, y, n)
        return rms_error(target, guess)    
    q=[0, 0]
    err=calc_error(math.exp(q[0]), math.exp(q[1]), n, p)
    bestq, besterr = q, err
    for i in range(generations):
        if besterr < tolerance:
            break
        q=list(bestq)
        if random.random() < 0.5:
            j=0
        else:
            j=1
        fac=((generations-i+1)/float(generations))**decay
        q[j]+=random.gauss(0, 1)*fac
        err=calc_error(math.exp(q[0]), math.exp(q[1]), n, p)
        if err < besterr:
            bestq, besterr = q, err
            # print (i, bestq, besterr)
    q, err = [math.exp(q) for q in bestq], besterr
    return (i, q, err)

哪个有效，但似乎需要相对大量的尝试才能返回一个不太优化的响应：

if __name__=="__main__":
    p, n = [0.5, 0.2, 0.3], 10
    q, err = solve_match_odds(p, n, 
                              generations=1000,
                              decay=2,
                              tolerance=1e-5)
    print q
    print simulate_match_odds(q[0], q[1], n)
    print (i, err)

和：

justin@justin-ThinkPad-X220:~/work/$ python solve.py 
[0.5, 0.2, 0.3]
[0.5000246335218251, 0.20006624338256798, 0.29990837191131686]
(999, 6.680993630511076e-05)
justin@justin-ThinkPad-X220:~/work/$

我不是 CS 专业的，我觉得我在这里错过了一系列搜索文献。有人可以建议一种更好的方法来搜索二维空间中的变量，例如这样吗？

谢谢。

4

1 回答 1

0

您可以使用scipy.optimize库来解决这个问题：

import math, random

def poisson(m, n):
    p=math.exp(-m)
    r=[p]
    for i in range(1, n):
        p*=m/float(i)
        r.append(p)
    return r

def simulate(mx, my, n):
    r=[0.0 for i in range(3)]
    px, py = (poisson(mx, n), 
              poisson(my, n))
    for i in range(n):
        for j in range(n):
            if i > j:
                k=0
            elif i < j:
                k=1
            else:
                k=2
            r[k]+=px[i]*py[j]
    return r

from scipy import optimize

Target, N = [0.5, 0.2, 0.3], 10

def error(p, target, n):
    r = simulate(p[0], p[1], n)
    return np.sum(np.subtract(target, r)**2)

r = optimize.fmin(error, (0, 0), args=(Target, N))
print simulate(r[0], r[1], n)

输出：

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 67
         Function evaluations: 125
[0.49999812501285623, 0.20000418001288445, 0.2999969464616799]
于 2013-03-21T07:56:27.210 回答