Sloping loudness

Colin Champion

Summary: this page contains the subroutine slopingbw to estimate the parameters of a sloping loudness model using the Baum-Welch algorithm. The test program illustrates its use.

doc : source : header : test : download & index

slopingbw has a header file which defines a slopingstats structure for the sufficient statistics needed by slopingbw. You may pass the statistics either through a structure as thus defined (recommended) or through several arrays. The prototypes are

double slopingbw(slopingstats s,double *mu,double **zb,double **zw,
                 double *muhat,double **zbhat,double **zwhat,int niter) ;
double slopingbw(int m,int n,double **S,double **ybar,double **ybardash,
                 int *ni,double *mu,double **zb,double **zw,
                 double *muhat,double **zbhat,double **zwhat,int niter) ;

where

s is a slopingstats structure;
m is the dimension of the feature vectors;
n is the number of sequences;
S is the matrix of scatter from sequence means;
ybar is the set of sequence means;
ybardash is the set of time-weighted sequence means;
ni is the set of sequence lengths;
mu is an initial estimate of the mean;
zb is an initial estimate of the covariance of sequence means between sequences;
zw is an initial estimate of the within-sequence covariance;
muhat will receive a revised estimate of the mean;
zbhat will receive a revised estimate of the covariance of sequence means between sequences;
zwhat will receive a revised estimate of the within-sequence covariance; and
niter is the number of iterations in the batch.

If niter>0 then that many iterations will be performed and the return value will correspond to the parameters on entry to the final iteration. If niter=0 then no reestimation will be performed: the log pdf will be computed for the supplied parameters and the ‘hat’ parameters will be ignored (set the corresponding pointers to 0 to avoid confusion). The ‘hat’ parameters may be the same arrays as the initial estimates if you want to overwrite.

