求助实现一下程序，小弟拜谢，很急。

shiyiming

The algorithm is based on Newton's method for solving the following system for A,
the vector of fitted numbers of cases at each nonzero exposure level. We have an equation
for each observed exposure level,
Lx + log(M1-A+) + log(Nx –Ax)- logAx -log (N0 – M1+A+)= 0,
where A+ is the sum of the elements of A (note that Ao is not in A, since Ao = M1 - A+). An
initial value A(0) may be the crude observed totals, if available, or the null expected value
M1N/n, where N is the vector of Nx for x≠ 0 and n is the total number of subjects in the
data. The algorithm may diverge from poor starting values; in our experience, convergence
was always achieved by starting with the crude observed totals rather than the null expected
values.
At iteration i, define
A0(i)=M1-A+(i)
Cx(i)=1/ Ax(i)+ 1/ (N-Ax(i)) for all x (including x = 0),
ex(i)=Lx+logA0(i)+log(Nx-Ax(i))-logAx(i)-log(N0-A0(i)) for x≠0，
e(i)=the vector of ex(i),
H(i)= the matrix with for Cx(i)+C0(i) for on-diagonal elements and C0(i) for all off-diagonal
elements, and
A(i+1)=A(i)+( H(i))-1 e(i)
Convergence is achieved when the increments become negligible relative to the Ax(i) and
N - Ax(i)  for all x. For person-time data, the equations become
Lx + log(M1- A+) +logNx - logAx – logN0 = 0;
the expression for ex(i)is similarly modified; and Cx(i) becomes 1/Ax(i).