Specification of vmeasure and emeasure Csnippets for enkf

[watermath]

I am wondering about the Csnippet specification of vmeasure and emeasure for enkf.

The documentation says the following:

"The goal of such a snippet is to fill variables named E_y with the conditional expectations of observables y. Accordingly, there should be one assignment of E_y for each observable y" and
"The goal of such a snippet is to fill variables named V_y_z with the conditional covariances of observables y, z. Accordingly, there should be one assignment of V_y_z and one assignment of V_z_y for each pair of observables y and z"

The mean of my multivariate n observables y are also parameters in my model called 'mu', so my POMP emeasure component is:

mod_emeasure<-Csnippet("
    const double *mu = &mu1;
    double E_y[n];
    int i;
    for (i=0; i<n i++) E_y[i]=mu[i];
  ")

and similarly, the variances of my n observables y are given by a model parameter Q which is diagonal (only variances, all covariances are 0), so my vmeasure component is:

mod_vmeasure<-Csnippet("
    const double *Q = &Q1;
    double V_y_z[n*n];
    double V_z_y[n*n];
    int i;
    
    for (i=0; i<n*n; i++) V_y_z[i] = 0;
    for (i=0; i<n; i++) V_y_z[(n*i)+i]=Q[i];
    for (i=0; i<n*n; i++) V_z_y[i] = V_y_z[i];
  ")

I also tried filling V_y_z and V_z_y with just the lower/upper triangle of my diagonal variance matrix as a vector.

When I run enkf, I get the error Error: in ‘enkf’: infinite or missing values in 'x'. And for eakf I get the error Error: in ‘eakf’: system is computationally singular: reciprocal condition number = 8.0338e-18.

I'm sure everything else if correctly specified as mif2 runs with no problem. So I think I must be filling the variable V_y_z incorrectly.

Is V_y_z supposed to be a vectors with elements for each observable y? Or should it have some other form?

Thank you, I appreciate your help as always. I have kept an ongoing document of my pomp lessons that I plan to add to this forum to pay forward help for others with similar questions.

[kingaa]

Just looking quickly at what you've done, I see that you've declared the E_ and V_ variables, which is a mistake. These are already defined in the context within which your snippet will be executed. Also, the y and z in the documentation you quote are variables that stand for the names of the variables in your model.
I don't have enough information to spell it out perfectly, but perhaps I can indicate the solution to you. Suppose your observables are named A and B. Then your emeasure snippet would have lines like

E_A = ...
E_B = ...

while your vmeasure snippet would have lines like

V_A_B = V_B_A = ...
V_A_A = ...
V_B_B = ...

Although you don't ask about this, I'm also wondering about your larger purpose. It seems that, since the probability distribution of the data does not depend on latent states, your model is really much simpler than a POMP. With no latent variables, there are much more efficient inference algorithms....

[watermath]

Ah I see, thank you Dr. King. I thought the E_y and V_y_z were built in variables for POMP, like when specifying lik in the dmeasure component. I wonder then how best to specify these for a general number of observables, A, B, C, ... without writing out E_A, ... and every combination of V_A_B = V_B_A = ... . This may be a trivial question but I am a novice at best when it comes to C coding.

I wanted to keep my question focused so I didn't include all the non-relevant details of the full model. My model does indeed have latent variables. Writing one POMP model has the advantage of being able to use the different algorithms implemented within POMP and the added speed of using Csnippets over other good alternatives like MARSS.

[kingaa]

@watermath : You can take the approach recommended in FAQ 3.1: extract a pointer to the beginning of the vector (for emeasure) or array (for vmeasure) and use it to index your expressions. This will be sensitive to the order in which the observables fall. You can control this, however, using the obsnames argument. You will have to make sure yourself that the order of entries in the expectation and covariance matrix corresponds to that of the entries in the data vector.