My question is simple. How we can add a function created via symbolic derivative from SymbolicC++ to the stan program in a mathematica interface? Can i load this library indirectly in Stan?
Thanks!
Stan isn’t able to use external packages, unfortunately. You would need to specify your function using Stan syntax in the functions{} block of your Stan model
Thats very bad. Then how can i perform a derivative of a function in stan? Cause i need to use a derivative of a complex function to fit my model
Another solution is to perform the symbolic derivative and just copy explicit the result in another archive. Can i at least load an external c++ archive to the stan interface in mathematica?
No that’s not possible. Stan uses custom data types and functions within its c++, so you can’t interface with external c++ libraries very easily.
You need to know the functional form of the derivative and specify that in Stan, unfortunately
Yes that was what im trying to say. I mean compute the derivative and obtain the explicit form of the function. Then just copy that and put in stan.
My question is just if i can copy the explicit function in an archive.txt and call this archive with stan, not using any external library just call a explicit function.
The function still needs to be specified in a language that is compatible with Stan. You can either specify the function using c++ and then include that function within Stan, or just write the function in pure Stan syntax;
But that is not always a straightforward process, if the derivative requires any functions outside of the c++ standard library, then you also need to specify the derivative of the derivative function for Stan to use it for sampling. If your derivative function doesn’t need any specialised functions, then it is generally much easier to respecify it in pure Stan.
Out of interest, what is the derivative function that you’re trying to use?
Im gonna copy here my program. Im trying to estimate parameters of a cosmological void model with supernovae data. This is the stan code with de definition of all functions that i need to use:
stanCode = "
functions {
real OmegaM(real r, // Comovil Position
real Omegaout, // Omega out
real Omegain, // Omega in
real r0, // Void size parameter
real Deltar // Change parameter
) {
return Omegaout+(Omegain-Omegaout)((1-tanh((r-r0)/(2.0
Deltar)))/(1+tanh(r0/(2.0* Deltar))));
}
real OmegaK(real r, // Comovil Position
real Omegaout, // Omega out
real Omegain, // Omega in
real r0, // Void size parameter
real Deltar // Change parameter
) {
return 1-OmegaM(r,Omegaout,Omegain,r0,Deltar);
}
real tBB(real r, // Comovil Position
real H0 // Local Hubble Rate
) {
real c = 1;
return c/H0;
}
real Hubble0(real r, // Comovil Position
real H0 , // Local Hubble Rate
real Omegaout, // Omega out
real Omegain, // Omega in
real r0, // Void size parameter
real Deltar // Change parameter
) {
return (1/tBB(r,H0))((1/OmegaK(r,Omegaout,Omegain,r0,Deltar))-(
OmegaM(r,Omegaout,Omegain,r0,Deltar)/sqrt(OmegaK(r,Omegaout,Omegain,
r0,Deltar))^3))
asinh(sqrt(OmegaK(r,Omegaout,Omegain,r0,Deltar)/OmegaM(r,Omegaout,
Omegain,r0,Deltar)));
}
real R0(real r // Comovil Position
) {
return r;
}
real M(real r, // Comovil Position
real H0 , // Local Hubble Rate
real Omegaout, // Omega out
real Omegain, // Omega in
real r0, // Void size parameter
real Deltar // Change parameter
) {
real G=1;
return (1/(2.0G))(Hubble0(r,H0,Omegaout,Omegain,r0,Deltar))^2 *
OmegaM(r,Omegaout,Omegain,r0,Deltar)*(R0(r))^3;
}
real k(real r, // Comovil Position
real H0 , // Local Hubble Rate
real Omegaout, // Omega out
real Omegain, // Omega in
real r0, // Void size parameter
real Deltar // Change parameter
) {
return (Hubble0(r,H0,Omegaout,Omegain,r0,Deltar))^2*(OmegaM(r,
Omegaout,Omegain,r0,Deltar)-1)*(R0(r))^3;
}
real InverseSerie(real y // variable
) {
return (3.0/2)^(2.0/3)* y^(2.0/3)+(3.0/10)* (3.0/2)^(1.0/3)*
y^(4.0/3)-(27.0/700)y^2+(23.0/3500)(3.0/2)^(2.0/3)y^(8.0/3)-(2841.
0/1247500)(3.0/2)^(1.0/3)*y^(10.0/3);
}
real contrast(real r, // Comovil Position
real Omegaout, // Omega out
real Omegain, // Omega in
real r0, // Void size parameter
real Deltar // Change parameter
) {
return OmegaK(r,Omegaout,Omegain,r0,Deltar)/OmegaM(r,Omegaout,
Omegain,r0,Deltar);
}
real a(real r, // Comovil Position
real t, // Comovil Time
real H0 , // Local Hubble Rate
real Omegaout, // Omega out
real Omegain, // Omega in
real r0, // Void size parameter
real Deltar // Change parameter
) {
return ((3.0/2)* Hubble0(r,H0,Omegaout,Omegain,r0,Deltar)*
sqrt(OmegaM(r,Omegaout,Omegain,r0,Deltar))* t)^(2.0/3);
}
real R(real r, // Comovil Position
real t, // Comovil Time
real H0 , // Local Hubble Rate
real Omegaout, // Omega out
real Omegain, // Omega in
real r0, // Void size parameter
real Deltar // Change parameter
) {
return (R0(r)/contrast(r,Omegaout,Omegain,r0,Deltar))InverseSerie((
2.0/3)(contrast(r,Omegaout,Omegain,r0,Deltar)*a(r,t,H0,Omegaout,
Omegain,r0,Deltar))^(3.0/2));
}
real dR_dr(real r, // Comovil Position
real t, // Comovil Time
real H0 , // Local Hubble Rate
real Omegaout, // Omega out
real Omegain, // Omega in
real r0, // Void size parameter
real Deltar // Change parameter
) {
return dR_dr;
}
real dR_drt(real r, // Comovil Position
real t, // Comovil Time
real H0 , // Local Hubble Rate
real Omegaout, // Omega out
real Omegain, // Omega in
real r0, // Void size parameter
real Deltar // Change parameter
) {
return dR_drt;
}
real dL(real z, // Redshift
real r, // Comovil Position
real t, // Comovil Time
real H0 , // Local Hubble Rate
real Omegaout, // Omega out
real Omegain, // Omega in
real r0, // Void size parameter
real Deltar // Change parameter
) {
return (1.0+z)^2 * R(r,t,H0,Omegaout,Omegain,r0,Deltar);
}
real Ev(real z // Redshift
) {
return z;
}
real mb(real z, // Redshift
real r, // Comovil Position
real t, // Comovil Time
real H0 , // Local Hubble Rate
real Omegaout, // Omega out
real Omegain, // Omega in
real r0, // Void size parameter
real Deltar // Change parameter
) {
return 5*log10(dL(z,r,t,H0,Omegaout,Omegain,r0,Deltar))+Ev(z);
}
real[] GeodesicEquation(real z, // Redshift
real[] c, // Comovil Coordinates
real[] theta, // Parameters
real[] x_r, // unused data
int[] x_i
) {
real r=c[1];
real t=c[2];
real H0 = theta[1];
real Omegaout = theta[2];
real Omegain = theta[3];
real r0 = theta[4];
real Deltar = theta[5];
real dt_dz=(dR_dr)/((1+z)*(dR_drt) // Not defined yet;
real dr_dz=sqrt(1-k(r,H0,Omegaout,Omegain,r0,Deltar))/((1+z)*(dR_drt)\
// Not defined yet;
return {dr_dz,dt_dz};
}
}
data {
int<lower = 0> N; // number of measurement
real zobs[N]; // Redshifts
real<lower = 0> msobs[N]; // Supernovaes magnitude
real<lower = 0> errorm[N]; // Supernovaes error
}
parameters {
real<lower = 0> H0; // Hubble local rate
real<lower = 0> Omegain; // Local density matter
real<lower = 0> r0; // Void size parameter
real<lower = 0> Deltar; // Change parameter
real<lower = 0> sigma; // Error
real<lower = 0> alpha; // Evolution parameter
real<lower = 0> beta; // Evolution parameter
}
transformed parameters {
real c=1;
real tG0=c/H0;
real rG0=0;
real ComovilCoordinates[N,2]
= integrate_ode_rk45(GeodesicEquation, {tG0,rG0}, 0, zobs,
{H0,1.0,Omegain,r0,Deltar},
rep_array(0.0, 0), rep_array(0, 0),
1e-5, 1e-3, 5e2);
}
model {
H0 ~ normal(1, 0.5);
Omegain ~ normal(0.05, 0.05);
r0 ~ lognormal(-1, 1);
Deltar ~ lognormal(log(10), 1);
sigma ~ lognormal(log(10), 1);
alpha ~ lognormal(log(10), 1);
beta ~ lognormal(log(10), 1); //not true priors
for(i in 1:N)
msobs[i] ~
normal(mb(zobs[i],ComovilCoordinates[i,2],ComovilCoordinates[i,1],H0,
1.0,Omegain,r0,Deltar), sigma);
}";
You can see that i need te derivatives of the function R (dR_dr, dR_drt) to define my differential equation in the z var and obtain the mb function that i need
Overall, all the quatities just depend on the functions:
OmegaM(r)=Omegaout+(Omegain-Omegaout) ((1-tanh((r-r0)/(2.0
Deltar)))/(1+tanh(r0/(2.0* Deltar))))
Hubble0(r)= (1/OmegaK(r,Omegaout,Omegain,r0,Deltar))-(
OmegaM(r,Omegaout,Omegain,r0,Deltar)/sqrt(OmegaK(r,Omegaout,Omegain,
r0,Deltar))^3))*
asinh(sqrt(OmegaK(r,Omegaout,Omegain,r0,Deltar)/OmegaM(r,Omegaout,
Omegain,r0,Deltar))
All other expressions are just algebraic manipulation of this
Just to close this issue, i found that the best solution is to perform the deruvative in mathematica and then write this in c languaje via Cform[ ] command. It will gives you the function and you can add it to the stan code exporting to a file “func.txt” and then just importing and pasting this text in the stan code like:
“… return” <> file <> ;
}
Otherfunction(…"
Hope it helps to someone