Skip to content
New issue

Have a question about this project? # for a free GitHub account to open an issue and contact its maintainers and the community.

#

By clicking “#”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? # to your account

Jump to bottom

Add more Gaussian process kernels #234

Closed
paul-buerkner opened this issue Jul 7, 2017 · 23 comments · Fixed by #1688
Closed

Add more Gaussian process kernels #234

paul-buerkner opened this issue Jul 7, 2017 · 23 comments · Fixed by #1688
Labels
feature formula good first issue
Milestone
brms 2.22.0

Comments

@paul-buerkner
Copy link
Owner

paul-buerkner commented Jul 7, 2017 •
edited
Loading

Right now, only the exponentiated-quadratic kernel is supported as it has a native implementation in Stan. However, there seem to be quite a few other kernels worth considering. This issue is ment to provide a list of kernels with results concerning their feasibility in Stan.
@jae0
Copy link

jae0 commented Nov 1, 2017 •
edited
Loading

Something like the following would cover the main forms discussed in :

[https://en.wikipedia.org/wiki/Mat%C3%A9rn_covariance_function]


functions {
            matrix matern_covariance( int N, matrix dist, real phi, real sigma_sq, int COVFN) {
              matrix[N,N] S;
              real dist_phi; 
              real sqrt3;
              real sqrt5;
              sqrt3=sqrt(3.0);
              sqrt5=sqrt(5.0);

              if (COVFN==1) { // exponential == Matern nu=1/2 , (p=0; nu=p+1/2)
                for(i in 1:(N-1)){
                  for(j in (i+1):N){
                    dist_phi = fabs(dist[i,j])/phi;
                    S[i,j] = sigma_sq * exp(- dist_phi ); 
                }}

              
              } else if (COVFN==2) { // Matern nu= 3/2 covariance
                for(i in 1:(N-1)){
                  for(j in (i+1):N){
                   dist_phi = fabs(dist[i,j])/phi;
                   S[i,j] = sigma_sq * (1 + sqrt3 * dist_phi) * exp(-sqrt3 * dist_phi);
                }}
              

              } else if (COVFN==3) { // Matern nu=5/2 covariance
                for(i in 1:(N-1)){
                  for(j in (i+1):N){
                    dist_phi = fabs(dist[i,j])/phi;
                    S[i,j] = sigma_sq * (1 + sqrt5 *dist_phi + 5* pow(dist_phi,2)/3) * exp(-sqrt5 *dist_phi);
                }}
            
              } else if (COVFN==4) { // Matern as nu->Inf become Gaussian (aka squared exponential cov)
                for(i in 1:(N-1)){
                  for(j in (i+1):N){
                    dist_phi = fabs(dist[i,j])/phi;
                    S[i,j] = sigma_sq * exp( -pow(dist_phi,2)/2 ) ;
                }}
              } 

                 
              for(i in 1:(N-1)){
              for(j in (i+1):N){
                S[j,i] = S[i,j];  // fill upper triangle
              }}

              // create diagonal: nugget(nonspatial) + spatial variance +  eps ensures positive definiteness
              for(i in 1:N) S[i,i] = sigma_sq ;            
              return(S)   ;
            }

      }

@paul-buerkner
Copy link
Owner Author

Thanks, this is very helpful!

@gabriuma
Copy link

gabriuma commented Dec 4, 2018

To use the Mattern covfun in approximate GP, we just to change the spectral density function (function 'spd' in the Stan code):

   #Spectral density for Mattern nu=3/2
real spd(real alpha, real rho, real w) {
	real S;
	S = 4*alpha^2 * (sqrt(3)/rho)^3 * 1/((sqrt(3)/rho)^2 + w^2)^2;
	return S;
}

    #Spectral density for Mattern nu=1/2
real spd(real alpha, real rho, real w) {
	real S;
	S = 2*alpha^2 * (1/rho) * 1/((1/rho)^2 + w^2);
	return S;
}

    #Spectral density for Mattern nu=5/2
real spd(real alpha, real rho, real w) {
	real S;
	S = 16/3*alpha^2 * (sqrt(5)/rho)^5 * 1/((sqrt(5)/rho)^2 + w^2)^3;
	return S;
}

All the rest Stan code remains the same.

@benearnthof
Copy link

Has there been an update on this? I'm new to this package and would like to use a Matern kernel in a model, but I'm not sure where to put the Stan code snippets.

@paul-buerkner
Copy link
Owner Author

I think(?) that matern kernels are not supported in Stan itself so supporting them in brms should not be a major issue. Perhaps I find the time in the upcoming weeks.

@benearnthof
Copy link

For future reference: It seems like the Matern covariance functions are available in The Stan Math library

@paul-buerkner paul-buerkner mentioned this issue May 20, 2020
Closed
@ASKurz ASKurz mentioned this issue Sep 11, 2020
Open
@jgradym
Copy link

jgradym commented Nov 9, 2020

Just wanted to voice support for adding a Ornstein–Uhlenbeck option for running a brm() function with a phylogeny. This seems to be the new standard evolutionary covariance structure when dealing with phylogenies, as its more realistic than others (ie, consistent with stabilizing selection). An option of model = "OU", "Brownian" etc would be super helpful. As is, I can't use brms for most phylo analyses. Most recent evol packages (using MLE instead of Bayesian) have about 5 standard options (eg see ?phylopars in Rphylopars). But the most important, by far is the OU one. Would be fantastic to get it into brms!

@potash
Copy link

potash commented Jun 12, 2024

Hello, I would like to see the Matern covariance functions (specifically gp_exponential_cov) in brms. Can I help make this happen? I have been browsing the source code and it seems fairly straightforward. Some questions I have:

  1. Do you have any documentation to help new contributors? I have never developed an R package. How do I install it from local source? Do I need to recompile/reinstall every time I modify the source?

  2. will the brms gp(x, sdgp, lscale, zgp) function in stan take cov as a parameter? Does stan have string parameters or will there be a set of constants int GP_EXP_QUAD_COV=1; int GP_EXPONENTIAL_COV=2; etc.?

  3. How can .predictor_gp() get passed the cov parameter?

  4. Is it necessary to implement the approximate gp for these covariances as well? I feel this should be done eventually but as I and many others do not necessarily need it, it would be OK and faster to start with only exact versions of these kernels and throw an error for approximate.

Thanks!

@potash
Copy link

potash commented Jun 12, 2024 •
edited
Loading

Edit: I now see @gabriuma's comment about approximate Matern covariances and see that it should be quite easy so perhaps my third question is moot.

@paul-buerkner
Copy link
Owner Author

I will add this issue to the list of issues for the next release to make sure I will at least implement the Matern kernels. I will have to think about the best way to implement different kernels in the Stan code so it is a bit hard for me to suggest concrete steps that you could already take to implement this.

@paul-buerkner paul-buerkner added this to the brms 2.22.0 milestone Jun 13, 2024
@paul-buerkner paul-buerkner mentioned this issue Sep 17, 2024
Merged
@paul-buerkner
Copy link
Owner Author

This has been implemented in #1688. I added matern and exponential kernels including HSGP versions of the matern kernels. If people need other kernels at some point, I suggest to open new issues for these specific kernels.

@potash
Copy link

potash commented Sep 23, 2024

Fantastic, thank you! I will install and test it soon, but just a quick note/question: the exponential kernel is a matern kernel (smoothness 1/2), so I was hoping an HSGP approximation is supported but from the way your last comment is worded it seems maybe not?

@paul-buerkner
Copy link
Owner Author

Good point. :-D

I didn't see the math of the HSGP for Matern 1/2 in our HSGP paper and I didn't have time to do it myself. If you have time to compute the spectral density for this kernel (if it exists) then it would be easy to implement it.

@avehtari was there a reason we didn't discuss the Matern 1/2 kernel in our paper?

@gabriuma
Copy link

gabriuma commented Sep 23, 2024 via email

@paul-buerkner
Copy link
Owner Author

Perfect! Thank you!

@avehtari
Copy link
Contributor

@avehtari was there a reason we didn't discuss the Matern 1/2 kernel in our paper?

Trying to get the paper out in finite time?

@paul-buerkner
Copy link
Owner Author

Fair point. :-D Thanks!

@gabriuma
Copy link

gabriuma commented Sep 23, 2024 via email

@potash
Copy link

potash commented Sep 23, 2024

@gabriuma Matern 1/2 is fairly common in soils and agriculture, maybe it's just a historical accident that should be corrected... or maybe it makes sense because these variables can be very rough with a lot of variation (not measurement error) on very short length scales. Maybe they would be better represented as a sum of two smoother kernels but matern 1/2 could be a decent choice for simplicity.

paul-buerkner added a commit that referenced this issue Sep 24, 2024
@paul-buerkner
enable HSGPs for the exponential kernel #234
7838419
@paul-buerkner
Copy link
Owner Author

paul-buerkner commented Sep 24, 2024 •
edited
Loading

HSGPs for the exponential kernel are now supported (via #1688)

@avehtari
Copy link
Contributor

How about HSGP with periodic kernel? I was going to test the Birthdays case study with brms + Pathfinder, but would need periodic HSGP

@paul-buerkner
Copy link
Owner Author

paul-buerkner commented Sep 28, 2024 via email

@gabriuma
Copy link

gabriuma commented Sep 29, 2024 via email

@paul-buerkner paul-buerkner mentioned this issue Sep 30, 2024
Open
# for free to join this conversation on GitHub. Already have an account? # to comment
Assignees
No one assigned
Labels
feature formula good first issue
Projects
None yet
Milestone
brms 2.22.0
Development

Successfully merging a pull request may close this issue.

Implement more GP kernels paul-buerkner/brms
7 participants
@potash @jae0 @avehtari @paul-buerkner @jgradym @gabriuma @benearnthof