For efficiency sort the sequences by length before calling slopingbw (using the sort function in slopingbw.h.

You could perform a Newton-Raphson hillclimb on the log pdf values, but – quite apart from the question of whether Baum-Welch would be faster – you should remember that the optimum mean is easily obtainable analytically from estimated variances, so the hillclimb can be performed on a reduced space.

slopingbw.h provides the following functions.

Call

  s = slopingstats(int m) ;

to construct an empty stats structure for dimension m. (s = slopingstats() ; generates a completely empty structure).

Call

  s.take(double **y,int n) ;

to accumulate information from a single group comprising n m-long observations for which the i^th dimension of the t^th observation is stored as y[t][i].

Call

  s.sort() ;

to sort the observations in increasing order of group size.

Call

  s.release() ;

to free all space allocated.

See the test program for an example of its use. Note that you need to include ij.h before slopingbw.h.

#include "memory.h"
#include <math.h>

double logcholesky(double **a,double **b,int n) ;
double    cholesky(double **a,double **b,int n) ;
static double log2pi=log(2*3.141592653589793) ;

double slopingbw(int m,int N,double **S,double **ybar,double **ybardash,
                 int *ni,double *mu,double **zb,double **zw,
                 double *muhat,double **zbhat,double **zwhat,int niter)
{ int i,j,k,t,mdash=m+1,ono,reest,dopdf,iterno,NL,NS,mlim ;
  double n0,n1,n2,ndot,*edot=vector(m),edotdot,qret,q,logalphadet ; 
  double **E=matrix(mdash,mdash),*ydash=vector(mdash),*a=vector(mdash) ; 
  double **alpha=matrix(mdash,mdash),*muzb=vector(mdash) ; 
  double **e=matrix(m,m),**zbinv=matrix(mdash,mdash),**alphainv=alpha,**Einv=E ;
  double logzwdet,logzbdet,logzbxxdet,*mux=vector(mdash) ;
  double **a0=matrix(mdash,mdash),**a1=matrix(mdash,mdash) ; 
  double **zbxxinv,**zbxxhat,*muzbxx,*muxhat ;

  for(ndot=NL=NS=t=0;t<N;t++) 
  { if(ni[t]==1) NS += 1 ; else NL += 1 ; ndot += ni[t] ; }
  if(NS>0) 
  { zbxxinv = matrix(m,m) ; zbxxhat = matrix(m,m) ;
    muzbxx = vector(m) ; muxhat = vector(m) ; 
  }

  for(iterno=0;iterno<niter||iterno==0;iterno++)
  { reest = (iterno<niter) ;
    dopdf = (iterno>=niter-1) ;
    if(iterno==0) 
    { for(i=0;i<mdash;i++) for(j=0;j<mdash;j++) zbinv[i][j] = zb[i][j] ;
      for(i=0;i<m;i++) for(j=0;j<m;j++) e[i][j] = zw[i][j] ;
    }
    else 
    { for(i=0;i<mdash;i++) for(j=0;j<mdash;j++) zbinv[i][j] = zbhat[i][j] ;
      for(i=0;i<m;i++) for(j=0;j<m;j++) e[i][j] = zwhat[i][j] ;
    }
    logzwdet = logcholesky(e,e,m) ; 
    if(NS>0) // block inversion of zb
    { logzbxxdet = logcholesky(zbinv,zbxxinv,m) ;
      for(i=0;i<m;i++) for(a[i]=j=0;j<m;j++) 
        a[i] += zbxxinv[i][j] * zbinv[j][m] ;
      for(q=zbinv[m][m],i=0;i<m;i++) q -= zbinv[i][m] * a[i] ;
      logzbdet = logzbxxdet + log(q) ;
      zbinv[m][m] = q = 1/q ;
      for(i=0;i<m;i++) for(j=0;j<m;j++) 
        zbinv[i][j] = zbxxinv[i][j] + q*a[i]*a[j] ;
      for(i=0;i<m;i++) zbinv[i][m] = zbinv[m][i] = -q*a[i] ;
    }
    else { logzbdet = logcholesky(zbinv,zbinv,mdash) ; logzbxxdet = 0 ; }

    /*----------------------------------------------------------------------- */
    /*page*/
    /*----------------------------------------------------------------------- */

    for(i=0;i<m;i++) for(edot[i]=j=0;j<m;j++) edot[i] += e[i][j] ; 
    for(edotdot=i=0;i<m;i++) edotdot += edot[i] ; 
    for(i=0;i<mdash;i++) mux[i] = mu[i] ; 

    if(dopdf) 
    { qret = NL*logzbdet + NS*logzbxxdet + ndot*(logzwdet+m*log2pi) ;
      for(i=0;i<m;i++) for(j=0;j<m;j++) qret += e[i][j] * S[i][j] ;
    }

    if(reest)
    { for(i=0;i<m;i++) for(j=0;j<m;j++) zwhat[i][j] = S[i][j] ; 
      for(i=0;i<mdash;i++) for(j=0;j<mdash;j++) zbhat[i][j] = 0 ; 
      for(i=0;i<mdash;i++) for(muzb[i]=j=0;j<mdash;j++)
        muzb[i] += mux[j] * zbinv[i][j] ; 
      for(i=0;i<mdash;i++) muhat[i] = 0 ; // may overwrite mu - beware!
      if(NS) 
      { for(i=0;i<m;i++) for(muzbxx[i]=j=0;j<m;j++)
          muzbxx[i] += mux[j] * zbxxinv[i][j] ; 
        for(i=0;i<m;i++) for(j=0;j<m;j++) zbxxhat[i][j] = 0 ; 
        for(i=0;i<m;i++) muxhat[i] = 0 ; 
      }
    }

    for(ono=-1,t=0;t<N;ono=ni[t],t++) if(ni[t]>0)
    { if(ono!=ni[t]&&ni[t]==1) 
      { n0 = 1 ; 
        for(i=0;i<m;i++) for(j=0;j<m;j++) 
          alphainv[i][j] = e[i][j] + zbxxinv[i][j] ;
        logalphadet = -logcholesky(alphainv,alpha,m) ; 
        if(dopdf)
        { for(i=0;i<m;i++) for(k=0;k<m;k++) 
          { for(q=j=0;j<m;j++) q += e[i][j] * alpha[j][k] ; a0[i][k] = q ; }
          for(i=0;i<m;i++) for(k=0;k<m;k++) 
          { for(q=j=0;j<m;j++) q += a0[i][j] * zbxxinv[j][k] ; 
            a1[i][k] = q ; 
          }
        }
      }
      if(ono!=ni[t]&&ni[t]>1) // terms whose cost is cubic in m are shared ...
      { n0 = ni[t] ;          // ... between all groups of the same length
        n1 = n0*(n0-1)/2 ; 
        n2 = n1*(2*n0-1)/3 ;
        // compute E
        for(i=0;i<m;i++) for(j=0;j<m;j++) E[i][j] = n0 * e[i][j] ; 
        for(i=0;i<m;i++) E[i][m] = E[m][i] = n1 * edot[i] ; 
        E[m][m] = n2 * edotdot ; 
        // compute alpha
        for(i=0;i<mdash;i++) for(j=0;j<mdash;j++) 
          alphainv[i][j] = E[i][j] + zbinv[i][j] ; 
        logalphadet = -logcholesky(alphainv,alpha,mdash) ; 
        if(dopdf)
        { for(i=0;i<mdash;i++) for(k=0;k<mdash;k++) 
          { for(q=j=0;j<mdash;j++) q += E[i][j] * alpha[j][k] ; a0[i][k] = q ; }
          for(i=0;i<mdash;i++) for(k=0;k<mdash;k++) 
          { for(q=j=0;j<mdash;j++) q += a0[i][j] * zbinv[j][k] ; 
            a1[i][k] = q ; 
          }
          cholesky(E,Einv,mdash) ; 
        }
      }
      /*--------------------------------------------------------------------- */
      /*page*/
      /*--------------------------------------------------------------------- */

      // compute ydash
      for(i=0;i<m;i++) for(ydash[i]=j=0;j<m;j++) 
        ydash[i] += n0 * e[i][j] * ybar[t][j] ;
      if(n0>1) for(ydash[m]=i=0;i<m;i++) 
        ydash[m] += n1 * edot[i] * ybardash[t][i] ; 

      // compute pdf
      if(dopdf) 
      { qret -= logalphadet  ;
        if(n0==1)
        { for(i=0;i<m;i++) a[i] = mux[i] - ybar[t][i] ;
          for(i=0;i<m;i++) for(j=0;j<m;j++) qret += a[i] * a1[i][j] * a[j] ; 
        }
        else
        { for(i=0;i<mdash;i++) for(a[i]=mux[i],j=0;j<mdash;j++) 
            a[i] -= ydash[j] * Einv[i][j] ;
          for(i=0;i<mdash;i++) for(j=0;j<mdash;j++) 
            qret += a[i] * a1[i][j] * a[j] ; 
          for(q=i=0;i<m;i++) q += (ybar[t][i]-ybardash[t][i]) * edot[i] ;
          qret -= 6*n1*q*q / ((n0+1)*edotdot) ;
        }
      }
      if(reest==0) continue ;

      // compute a
      if(n0==1) for(i=0;i<m;i++) for(a[i]=j=0;j<m;j++)
        a[i] += (ydash[j]+muzbxx[j]) * alpha[i][j] ; 
      else for(i=0;i<mdash;i++) for(a[i]=j=0;j<mdash;j++)
        a[i] += (ydash[j]+muzb[j]) * alpha[i][j] ; 

      // update zwhat (first step, for all sequence lengths)
      for(i=0;i<m;i++) for(j=0;j<m;j++) zwhat[i][j] += 
        n0 * ( (a[i]-ybar[t][i])*(a[j]-ybar[t][j]) + alpha[i][j] ) ;

      // reestimation: update muhat, zbhat
      if(n0==1) 
      { for(i=0;i<m;i++) muxhat[i] += a[i] ; 
        for(i=0;i<m;i++) for(j=0;j<m;j++) 
          zbxxhat[i][j] += a[i]*a[j] + alpha[i][j] ; 
      }
      else 
      { for(i=0;i<mdash;i++) muhat[i] += a[i] ; 
        for(i=0;i<mdash;i++) for(j=0;j<mdash;j++) 
          zbhat[i][j] += a[i]*a[j] + alpha[i][j] ; 
        for(i=0;i<m;i++) a[i] = a[m] * ( a[i]-ybardash[t][i] ) ;
        // update zwhat (second step)
        for(i=0;i<m;i++) for(j=0;j<m;j++) 
          zwhat[i][j] += n1 * ( a[i] + a[j] + alpha[i][m] + alpha[j][m] ) +
                         n2 * ( a[m]*a[m] + alpha[m][m] ) ;
      }
    }
    /*----------------------------------------------------------------------- */
    /*page*/
    /*----------------------------------------------------------------------- */

    if(reest==0) continue ;

    // put new estimates together
    for(i=0;i<m;i++) for(j=0;j<m;j++) zwhat[i][j] /= ndot ; 
    for(i=0;i<mdash;i++) muhat[i] /= NL ; 
    for(i=0;i<mdash;i++) for(j=0;j<mdash;j++) 
      zbhat[i][j] -= NL*muhat[i]*muhat[j] ;

    if(NS==0) 
    { for(i=0;i<mdash;i++) for(j=0;j<mdash;j++) zbhat[i][j] /= N ; continue ; }
    
    for(i=0;i<m;i++) { muxhat[i] /= NS ; a[i] = muxhat[i] - muhat[i] ; }
    for(i=0;i<m;i++) for(j=0;j<m;j++) 
      zbxxhat[i][j] += NL*NS*a[i]*a[j]/N - NS*muxhat[i]*muxhat[j] ;
    // zbhat is V, zbxxhat is W

    cholesky(zbhat,zbinv,m) ;             // compute zbinv as 1/Vxx
    for(i=0;i<m;i++) for(j=0;j<m;j++)     // new zetaBhat|xx
      zbhat[i][j] = (zbhat[i][j]+zbxxhat[i][j]) / N ;

    for(i=0;i<m;i++) for(a[i]=j=0;j<m;j++) 
      a[i] += zbhat[j][m] * zbinv[i][j] ; // Vsx*(1/Vxx)
    for(i=0;i<m;i++)                      // (Vsx/Vxx)*zetaBxxhat
    { for(q=j=0;j<m;j++) q += a[j] * zbhat[i][j] ; zbhat[i][m] = q ; }

    for(q=i=0;i<m;i++) for(j=0;j<m;j++) q += a[i] * zbxxhat[i][j] * a[j] ;
    zbhat[m][m] = zbhat[m][m]/NL + q/N ;
    for(q=i=0;i<m;i++) q += a[i] * zbhat[m][i] ; //zbhat[m][] is Vsx
    zbhat[m][m] -= NS*q/(N*(double)NL) ;
    for(i=0;i<m;i++) zbhat[m][i] = zbhat[i][m] ;

    // zbhat now contains new estimate of zetab; sort out the mean
    cholesky(zbhat,zbinv,m) ;             // compute zbinv as 1/zetaBxxhat
    for(i=0;i<m;i++) muhat[i] = (NL*muhat[i]+NS*muxhat[i]) / N ;
    for(i=0;i<m;i++) for(j=0;j<m;j++) 
      muhat[m] -= NS * (muhat[i]-muxhat[i]) * zbinv[i][j] * zbhat[j][m] / NL ;
  }

  // free workspace
  free(edot,ydash,a,muzb,mux) ; freematrix(e,zbinv,E,alpha,a0,a1) ; 
  if(NS>0) { freematrix(zbxxinv,zbxxhat) ; free(muzbxx,muxhat) ; }
  return -qret / (2*ndot) ; 
}
/*--------------------------------------------------------------------------- */

• slopingstats • take • release • sort • slopingbw

double slopingbw(int m,int n,double **S,double **ybar,double **ybardash,
                 int *ni,double *mu,double **zb,double **zw,
                 double *muhat,double **zbhat,double **zwhat,int niter) ;

struct slopingstats
{ int m,N,len,*n ; double **S,**ybar,**ybardash ;
  slopingstats() { N = len = m = 0 ; n = 0 ; S = ybar = ybardash = 0 ; }
  slopingstats(int p) 
  { this[0] = slopingstats() ; m = p ; S = matrix(m+1,m+1) ; }
  void take(double **y,int ni)
  { if(ni<=0) return ;
    double q=1.0/ni ;
    int t,j,k ;
    if(len==N)
    { len += 20 + len ;
      n = ivector(n,len) ;
      ybar = matrix(ybar,len,m) ;
      ybardash = matrix(ybardash,len,m) ;
      for(t=N;t<len;t++)
      { n[t] = 0 ; for(j=0;j<m;j++) ybar[t][j] = ybardash[t][j] = 0 ; }
    }
    for(t=0;t<ni;t++) for(j=0;j<m;j++) ybar[N][j] += q * y[t][j] ;
    if(ni>1) for(q=2/(ni*(ni-1.0)),t=0;t<ni;t++) for(j=0;j<m;j++) 
      ybardash[N][j] += q * t * y[t][j] ;
    for(j=0;j<m;j++) for(k=0;k<m;k++) for(t=0;t<ni;t++)
      S[j][k] += (y[t][j]-ybar[N][j]) * (y[t][k]-ybar[N][k]) ;
    n[N++] = ni ;
  }
  void release() 
  { free(n) ; freematrix(S,ybar,ybardash) ; this[0] = slopingstats() ; }
  void sort()
  { int t,j ;
    double **x=matrix(N,m) ;
    ij *list = ijvector(n,N) ;
    ijsort(list,N) ;
    for(t=0;t<N;t++) for(j=0;j<m;j++) x[t][j] = ybar[list[t].j][j] ;
    for(t=0;t<N;t++) for(j=0;j<m;j++) ybar[t][j] = ybardash[list[t].j][j] ;
    for(t=0;t<N;t++) n[t] = list[t].i ;
    freematrix(ybardash) ;
    ybardash = ybar ; 
    ybar = x ;
    free(list) ;
  }
} ;
double slopingbw(slopingstats s,double *mu,double **zb,double **zw,
                 double *muhat,double **zbhat,double **zwhat,int niter) 
{ return slopingbw(s.m,s.N,s.S,s.ybar,s.ybardash,s.n,mu,zb,zw,
                   muhat,zbhat,zwhat,niter) ; }

doc : source

slopingtest generates random data under a sloping loudness model and estimates parameters using slopingbw in two different ways: firstly it hillclimbs on the log pdf, the it performs the Baum-Welch algorithm until convergence. The data sizes wired in are small and many of the sequences are singletons, making the test quite demanding. The hillclimb returns a log pdf much greater than is assigned to the true parameters, and the Baum-Welch procedure (in many fewer iterations) does almost as well.

You’re meant to just run it without parameters, but if you insist the calling sequence is

slopingtest m seed

where m is the dimension of the space and seed seeds the random number generator.

Run without parameters you get the following output:

20 groups with avge size = 2.0
3235 function evals
logpdf=3.301902 (10 iterations)
25 Baum-Welch iterations
pdf       : 2.6481
pdf  (ML) : 3.3769
pdf  (BW) : 3.3754

mu:
  -0.203  -1.009   1.156   0.454  -0.491   5.813

mu MLE:
  -0.224  -0.941   1.068   0.393  -0.461  16.136

mu Baum-Welch estimator:
  -0.224  -0.941   1.068   0.393  -0.460  16.116

-----------------------------------------

zeta[W]*1000:
   7.412   0.469  -1.569   1.433   0.118
   0.469   5.678   0.251   0.149   1.220
  -1.569   0.251   6.750   1.445  -0.538
   1.433   0.149   1.445   6.758  -0.704
   0.118   1.220  -0.538  -0.704   5.626

zeta[W] MLE:
   7.840   1.627  -0.485   3.835  -2.312
   1.627   3.998  -0.387  -1.419  -0.441
  -0.485  -0.387   6.988   1.601   1.038
   3.835  -1.419   1.601   9.948   0.004
  -2.312  -0.441   1.038   0.004   2.775

zeta[W] Baum-Welch:
   7.706   1.637  -0.578   3.696  -2.253
   1.637   4.012  -0.510  -1.394  -0.432
  -0.578  -0.510   6.701   1.521   0.885
   3.696  -1.394   1.521   9.789   0.077
  -2.253  -0.432   0.885   0.077   2.756

-----------------------------------------

zeta[B]*10:
   0.250   0.030  -0.163  -0.093   0.042   0.075
   0.030   0.524  -0.369  -0.204   0.309   0.345
  -0.163  -0.369   0.764   0.362  -0.336  -0.459
  -0.093  -0.204   0.362   0.370  -0.255  -0.280
   0.042   0.309  -0.336  -0.255   0.334   0.301
   0.075   0.345  -0.459  -0.280   0.301   0.413

zeta[B] MLE:
   0.103   0.088  -0.082  -0.080   0.079   0.154
   0.088   0.312  -0.123   0.018   0.117   0.203
  -0.082  -0.123   0.648   0.461  -0.205  -0.403
  -0.080   0.018   0.461   0.532  -0.140  -0.378
   0.079   0.117  -0.205  -0.140   0.146   0.251
   0.154   0.203  -0.403  -0.378   0.251   0.516

zeta[B] Baum-Welch:
   0.102   0.087  -0.082  -0.080   0.077   0.154
   0.087   0.311  -0.123   0.016   0.117   0.204
  -0.082  -0.123   0.649   0.458  -0.203  -0.402
  -0.080   0.016   0.458   0.531  -0.143  -0.375
   0.077   0.117  -0.203  -0.143   0.147   0.250
   0.154   0.204  -0.402  -0.375   0.250   0.514

doc : source

• ranvec • vmatrix • vvector • f • main

#include "memory.h"
#include <math.h>
#include "random.h"
#include "ij.h"
#include "slopingbw.h"

double uoamax(double (*func)(double*),double *x,int n,
              double xinc,double tol,int maxiter) ;
int uoacond() ;
double cholsqrt(double **a,double **b,int n) ;

static slopingstats stats ;
static int ncall ;

/*--------------------------------------------------------------------------- */

static void ranvec(double *x,int m,double *sigma,double **rot,double *mu)
{ int i,j ; 
  double q ;
  for(i=0;i<m;i++) x[i] = mu[i] ; 
  for(i=0;i<m;i++) 
  { q = sigma[i] * gaussv() ; for(j=0;j<m;j++) x[j] += q * rot[i][j] ; }
}
/*--------------------------------------------------------------------------- */

static int vmatrix(double *v,double **z0,int m) 
{ int i,j,k,t ; 
  double q,**z=matrix(m,m) ;

  for(t=i=0;i<m;i++) for(j=0;j<=i;j++,t++) z[i][j] = v[t] ; 
  for(i=0;i<m;i++) for(k=0;k<=i;k++) 
  { for(q=j=0;j<=k;j++) q += z[i][j] * z[k][j] ; z0[i][k] = z0[k][i] = q ; }
  freematrix(z) ; 
  return t ;
}
static void vvector(double *mu,double **zb,double **zw,int m,double *v)
{ int i,j,t,mdash=m+1 ;
  double **z=matrix(m,m) ;
  for(t=0;t<mdash;t++) v[t] = mu[t] ;
  cholsqrt(zw,z,m) ; 
  for(j=0;j<m;j++) for(i=0;i<=j;i++) v[t++] = z[i][j] ; 
  freematrix(z) ;
  z = matrix(mdash,mdash) ;
  cholsqrt(zb,z,mdash) ; 
  for(j=0;j<mdash;j++) for(i=0;i<=j;i++) v[t++] = z[i][j] ; 
}
/*--------------------------------------------------------------------------- */
/*page*/
/*--------------------------------------------------------------------------- */

static double f(double *v)
{ int i,t,m=stats.m,mdash=m+1 ;
  double q,*mu=vector(mdash),**zw=matrix(m,m),**zb=matrix(mdash,mdash) ;
  for(t=0;t<mdash;t++) mu[t] = v[t] ;
  t += vmatrix(v+t,zw,m) ;
  t += vmatrix(v+t,zb,mdash) ;
  q = slopingbw(stats,mu,zb,zw,0,0,0,0) ;
  free(mu) ;
  freematrix(zw,zb) ;
  ncall += 1 ;
  return q ;
}
/*--------------------------------------------------------------------------- */

int main(int argc,char **argv)
{ if(argc>2) ranset(atoi(argv[2])) ; 
  
  int m=argc>1?atoi(argv[1]):5,mdash=m+1,ndash=20 ;
  int groupno,i,j,k,t,nneg,iter,targ,ndot,cond,*ni=ivector(ndash) ; ; 
  int vlen = mdash + (mdash*(mdash+1))/2 + (m*(m+1))/2 ; 
  double ***y=(double ***) cjcalloc(ndash,sizeof(double**)) ;
  double q,q0,*eb=vector(mdash),*ew=vector(m),**xs=matrix(ndash,mdash),scale ;
  double **rotw=matrix(m,m),**rotb=matrix(mdash,mdash),*mu=vector(mdash) ;
  double *mubw=vector(mdash),**zbbw=matrix(mdash,mdash),**zwbw=matrix(m,m) ;
  double *muml=vector(mdash),**zbml=matrix(mdash,mdash),**zwml=matrix(m,m) ;
  double *v=vector(vlen),**zw=matrix(m,m),**zb=matrix(mdash,mdash) ;
  double *a=vector(m),qpdf,qmlpdf,q0pdf,qbwpdf ;

  stats = slopingstats(m) ;

  // generate group sizes
  for(ndot=t=0;t<ndash/2;t++) ndot += ( ni[t] = 1 ) ; 
  for(;t<ndash;t++) ndot += ( ni[t] = 3 ) ; 
  isort(ni,ndash) ; 
  scale = ndot / ndash ; // avge group size
  printf("%d groups with avge size = %.1f\n",ndash,scale) ;

  // generate random global parameters
  ranrot(rotw,m) ;
  for(i=0;i<m;i++) ew[i] = 1/(20*sqrt(ranf())) ; 
  for(i=0;i<m;i++) for(k=0;k<m;k++) for(zw[i][k]=j=0;j<m;j++)
    zw[i][k] += rotw[j][i] * rotw[j][k] * ew[j] * ew[j] ; 

  for(i=0;i<mdash;i++) mu[i] = gaussv() ; 
  ranrot(rotb,mdash) ; 
  for(i=0;i<mdash;i++) eb[i] = 1/(20*sqrt(ranf()/scale)) ; 
  for(i=0;i<mdash;i++) for(k=0;k<mdash;k++) for(zb[i][k]=j=0;j<mdash;j++)
    zb[i][k] += rotb[j][i] * rotb[j][k] * eb[j] * eb[j] ; 

  /*------------------------------------------------------------------------- */
  /*page*/
  /*------------------------------------------------------------------------- */

  // generate group (x,s) vectors and observations
  for(groupno=0;groupno<ndash;groupno++) 
  { ranvec(xs[groupno],mdash,eb,rotb,mu) ; 
    xs[groupno][m] /= scale ; 
    y[groupno] = matrix(ni[groupno],m) ; 
    for(t=0;t<ni[groupno];t++) 
    { for(i=0;i<m;i++) a[i] = xs[groupno][i] + t*xs[groupno][m] ;
      ranvec(y[groupno][t],m,ew,rotw,a) ; 
    } 
    stats.take(y[groupno],ni[groupno]) ;
  }
  free(ni) ;
  stats.sort() ;

  mu[m] /= scale ;
  for(i=0;i<mdash;i++) { zb[i][m] /= scale ; zb[m][i] /= scale ; }
  qpdf = slopingbw(stats,mu,zb,zw,0,0,0,0) ;

  /*----------------------------------- ML ---------------------------------- */

  for(i=0;i<m;i++) for(j=0;j<m;j++) zwml[i][j] = 0 ;
  for(i=0;i<m;i++) { muml[i] = 0 ; zwml[i][i] = 1 ; }
  for(i=0;i<mdash;i++) for(j=0;j<mdash;j++) zbml[i][j] = 0 ;
  for(i=0;i<mdash;i++) zbml[i][i] = 1 ;
  vvector(muml,zbml,zwml,m,v) ;
  ncall = 0 ;

  qmlpdf = uoamax(f,v,vlen,0.1,1e-4,500000) ;
  cond = uoacond() ;
  if(cond!=0) printf("*** cond=%d ***: ",cond) ;
  printf("%d function evals\n",ncall) ;

  for(t=0;t<mdash;t++) muml[t] = v[t] ;
  t += vmatrix(v+t,zwml,m) ;
  t += vmatrix(v+t,zbml,mdash) ;
  
  /*------------------------------- Baum Welch ------------------------------ */

  for(i=0;i<mdash;i++) for(j=0;j<mdash;j++) zbbw[i][j] = 0 ;
  for(i=0;i<mdash;i++) { mubw[i] = 0 ; zbbw[i][i] = 1 ; }
  for(i=0;i<m;i++) for(j=0;j<m;j++) zwbw[i][j] = 0 ;
  for(i=0;i<m;i++) zwbw[i][i] = 1 ;

  for(targ=10,iter=5;iter<1000;q=qbwpdf,iter+=5)
  { qbwpdf = slopingbw(stats,mubw,zbbw,zwbw,mubw,zbbw,zwbw,5) ;
    if(iter==targ)
    { printf("logpdf=%.6f (%d iterations)\n",qbwpdf,iter) ; targ *= 10 ; } 
    if(iter==5||qbwpdf-q>q0) q0 = qbwpdf - q ;
    if(iter>1&&qbwpdf-q<q0/1000) break ;
  } 
  printf("%d Baum-Welch iterations\n",iter) ;
  stats.release() ;

  /*---------------------------- print estimates ---------------------------- */

  printf("pdf       : %.4f\n",qpdf) ;
  printf("pdf  (ML) : %.4f\n",qmlpdf) ;
  printf("pdf  (BW) : %.4f\n\n",qbwpdf) ;

  printf("mu:\n") ;
  for(i=0;i<m;i++) printf("%8.3f",mu[i]) ; 
  printf("%8.3f\n\n",100*mu[m]) ;

  printf("mu MLE:\n") ;
  for(i=0;i<m;i++) printf("%8.3f",muml[i]) ; 
  printf("%8.3f\n\n",100*muml[m]) ; 

  printf("mu Baum-Welch estimator:\n") ;
  for(i=0;i<m;i++) printf("%8.3f",mubw[i]) ; 
  printf("%8.3f\n",100*mubw[m]) ; 
  printf("\n-----------------------------------------\n\n") ; 

  printf("zeta[W]*1000:\n") ;
  for(i=0;i<m;i++)
  { for(j=0;j<m;j++) printf("%8.3f",1000*zw[i][j]) ; printf("\n") ; }
  printf("\n") ;

  printf("zeta[W] MLE:\n") ;
  for(i=0;i<m;i++) 
  { for(j=0;j<m;j++) printf("%8.3f",1000*zwml[i][j]) ; printf("\n") ; }
  printf("\n") ; 

  printf("zeta[W] Baum-Welch:\n") ;
  for(i=0;i<m;i++) 
  { for(j=0;j<m;j++) printf("%8.3f",1000*zwbw[i][j]) ; printf("\n") ; }
  printf("\n-----------------------------------------\n\n") ; 

  printf("zeta[B]*10:\n") ;
  for(i=0;i<mdash;i++)
  { for(j=0;j<mdash;j++) printf("%8.3f",10*zb[i][j]) ; printf("\n") ; }
  printf("\n") ;

  printf("zeta[B] MLE:\n") ;
  for(i=0;i<mdash;i++) 
  { for(j=0;j<mdash;j++) printf("%8.3f",10*zbml[i][j]) ; printf("\n") ; }
  printf("\n") ; 

  printf("zeta[B] Baum-Welch:\n") ;
  for(i=0;i<mdash;i++) 
  { for(j=0;j<mdash;j++) printf("%8.3f",10*zbbw[i][j]) ; printf("\n") ; }
  printf("\n") ; 
}

doc : source

This download procedure works under Firefox (so long as you enable popup windows for the site) but probably not under any other browser.

slopingbw.c slopingbw.h slopingtest.c       # this file
utils: memory.h ij.h random.h ranf.c cholesky.c qr.c ranrot.c
optim: newuoa.c
---
g++ -o slopingtest -O2 -g (*.c) ; slopingtest