Created
November 30, 2023 02:28
-
-
Save hrlai/9d55c7955337839d709dd187f499160a to your computer and use it in GitHub Desktop.
greta_clique
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# I am translating the code from https://github.com/mrcbarbier/diffuseclique/blob/master/codes/Clique%20Coexistence%20Matrix%20Generation.R to greta | |
# using the original function | |
MatrixGen=function(etai, beta.m, beta.sd) | |
{ | |
S=length(etai) | |
mat=diag(S) | |
for(sp in 1:S) | |
{ | |
Nmi=etai[-sp] | |
rn=matrix(rnorm((S-2)^2),nrow=S-2, ncol=S-2) | |
tmp=rbind(Nmi[-1],rn) | |
AA=rbind(Nmi,t(tmp)) | |
QQ=t(qr.Q(qr(AA))) | |
QQ=QQ*sign(QQ[1,1]) | |
x=QQ%*%rep(beta.m,S-1) | |
x[1]=(1-etai[sp])/sqrt(sum(Nmi^2)) | |
x[2:(S-1)]=x[2:(S-1)]+beta.sd*rnorm(S-2) | |
row=t(QQ)%*%x | |
mat[sp,-sp]=row | |
} | |
return(mat) | |
} | |
n_sp <- 5 | |
beta.m <- 0.5 | |
beta.sd <- 0.4 | |
# the original authors chose arbitrary values of etai | |
# but I'm gonna simulate them from logits | |
(etai = plogis(rnorm(n_sp, -2, 1))) | |
m1 = MatrixGen(etai, beta.m, beta.sd) | |
image(m1) | |
# this could be how we generate N random interaction coef matrices | |
# I will call them "prior distributions" | |
m2 <- replicate(100, MatrixGen(etai, beta.m, beta.sd)) | |
# this could be how we look at the correlations between priors | |
library(tidyverse) | |
m3 <- | |
reshape2::melt(m2) %>% | |
filter(Var1 != Var2) %>% | |
pivot_wider(names_from = c(Var1, Var2), | |
values_from = value) | |
pairs(m3[, -1]) | |
# now we code it in greta | |
library(greta) | |
S <- length(etai) | |
mat <- zeros(S, S) | |
diag(mat) <- 1 | |
for (sp in 1:S) { | |
Nmi = etai[-sp] | |
rn = normal(0, 1, dim = c(S - 2, S - 2)) | |
tmp = rbind(t(Nmi[-1]), rn) | |
AA = rbind(t(Nmi), t(tmp)) | |
# QQ = t(qr.Q(qr(AA))) | |
r <- chol(t(AA) %*% AA) | |
QQ <- t(AA %*% solve(r)) | |
QQ = QQ * sign(QQ[1, 1]) | |
x = QQ %*% greta_array(beta.m, c(S - 1, 1)) | |
x[1] = (1 - etai[sp]) / sqrt(sum(Nmi ^ 2)) | |
x[2:(S - 1)] = x[2:(S - 1)] + beta.sd * normal(0, 1, dim = S - 2) | |
row = t(QQ) %*% x | |
mat[sp, -sp] = row | |
} | |
# prior prediction using greta | |
mat_sim <- calculate(mat, nsim = 100) | |
dim(mat_sim$mat) | |
mat_sim_flat <- apply(mat_sim$mat, 1, as.vector) | |
# again looking at the prior correlations | |
library(bayesplot) | |
library(posterior) | |
mcmc_pairs(as_draws(t(mat_sim_flat))) | |
m4 <- | |
reshape2::melt(mat_sim$mat) %>% | |
filter(Var2 != Var3) %>% | |
pivot_wider(names_from = c(Var2, Var3), | |
values_from = value) | |
pairs(m4[, -1]) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment