带cvxpy凸似然函数的极大似然

2024-05-14 06:34:57 发布

您现在位置:Python中文网/ 问答频道 /正文

我试图用cvxpy找到凸似然函数的最优值。变量表示分段线性函数在预先指定的x值处的y值,x值表示一天内的时间。我想知道是否有经验的人对如何正确实现下面的似然函数有意见

The log-likelihood function is defined like this

The function ^{} in the likelihood function is defined as this

Here is a link to the full research paper that I am trying to implement,但我怀疑你可能需要某种大学/研究人员的权限来阅读它

这是我写的代码。下面是我目前的错误和堆栈跟踪,但我怀疑还有更大的概念错误-但我现在很迷茫,所以请问一些不清楚的问题,我将尝试进一步说明

def ll(values, times, y):
  N = pd.DataFrame({
    'count_at_time': values,
    'timestamp': pd.to_datetime(times.flatten()),
  })
  N = N.set_index('timestamp')
  dates = N.index.normalize().drop_duplicates().strftime('%Y-%m-%d').to_list()
  f = np.vectorize(lambda date: sum_term(N, y, date))
  terms = f(dates)
  n = len(dates)
  ll = (1.0/n) * cp.sum(terms)
  return ll

# N[i](s) is the cumulative arrivals on day i at time s
def sum_term(N, y, date):
  t = Intervals.get(date)
  # hack: add 1 second to simulate open interval
  t0 = t[0].start_time.to_numpy().astype(int) + 1000000000
  tp = t[-1].end_time.to_numpy().astype(int)

  def t1(s):
    v1 = rate(y, t, s)
    v2 = at(N, date, s)

    return cp.log(v1 * v2)

  xvals = np.linspace(t0, tp, num=50, endpoint=False)

  yvals = list(map(lambda s: t1(s), xvals))
  v1 = np.trapz(yvals, x=xvals)

  yvals2 = list(map(lambda s: rate(y, t, s), xvals))
  v2 = np.trapz(yvals2, x=xvals)

  return v1 - v2

def rate(y, t, s):
  thing = t.start_time.astype(np.int64).values

  # end of day hack
  thing = np.append(thing, t[23].end_time.to_numpy().astype(int))

  i = np.searchsorted(thing, s, side='left')
  m = y[i-1]

  k = (y[i]- y[i-1]) / (thing[i] - thing[i-1])
  x = (s - thing[i-1]) # difference between time s and time at interval point t[i-1]

  return m + k * x

# data is:
#   - arrival time grouped by day
#   - count of arrivals at that time that day

data = ExperimentsService.run().data

N = data.count_at_time.values

X = np.array(data.index.astype(np.int64))
print("X", X, type(X), np.shape(X))

y = cp.Variable(25)
print(y.shape)

objective = cp.Minimize(-ll(N, X, y))
constraints = [ y >= 0, y[0] == y[-1] ]
prob = cp.Problem(objective, constraints)
val = prob.solve()
solution = y.value
print("Optimal value", val)
print("Optimal variables", solution)
Traceback (most recent call last):
  File "convex.py", line 132, in <module>
    objective = cp.Minimize(-ll(N, X, y))
  File "convex.py", line 34, in ll
    ll = (1.0/n) * cp.sum(terms)
  File "/Users/aron/.local/share/virtualenvs/project-nFRdcLaN/lib/python3.7/site-packages/cvxpy/atoms/affine/sum.py", line 109, in sum
    return Sum(expr, axis, keepdims)
  File "/Users/aron/.local/share/virtualenvs/project-nFRdcLaN/lib/python3.7/site-packages/cvxpy/atoms/affine/sum.py", line 39, in __init__
    super(Sum, self).__init__(expr, axis=axis, keepdims=keepdims)
  File "/Users/aron/.local/share/virtualenvs/project-nFRdcLaN/lib/python3.7/site-packages/cvxpy/atoms/axis_atom.py", line 31, in __init__
    super(AxisAtom, self).__init__(expr)
  File "/Users/aron/.local/share/virtualenvs/project-nFRdcLaN/lib/python3.7/site-packages/cvxpy/atoms/atom.py", line 43, in __init__
    self.args = [Atom.cast_to_const(arg) for arg in args]
  File "/Users/aron/.local/share/virtualenvs/project-nFRdcLaN/lib/python3.7/site-packages/cvxpy/atoms/atom.py", line 43, in <listcomp>
    self.args = [Atom.cast_to_const(arg) for arg in args]
  File "/Users/aron/.local/share/virtualenvs/project-nFRdcLaN/lib/python3.7/site-packages/cvxpy/expressions/expression.py", line 446, in cast_to_const
    return expr if isinstance(expr, Expression) else cvxtypes.constant()(expr)
  File "/Users/aron/.local/share/virtualenvs/project-nFRdcLaN/lib/python3.7/site-packages/cvxpy/expressions/constants/constant.py", line 44, in __init__
    self._value = intf.DEFAULT_INTF.const_to_matrix(value)
  File "/Users/aron/.local/share/virtualenvs/project-nFRdcLaN/lib/python3.7/site-packages/cvxpy/interface/numpy_interface/ndarray_interface.py", line 50, in const_to_matrix
    return result.astype(numpy.float64)
ValueError: setting an array element with a sequence.

Tags: toinpyprojectsharetimelocalnp