Re: Problems with mkoctfile

[Henry F. Mollet]

Was NegBinSNP.cc supposed to be included with Octave 2.1.71?

octave:3> NegBinSNP(theta,data,{2})
error: `NegBinSNP' undefined near line 3 column 1

[~] -bash-2.05b 502$ locate NegBinSNP.cc
Produced zero results

The 'problem' was that my Octave 2.1.71 did not include NegBinSNP.cc given
below. After saving it and compiling it to produce (NegBinSNP.o) and
NegBinSNP.oct, it worked:
octave:3> NegBinSNP(theta,data,{2})
ans =

  -1.05123
  -1.25575
  -1.35181
  -1.45575
  -1.31630
  -1.15245
  -0.87225
  -1.36387
  -1.03389
  -1.28454

octave:4> pwd
/usr/local/include/octave-2.1.71/octave

Henry 


on 5/4/06 1:40 AM, Michael Creel at address@hidden wrote:

> Hello all,
> I'm having some trouble compiling a DLD file. I'm using octave 2.1.73 from
> Debian. The NegBinSNP.cc
> file compiles (mkoctfile NegBinSNP.cc) without apparent trouble, but doesn't
> work properly. This did
> work with Octave 2.1.71. To give an example of use and the problem:
> 
> octave:1> data = rand(10,3)
> data =
> 
>    0.350677  0.230081  0.861102
>    0.474211  0.608506  0.401342
>    0.914846  0.299243  0.147354
>    0.755405  0.468921  0.673132
>    0.654121  0.247208  0.013554
>    0.472270  0.889212  0.097541
>    0.968990  0.240560  0.959222
>    0.063041  0.171184  0.468724
>    0.921873  0.759646  0.639991
>    0.816379  0.137021  0.502845
> 
> octave:2> theta = rand(5,1)
> theta =
> 
>    0.60613
>    0.81519
>    0.34988
>    0.91822
>    0.91403
> 
> octave:3> NegBinSNP(theta,data,{2})
> *** glibc detected *** double free or corruption (!prev): 0x08a8c4f8 ***
> panic: Aborted -- stopping myself...
> attempting to save variables to `octave-core'...
> save to `octave-core' complete
> Aborted
> 
> The .cc file is attached. Any pointers greatly appreciated. Thanks, Michael
> // Copyright (C) 2004   Michael Creel   <address@hidden>
> //
> //  This program is free software; you can redistribute it and/or modify
> //  it under the terms of the GNU General Public License as published by
> //  the Free Software Foundation; either version 2 of the License, or
> //  (at your option) any later version.
> // 
> //  This program is distributed in the hope that it will be useful,
> //  but WITHOUT ANY WARRANTY; without even the implied warranty of
> //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
> //  GNU General Public License for more details.
> // 
> //  You should have received a copy of the GNU General Public License
> //  along with this program; if not, write to the Free Software
> //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
> 
> // This is the NegBinSNP loglikelihood function
> 
> #include <octave/oct.h>
> #include <octave/Cell.h>
> #include <octave/parse.h>
> #include <float.h>
> 
> 
> // define argument checks
> static bool
> any_bad_argument(const octave_value_list& args)
> {
> 
> int nargin = args.length();
> if (!(nargin == 3))
> {
> error("NegBinSNP: you must supply 3 arguments");
> return true;
> }
> if (!args(0).is_matrix_type())
> {
> error(": NegBinSNP: first argument vector theta (parameters)");
> return true;
> }
> if (!args(1).is_matrix_type())
> {
> error("NegBinSNP: second argument must be a matrix of data");
> return true;
> } 
> Cell otherargs (args(2).cell_value());
> if (!(otherargs.length() == 1))
> {
> error("NegBinSNP: the third argument must be cell");
> return true;
> }
> if (!otherargs(0).is_real_scalar())
> {
> error("NegBinSNP: third arg must cell containing negbin_type (1/2)");
> return true;
> }
> return false;
> }
> 
> DEFUN_DLD(NegBinSNP, args, ,"NebBinSNP likelihood function")
> {
> 
> // check the arguments
> if (any_bad_argument(args)) return octave_value_list();
> 
> ColumnVector theta = args(0).column_vector_value();
> Matrix data (args(1).matrix_value());
> 
> const int n = data.rows();
> const int k = data.columns() - 1;
> ColumnVector y = data.column(0);
> 
> Matrix x = data.extract_n(0, 1, n, k);
> 
> Cell otherargs (args(2).cell_value());
> int nb_type (otherargs(0).int_value());
> octave_value_list f_return;
> 
> int i, j;
> ColumnVector m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11;
> double t1, t6, t7, t8, t9, t10, t11, t12, t13, t14, t15, t16;
> double t17, t18, t19, t20, t21, t22, t23, t24, t25, t26, t27;
> double t28, t29, t30, t31, t32, t33, t34, t35, t36, t37, t38, t39;
> double t40, t41, t42, t43, t44;
> 
> const int kk = theta.rows();
> const int order = kk - k - 1; // that "1" is for alpha
> 
> if (order > 4)
> {
>         error("NegBinSNP: polynomial order is limited to 4 at most, at
> present");
>         return octave_value_list();
> }
> 
> 
> double alpha = DBL_EPSILON + exp(theta(k));
> ColumnVector beta(k);
> ColumnVector gam(order + 1); // 2nd degree has 3 coefficients, e.g.
> ColumnVector lam(n);
> ColumnVector psi(n);
> ColumnVector lampsi(n);
> ColumnVector logdensity(n);
> ColumnVector xbeta;
> ColumnVector normfactor(n);
> ColumnVector prediction(n);
> ColumnVector poly(n);
> 
> // how many moments are needed?
> // two times order for density,
> // two times order  plus 1 for prediction
> Matrix moments(n,2*order+2);
> moments.fill(1.0);
> 
> // extract beta as leading k elements of theta
> for(i = 0; i < k; i++)
> {
> beta(i) = theta(i);
> }
> 
> gam(0) = 1; // normalization
> // extract gam as last kk - k  - 1 elements of theta
> for(i = 0; i < order; i++)
> {
> gam(i+1) = theta(k+1+i) / pow(10.0, ((double) i) +1.0);
> }
> 
> xbeta = x*beta;
> 
> for (i=0; i < n; i++)
> {
> lam(i) = DBL_EPSILON + exp(xbeta(i));
> }
> 
> if (nb_type==1)
> {
> psi = lam / (alpha);
> } 
> if (nb_type==2)
> {
> psi = psi.fill(DBL_EPSILON + 1.0/alpha);
> } 
> 
> 
> // top-bound psi to bulletproof
> for (i=0; i < n; i++)
> {
> 
> psi(i) = (psi(i) < 10000) * psi(i) + (psi(i) > 10000)*10000;
> }
> 
> 
> // moments
> if (order==1)
> {
> for(i = 0; i < n; i++)
> {
> moments(i,1) = lam(i);
> 
> t11 = lam(i) + psi(i) ;
> t12 = psi(i)/t11 ;
> t13 = pow(t12, psi(i)) ;
> t14 = -psi(i) - 1.0 ;
> t15 = -(lam(i)/t11) + 1.0 ;
> t16 = log(t15) ;
> moments(i,2) = (t13*lam(i)*psi(i)*exp(t14*t16))/t11 - pow(t11,
> -2.0)*t13*t14*(lam(i)*lam(i)) \
> *psi(i)*exp(t16*(-psi(i) - 2.0)) ;
> 
> 
> t15 = lam(i) + psi(i) ;
> t16 = psi(i)/t15 ;
> t17 = pow(t16, psi(i)) ;
> t18 = -psi(i) - 1.0 ;
> t24 = lam(i)/t15 ;
> t19 = -t24 + 1.0 ;
> t20 = log(t19) ;
> t21 = lam(i)*lam(i) ;
> t22 = t15*t15 ;
> t23 = -psi(i) - 2.0 ;
> moments(i,3) = (t17*lam(i)*psi(i)*exp(t20*t18))/t15 - 3.0*t21*pow(t15,
> -2.0)*t17*t18 \
> * psi(i)*exp(t20*t23) + (t21*t23*t17*t18*lam(i)*psi(i)*exp(t20*(-psi(i) -
> 3.0)))/(t22*t15) ;
> }
> } 
> 
> if (order==2)
> {
> for(i = 0; i < n; i++)
> {
> moments(i,1) = lam(i);
> 
> t11 = lam(i) + psi(i) ;
> t12 = psi(i)/t11 ;
> t13 = pow(t12, psi(i)) ;
> t14 = -psi(i) - 1.0 ;
> t15 = -(lam(i)/t11) + 1.0 ;
> t16 = log(t15) ;
> moments(i,2) = (t13*lam(i)*psi(i)*exp(t14*t16))/t11 - pow(t11,
> -2.0)*t13*t14*(lam(i)*lam(i)) \
> *psi(i)*exp(t16*(-psi(i) - 2.0)) ;
> 
> 
> t15 = lam(i) + psi(i) ;
> t16 = psi(i)/t15 ;
> t17 = pow(t16, psi(i)) ;
> t18 = -psi(i) - 1.0 ;
> t24 = lam(i)/t15 ;
> t19 = -t24 + 1.0 ;
> t20 = log(t19) ;
> t21 = lam(i)*lam(i) ;
> t22 = t15*t15 ;
> t23 = -psi(i) - 2.0 ;
> moments(i,3) = (t17*lam(i)*psi(i)*exp(t20*t18))/t15 - 3.0*t21*pow(t15,
> -2.0)*t17*t18 \
> * psi(i)*exp(t20*t23) + (t21*t23*t17*t18*lam(i)*psi(i)*exp(t20*(-psi(i) -
> 3.0)))/(t22*t15) ;
> 
> t16 = lam(i) + psi(i) ;
> t17 = psi(i)/t16 ;
> t18 = pow(t17, psi(i)) ;
> t19 = -psi(i) - 1.0 ;
> t25 = lam(i)/t16 ;
> t20 = -t25 + 1.0 ;
> t21 = log(t20) ;
> t22 = lam(i)*lam(i) ;
> t23 = t16*t16 ;
> t24 = -psi(i) - 2.0 ;
> t26 = -psi(i) - 3.0 ;
> moments(i,4) = (t18*lam(i)*psi(i)*exp(t21*t19))/t16 - 7.0*t22*pow(t16,
> -2.0)*t18*t19 \
> *psi(i)*exp(t21*t24) +
> (6.0*t22*t24*t18*t19*lam(i)*psi(i)*exp(t21*t26))/(t23*t16) \
> -(t22*t22)*pow(t23, -2.0)*t24*t26*t18*t19*psi(i)*exp(t21*(-psi(i) - 4.0)) ;
> 
> t19 = lam(i) + psi(i) ;
> t20 = psi(i)/t19 ;
> t21 = pow(t20, psi(i)) ;
> t22 = -psi(i) - 1.0 ;
> t28 = lam(i)/t19 ;
> t23 = -t28 + 1.0 ;
> t24 = log(t23) ;
> t25 = lam(i)*lam(i) ;
> t26 = t19*t19 ;
> t27 = -psi(i) - 2.0 ;
> t29 = -psi(i) - 3.0 ;
> t30 = t25*t25 ;
> t31 = t26*t26 ;
> t32 = -psi(i) - 4.0 ;
> moments(i,5) = (t21*lam(i)*psi(i)*exp(t22*t24))/t19 -
> 15.0*t21*t22*t25*pow(t19, -2.0)\
> *psi(i)*exp(t24*t27) - 10.0*t21*t30*t22*pow(t26, -2.0)*t27*t29*psi(i) \
> *exp(t32*t24) + (25.0*t21*t22*t25*t27*lam(i)*psi(i)*exp(t24*t29))/(t26*t19) \
> + (t21*t30*t22*t32*t27*t29*lam(i)*psi(i)*exp(t24*(-psi(i) - 5.0)))/(t31*t19) ;
> }
> } 
> 
> if (order==3)
> {
> for(i = 0; i < n; i++)
> {
> moments(i,1) = lam(i);
> 
> t11 = lam(i) + psi(i) ;
> t12 = psi(i)/t11 ;
> t13 = pow(t12, psi(i)) ;
> t14 = -psi(i) - 1.0 ;
> t15 = -(lam(i)/t11) + 1.0 ;
> t16 = log(t15) ;
> moments(i,2) = (t13*lam(i)*psi(i)*exp(t14*t16))/t11 - pow(t11,
> -2.0)*t13*t14*(lam(i)*lam(i)) \
> *psi(i)*exp(t16*(-psi(i) - 2.0)) ;
> 
> 
> t15 = lam(i) + psi(i) ;
> t16 = psi(i)/t15 ;
> t17 = pow(t16, psi(i)) ;
> t18 = -psi(i) - 1.0 ;
> t24 = lam(i)/t15 ;
> t19 = -t24 + 1.0 ;
> t20 = log(t19) ;
> t21 = lam(i)*lam(i) ;
> t22 = t15*t15 ;
> t23 = -psi(i) - 2.0 ;
> moments(i,3) = (t17*lam(i)*psi(i)*exp(t20*t18))/t15 - 3.0*t21*pow(t15,
> -2.0)*t17*t18 \
> * psi(i)*exp(t20*t23) + (t21*t23*t17*t18*lam(i)*psi(i)*exp(t20*(-psi(i) -
> 3.0)))/(t22*t15) ;
> 
> t16 = lam(i) + psi(i) ;
> t17 = psi(i)/t16 ;
> t18 = pow(t17, psi(i)) ;
> t19 = -psi(i) - 1.0 ;
> t25 = lam(i)/t16 ;
> t20 = -t25 + 1.0 ;
> t21 = log(t20) ;
> t22 = lam(i)*lam(i) ;
> t23 = t16*t16 ;
> t24 = -psi(i) - 2.0 ;
> t26 = -psi(i) - 3.0 ;
> moments(i,4) = (t18*lam(i)*psi(i)*exp(t21*t19))/t16 - 7.0*t22*pow(t16,
> -2.0)*t18*t19 \
> *psi(i)*exp(t21*t24) +
> (6.0*t22*t24*t18*t19*lam(i)*psi(i)*exp(t21*t26))/(t23*t16) \
> -(t22*t22)*pow(t23, -2.0)*t24*t26*t18*t19*psi(i)*exp(t21*(-psi(i) - 4.0)) ;
> 
> t19 = lam(i) + psi(i) ;
> t20 = psi(i)/t19 ;
> t21 = pow(t20, psi(i)) ;
> t22 = -psi(i) - 1.0 ;
> t28 = lam(i)/t19 ;
> t23 = -t28 + 1.0 ;
> t24 = log(t23) ;
> t25 = lam(i)*lam(i) ;
> t26 = t19*t19 ;
> t27 = -psi(i) - 2.0 ;
> t29 = -psi(i) - 3.0 ;
> t30 = t25*t25 ;
> t31 = t26*t26 ;
> t32 = -psi(i) - 4.0 ;
> moments(i,5) = (t21*lam(i)*psi(i)*exp(t22*t24))/t19 -
> 15.0*t21*t22*t25*pow(t19, -2.0)\
> *psi(i)*exp(t24*t27) - 10.0*t21*t30*t22*pow(t26, -2.0)*t27*t29*psi(i) \
> *exp(t32*t24) + (25.0*t21*t22*t25*t27*lam(i)*psi(i)*exp(t24*t29))/(t26*t19) \
> + (t21*t30*t22*t32*t27*t29*lam(i)*psi(i)*exp(t24*(-psi(i) - 5.0)))/(t31*t19) ;
> 
> t20 = lam(i) + psi(i) ;
> t21 = psi(i)/t20 ;
> t22 = pow(t21, psi(i)) ;
> t23 = -psi(i) - 1.0 ;
> t29 = lam(i)/t20 ;
> t24 = -t29 + 1.0 ;
> t25 = log(t24) ;
> t26 = lam(i)*lam(i) ;
> t27 = t20*t20 ;
> t28 = -psi(i) - 2.0 ;
> t30 = -psi(i) - 3.0 ;
> t31 = t26*t26 ;
> t32 = t27*t27 ;
> t33 = -psi(i) - 4.0 ;
> t34 = -psi(i) - 5.0 ;
> moments(i,6) = (t22*lam(i)*psi(i)*exp(t23*t25))/t20 - 31.0*pow(t20,
> -2.0)*t22*t23*t26\
> *psi(i)*exp(t25*t28) - 65.0*t30*t22*t31*t23*pow(t27, -2.0)*t28*psi(i)*exp(t33
> \
> *t25) + (90.0*t22*t23*t26*t28*lam(i)*psi(i)*exp(t30*t25))/(t20*t27) \
> + (15.0*t30*t22*t31*t23*t33*t28*lam(i)*psi(i)*exp(t25*t34))/(t20*t32) \
> - (t30*t22*t31*t23*t33*t34*t26*t28*psi(i)*exp(t25*(-psi(i) - 6.0)))/(t32*t27)
> ;
> 
> t21 = lam(i) + psi(i) ;
> t22 = psi(i)/t21 ;
> t23 = pow(t22, psi(i)) ;
> t24 = -psi(i) - 1.0 ;
> t30 = lam(i)/t21 ;
> t25 = -t30 + 1.0 ;
> t26 = log(t25) ;
> t27 = lam(i)*lam(i) ;
> t28 = t21*t21 ;
> t29 = -psi(i) - 2.0 ;
> t31 = -psi(i) - 3.0 ;
> t32 = t27*t27 ;
> t33 = t28*t28 ;
> t34 = -psi(i) - 4.0 ;
> t35 = -psi(i) - 5.0 ;
> t36 = -psi(i) - 6.0 ;
> moments(i,7) = (t23*lam(i)*psi(i)*exp(t24*t26))/t21 - 63.0*pow(t21,
> -2.0)*t23*t24*t27 \
> *psi(i)*exp(t26*t29) - 350.0*t31*t23*t32*t24*pow(t28, -2.0)*t29*psi(i) \
> *exp(t34*t26) + (301.0*t23*t24*t27*t29*lam(i)*psi(i)*exp(t31*t26))/(t21*t28) \
> + (140.0*t31*t23*t32*t24*t34*t29*lam(i)*psi(i)*exp(t26*t35))/(t21*t33) \
> - (21.0*t31*t23*t32*t24*t34*t35*t27*t29*psi(i)*exp(t26*t36))/(t33*t28) \
> + (t31*t23*t32*t24*t34*t35*t27*t36*t29*lam(i)*psi(i)*exp(t26*(-psi(i) -
> 7.0)))/(t21*t33*t28) ;
> }
> } 
> 
> 
> if (order==4)
> {
> for(i = 0; i < n; i++)
> {
> moments(i,1) = lam(i);
> 
> t11 = lam(i) + psi(i) ;
> t12 = psi(i)/t11 ;
> t13 = pow(t12, psi(i)) ;
> t14 = -psi(i) - 1.0 ;
> t15 = -(lam(i)/t11) + 1.0 ;
> t16 = log(t15) ;
> moments(i,2) = (t13*lam(i)*psi(i)*exp(t14*t16))/t11 - pow(t11,
> -2.0)*t13*t14*(lam(i)*lam(i)) \
> *psi(i)*exp(t16*(-psi(i) - 2.0)) ;
> 
> 
> t15 = lam(i) + psi(i) ;
> t16 = psi(i)/t15 ;
> t17 = pow(t16, psi(i)) ;
> t18 = -psi(i) - 1.0 ;
> t24 = lam(i)/t15 ;
> t19 = -t24 + 1.0 ;
> t20 = log(t19) ;
> t21 = lam(i)*lam(i) ;
> t22 = t15*t15 ;
> t23 = -psi(i) - 2.0 ;
> moments(i,3) = (t17*lam(i)*psi(i)*exp(t20*t18))/t15 - 3.0*t21*pow(t15,
> -2.0)*t17*t18 \
> * psi(i)*exp(t20*t23) + (t21*t23*t17*t18*lam(i)*psi(i)*exp(t20*(-psi(i) -
> 3.0)))/(t22*t15) ;
> 
> t16 = lam(i) + psi(i) ;
> t17 = psi(i)/t16 ;
> t18 = pow(t17, psi(i)) ;
> t19 = -psi(i) - 1.0 ;
> t25 = lam(i)/t16 ;
> t20 = -t25 + 1.0 ;
> t21 = log(t20) ;
> t22 = lam(i)*lam(i) ;
> t23 = t16*t16 ;
> t24 = -psi(i) - 2.0 ;
> t26 = -psi(i) - 3.0 ;
> moments(i,4) = (t18*lam(i)*psi(i)*exp(t21*t19))/t16 - 7.0*t22*pow(t16,
> -2.0)*t18*t19 \
> *psi(i)*exp(t21*t24) +
> (6.0*t22*t24*t18*t19*lam(i)*psi(i)*exp(t21*t26))/(t23*t16) \
> -(t22*t22)*pow(t23, -2.0)*t24*t26*t18*t19*psi(i)*exp(t21*(-psi(i) - 4.0)) ;
> 
> t19 = lam(i) + psi(i) ;
> t20 = psi(i)/t19 ;
> t21 = pow(t20, psi(i)) ;
> t22 = -psi(i) - 1.0 ;
> t28 = lam(i)/t19 ;
> t23 = -t28 + 1.0 ;
> t24 = log(t23) ;
> t25 = lam(i)*lam(i) ;
> t26 = t19*t19 ;
> t27 = -psi(i) - 2.0 ;
> t29 = -psi(i) - 3.0 ;
> t30 = t25*t25 ;
> t31 = t26*t26 ;
> t32 = -psi(i) - 4.0 ;
> moments(i,5) = (t21*lam(i)*psi(i)*exp(t22*t24))/t19 -
> 15.0*t21*t22*t25*pow(t19, -2.0)\
> *psi(i)*exp(t24*t27) - 10.0*t21*t30*t22*pow(t26, -2.0)*t27*t29*psi(i) \
> *exp(t32*t24) + (25.0*t21*t22*t25*t27*lam(i)*psi(i)*exp(t24*t29))/(t26*t19) \
> + (t21*t30*t22*t32*t27*t29*lam(i)*psi(i)*exp(t24*(-psi(i) - 5.0)))/(t31*t19) ;
> 
> t20 = lam(i) + psi(i) ;
> t21 = psi(i)/t20 ;
> t22 = pow(t21, psi(i)) ;
> t23 = -psi(i) - 1.0 ;
> t29 = lam(i)/t20 ;
> t24 = -t29 + 1.0 ;
> t25 = log(t24) ;
> t26 = lam(i)*lam(i) ;
> t27 = t20*t20 ;
> t28 = -psi(i) - 2.0 ;
> t30 = -psi(i) - 3.0 ;
> t31 = t26*t26 ;
> t32 = t27*t27 ;
> t33 = -psi(i) - 4.0 ;
> t34 = -psi(i) - 5.0 ;
> moments(i,6) = (t22*lam(i)*psi(i)*exp(t23*t25))/t20 - 31.0*pow(t20,
> -2.0)*t22*t23*t26\
> *psi(i)*exp(t25*t28) - 65.0*t30*t22*t31*t23*pow(t27, -2.0)*t28*psi(i)*exp(t33
> \
> *t25) + (90.0*t22*t23*t26*t28*lam(i)*psi(i)*exp(t30*t25))/(t20*t27) \
> + (15.0*t30*t22*t31*t23*t33*t28*lam(i)*psi(i)*exp(t25*t34))/(t20*t32) \
> - (t30*t22*t31*t23*t33*t34*t26*t28*psi(i)*exp(t25*(-psi(i) - 6.0)))/(t32*t27)
> ;
> 
> t21 = lam(i) + psi(i) ;
> t22 = psi(i)/t21 ;
> t23 = pow(t22, psi(i)) ;
> t24 = -psi(i) - 1.0 ;
> t30 = lam(i)/t21 ;
> t25 = -t30 + 1.0 ;
> t26 = log(t25) ;
> t27 = lam(i)*lam(i) ;
> t28 = t21*t21 ;
> t29 = -psi(i) - 2.0 ;
> t31 = -psi(i) - 3.0 ;
> t32 = t27*t27 ;
> t33 = t28*t28 ;
> t34 = -psi(i) - 4.0 ;
> t35 = -psi(i) - 5.0 ;
> t36 = -psi(i) - 6.0 ;
> moments(i,7) = (t23*lam(i)*psi(i)*exp(t24*t26))/t21 - 63.0*pow(t21,
> -2.0)*t23*t24*t27 \
> *psi(i)*exp(t26*t29) - 350.0*t31*t23*t32*t24*pow(t28, -2.0)*t29*psi(i) \
> *exp(t34*t26) + (301.0*t23*t24*t27*t29*lam(i)*psi(i)*exp(t31*t26))/(t21*t28) \
> + (140.0*t31*t23*t32*t24*t34*t29*lam(i)*psi(i)*exp(t26*t35))/(t21*t33) \
> - (21.0*t31*t23*t32*t24*t34*t35*t27*t29*psi(i)*exp(t26*t36))/(t33*t28) \
> + (t31*t23*t32*t24*t34*t35*t27*t36*t29*lam(i)*psi(i)*exp(t26*(-psi(i) -
> 7.0)))/(t21*t33*t28) ;
> 
> t22 = lam(i) + psi(i) ;
> t23 = psi(i)/t22 ;
> t24 = pow(t23, psi(i)) ;
> t25 = -psi(i) - 1.0 ;
> t31 = lam(i)/t22 ;
> t26 = -t31 + 1.0 ;
> t27 = log(t26) ;
> t28 = lam(i)*lam(i) ;
> t29 = t22*t22 ;
> t30 = -psi(i) - 2.0 ;
> t32 = -psi(i) - 3.0 ;
> t33 = t28*t28 ;
> t34 = t29*t29 ;
> t35 = -psi(i) - 4.0 ;
> t36 = -psi(i) - 5.0 ;
> t37 = -psi(i) - 6.0 ;
> t38 = -psi(i) - 7.0 ;
> moments(i,8) = (t24*lam(i)*psi(i)*exp(t25*t27))/t22 - 127.0*pow(t22,
> -2.0)*t24*t25*t28 \
> *psi(i)*exp(t30*t27) - 1701.0*t30*t32*t24*t33*t25*pow(t29,
> -2.0)*psi(i)*exp(t35\
> *t27) + (966.0*t30*t24*t25*t28*lam(i)*psi(i)*exp(t32*t27))/(t22*t29) +
> (1050.0*t30 \
> *t32*t24*t33*t25*t35*lam(i)*psi(i)*exp(t27*t36))/(t22*t34) -
> (266.0*t30*t32*t24 \
> *t33*t25*t35*t36*t28*psi(i)*exp(t27*t37))/(t34*t29) + (28.0*t30*t32*t24*t33 \
> *t25*t35*t36*t28*t37*lam(i)*psi(i)*exp(t27*t38))/(t22*t34*t29) -
> t30*t32*t24*(t33\
> *t33)*t25*pow(t34, -2.0)*t35*t36*t37*t38*psi(i)*exp(t27*(-psi(i) - 8.0)) ;
> 
> t25 = lam(i) + psi(i) ;
> t26 = psi(i)/t25 ;
> t27 = pow(t26, psi(i)) ;
> t28 = -psi(i) - 1.0 ;
> t34 = lam(i)/t25 ;
> t29 = -t34 + 1.0 ;
> t30 = log(t29) ;
> t31 = lam(i)*lam(i) ;
> t32 = t25*t25 ;
> t33 = -psi(i) - 2.0 ;
> t35 = -psi(i) - 3.0 ;
> t36 = t31*t31 ;
> t37 = t32*t32 ;
> t38 = -psi(i) - 4.0 ;
> t39 = -psi(i) - 5.0 ;
> t40 = -psi(i) - 6.0 ;
> t41 = -psi(i) - 7.0 ;
> t42 = t36*t36 ;
> t43 = t37*t37 ;
> t44 = -psi(i) - 8.0 ;
> moments(i,9) = (t27*lam(i)*psi(i)*exp(t30*t28))/t25 - 255.0*t31*pow(t25,
> -2.0)*t27*t28 \
> *psi(i)*exp(t30*t33) - 7770.0*pow(t32,
> -2.0)*t33*t35*t27*t36*t28*psi(i)*exp(t30\
> *t38) - 36.0*t40*t41*t33*t42*t35*t27*t28*pow(t37, -2.0)*t38*t39*psi(i)*exp(t30
> \
> *t44) + (3025.0*t31*t33*t27*t28*lam(i)*psi(i)*exp(t30*t35))/(t32*t25) +
> (6951.0\
> *t33*t35*t27*t36*t28*t38*lam(i)*psi(i)*exp(t30*t39))/(t25*t37) -
> (2646.0*t31*t33\
> *t35*t27*t36*t28*t38*t39*psi(i)*exp(t30*t40))/(t32*t37) + (462.0*t31*t40*t33 \
> *t35*t27*t36*t28*t38*t39*lam(i)*psi(i)*exp(t30*t41))/(t32*t25*t37) +
> (t40*t41*t33\
> *t42*t35*t44*t27*t28*t38*t39*lam(i)*psi(i)*exp(t30*(-psi(i) - 9.0)))/(t25*t43)
> ;
> }
> } 
> 
> 
> 
> // normalization factor and shaping polynomial
> normfactor = normfactor.fill(0.0);
> poly = poly.fill(0.0);
> prediction = prediction.fill(0.0);
> for(i = 0; i <= order; i++)
> {
> for(j = 0; j <= order; j++)
> {
> normfactor = normfactor + gam(i)*gam(j)*moments.column(i+j);
> prediction = prediction + gam(i)*gam(j)*moments.column(i+j+1);
> }
> 
> for(j = 0; j < n; j++)
> {
> poly(j) = poly(j) + gam(i)*pow(y(j), (double) i);
> }
> } 
> 
> 
> // log density and prediction
> for(i = 0; i < n; i++)
> {
> poly(i) = log(DBL_EPSILON + poly(i)*poly(i));
> t1 = DBL_EPSILON + psi(i)+lam(i);
> logdensity(i) = lgamma(y(i)+psi(i)) - lgamma(psi(i)) - lgamma(y(i)+1.0) \
> + psi(i)*log(DBL_EPSILON + psi(i) / t1) + y(i)*log(DBL_EPSILON + lam(i) /
> t1); 
> logdensity(i) = poly(i) + logdensity(i) - log(DBL_EPSILON+normfactor(i));
> prediction(i) = prediction(i) / normfactor(i);
> }
> 
> f_return(0) = logdensity;
> f_return(1) = "na";
> f_return(2) = prediction;
> return octave_value_list(f_return);
> 
> }
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://www.cae.wisc.edu/mailman/listinfo/help-octave