R Under development (unstable) (2025-09-11 r88813 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(testthat)
> library(causact)
WARNING: The 'r-causact' Conda environment does not exist. To use the 'dag_numpyro()' function, you need to set up the 'r-causact' environment. Run install_causact_deps() when ready to set up the 'r-causact' environment.

Attaching package: 'causact'

The following objects are masked from 'package:stats':

    binomial, poisson

The following objects are masked from 'package:base':

    beta, gamma

> 
> test_check("causact")

## The below code will return a posterior distribution 
## for the given DAG. Use dag_numpyro(mcmc=TRUE) to return a
## data frame of the posterior distribution: 
reticulate::py_run_string("
import numpy as np
import numpyro as npo
import numpyro.distributions as dist
import pandas as pd
from jax import random
from numpyro.infer import MCMC, NUTS
from jax.numpy import transpose as t
from jax.numpy import (exp, log, log1p, expm1, abs, mean,
                 sqrt, sign, round, concatenate, atleast_1d,
                 cos, sin, tan, cosh, sinh, tanh,
                 sum, prod, min, max, cumsum, cumprod )
## note that above is from JAX numpy package, not numpy.

nTrials = np.array(r.gymDF.nTrialCustomers)   #DATA
x = np.array(r.gymDF.yogaStretch)             #DATA
k = np.array(r.gymDF.nSigned)                 #DATA
j      = pd.factorize(np.array(r.gymDF.gymID),use_na_sentinel=True)[0]   #DIM
j_dim  = len(np.unique(j))   #DIM
j_crd  = pd.factorize(np.array(r.gymDF.gymID),use_na_sentinel=True)[1]   #DIM
def graph2_model(nTrials,x,k,j):
	## Define random variables and their relationships
	mu_alpha = npo.sample('mu_alpha', dist.Normal(-1,1.5))

	mu_beta = npo.sample('mu_beta', dist.Normal(0,0.75))

	sd_alpha = npo.sample('sd_alpha', dist.Uniform(0,3))

	sd_beta = npo.sample('sd_beta', dist.Uniform(0,1.5))

	with npo.plate('j_dim',j_dim):
		alpha = npo.sample('alpha', dist.Normal(mu_alpha,sd_alpha))
		beta = npo.sample('beta', dist.Normal(mu_beta,sd_beta))

	y = npo.deterministic('y', alpha[j] + beta[j] * x)

	theta = npo.deterministic('theta', 1 / (1 + exp(-y)))

	k = npo.sample('k', dist.Binomial(nTrials,theta),obs=k)


# computationally get posterior
mcmc = MCMC(NUTS(graph2_model), num_warmup = 1000, num_samples = 4000)
rng_key = random.PRNGKey(seed = 1234567)
mcmc.run(rng_key,nTrials,x,k,j)
axis_labels = {'j_dim': j_crd}
rv_to_axes = {'alpha': ['j_dim'], 'beta': ['j_dim']}
drop_names = {'y', 'theta'}
samples = mcmc.get_samples(group_by_chain=True)  # dict: name -> [chains, draws, *axes]
import numpy as np, pandas as pd, string
# Flatten (chains*draws) and expand RVs using rv_to_axes + axis_labels
flat = {name: np.reshape(val, (-1,) + val.shape[2:]) for name, val in list(samples.items())}
out = {}
allowed = set(string.ascii_letters + string.digits + '._')
def sanitize_label(s):
    s = str(s).strip().replace(' ', '.')          # spaces -> dots
    s = ''.join((ch if ch in allowed else '.') for ch in s)
    # collapse repeated dots without regex
    while '..' in s:
        s = s.replace('..', '.')
    s = s.strip('.')
    return s if s else '1'

for name, arr in list(flat.items()):
    # Drop deterministic RVs entirely
    if name in drop_names:
        continue
    axes = rv_to_axes.get(name, [])  # [] if not present
    # If scalar per draw, keep as a single column
    if arr.ndim == 1:
        out[name] = arr
        continue
    # Otherwise (vector/matrix/...): expand to separate columns
    trailing = arr.shape[1:]
    arr2 = arr.reshape(arr.shape[0], int(np.prod(trailing)))
    for j in range(arr2.shape[1]):
        idx = np.unravel_index(j, trailing)
        parts = []
        for axis, i in enumerate(idx):
            if axis < len(axes):
                axis_name = axes[axis]
                labels = axis_labels.get(axis_name)
                if labels is not None:
                    labs = np.asarray(labels).astype(str)
                    lab = labs[i] if i < labs.shape[0] else str(i + 1)
                    parts.append(sanitize_label(lab))
                else:
                    parts.append(str(i + 1))
            else:
                parts.append(str(i + 1))
        # No brackets; join labels with '_' so tibble won't append trailing dots
        col = name if not parts else name + '_' + '_'.join(parts)
        out[col] = arr2[:, j]
drawsDF = pd.DataFrame(out)"
) ## END PYTHON STRING
drawsDF = reticulate::py$drawsDF
This function is currently defunct.  It has been superseded by dag_numpyro() because of tricky installation issues related to the greta package's use of tensorflow.  If the greta package resolves those issues, this function may return, but please use dag_numpyro() as a direct replacement.

## The below code will return a posterior distribution 
## for the given DAG. Use dag_numpyro(mcmc=TRUE) to return a
## data frame of the posterior distribution: 
reticulate::py_run_string("
import numpy as np
import numpyro as npo
import numpyro.distributions as dist
import pandas as pd
from jax import random
from numpyro.infer import MCMC, NUTS
from jax.numpy import transpose as t
from jax.numpy import (exp, log, log1p, expm1, abs, mean,
                 sqrt, sign, round, concatenate, atleast_1d,
                 cos, sin, tan, cosh, sinh, tanh,
                 sum, prod, min, max, cumsum, cumprod )
## note that above is from JAX numpy package, not numpy.

x = np.array(r.carModelDF.getCard)   #DATA
y      = pd.factorize(np.array(r.carModelDF.carModel),use_na_sentinel=True)[0]   #DIM
y_dim  = len(np.unique(y))   #DIM
y_crd  = pd.factorize(np.array(r.carModelDF.carModel),use_na_sentinel=True)[1]   #DIM
def graph_model(x,y):
	## Define random variables and their relationships
	with npo.plate('y_dim',y_dim):
		theta = npo.sample('theta', dist.Uniform(0,1))

	x = npo.sample('x', dist.Bernoulli(theta[y]),obs=x)


# computationally get posterior
mcmc = MCMC(NUTS(graph_model), num_warmup = 1000, num_samples = 4000)
rng_key = random.PRNGKey(seed = 1234567)
mcmc.run(rng_key,x,y)
axis_labels = {'y_dim': y_crd}
rv_to_axes = {'theta': ['y_dim']}
drop_names = set()
samples = mcmc.get_samples(group_by_chain=True)  # dict: name -> [chains, draws, *axes]
import numpy as np, pandas as pd, string
# Flatten (chains*draws) and expand RVs using rv_to_axes + axis_labels
flat = {name: np.reshape(val, (-1,) + val.shape[2:]) for name, val in list(samples.items())}
out = {}
allowed = set(string.ascii_letters + string.digits + '._')
def sanitize_label(s):
    s = str(s).strip().replace(' ', '.')          # spaces -> dots
    s = ''.join((ch if ch in allowed else '.') for ch in s)
    # collapse repeated dots without regex
    while '..' in s:
        s = s.replace('..', '.')
    s = s.strip('.')
    return s if s else '1'

for name, arr in list(flat.items()):
    # Drop deterministic RVs entirely
    if name in drop_names:
        continue
    axes = rv_to_axes.get(name, [])  # [] if not present
    # If scalar per draw, keep as a single column
    if arr.ndim == 1:
        out[name] = arr
        continue
    # Otherwise (vector/matrix/...): expand to separate columns
    trailing = arr.shape[1:]
    arr2 = arr.reshape(arr.shape[0], int(np.prod(trailing)))
    for j in range(arr2.shape[1]):
        idx = np.unravel_index(j, trailing)
        parts = []
        for axis, i in enumerate(idx):
            if axis < len(axes):
                axis_name = axes[axis]
                labels = axis_labels.get(axis_name)
                if labels is not None:
                    labs = np.asarray(labels).astype(str)
                    lab = labs[i] if i < labs.shape[0] else str(i + 1)
                    parts.append(sanitize_label(lab))
                else:
                    parts.append(str(i + 1))
            else:
                parts.append(str(i + 1))
        # No brackets; join labels with '_' so tibble won't append trailing dots
        col = name if not parts else name + '_' + '_'.join(parts)
        out[col] = arr2[:, j]
drawsDF = pd.DataFrame(out)"
) ## END PYTHON STRING
drawsDF = reticulate::py$drawsDF

## The below code will return a posterior distribution 
## for the given DAG. Use dag_numpyro(mcmc=TRUE) to return a
## data frame of the posterior distribution: 
reticulate::py_run_string("
import numpy as np
import numpyro as npo
import numpyro.distributions as dist
import pandas as pd
from jax import random
from numpyro.infer import MCMC, NUTS
from jax.numpy import transpose as t
from jax.numpy import (exp, log, log1p, expm1, abs, mean,
                 sqrt, sign, round, concatenate, atleast_1d,
                 cos, sin, tan, cosh, sinh, tanh,
                 sum, prod, min, max, cumsum, cumprod )
## note that above is from JAX numpy package, not numpy.

x = np.array(r.attitude.complaints)   #DATA
y = np.array(r.attitude.rating)       #DATA
def graph_model(x,y):
	## Define random variables and their relationships
	int = npo.sample('int', dist.Normal(0,10))

	coef = npo.sample('coef', dist.Normal(0,10))

	sd = npo.sample('sd', dist.TruncatedCauchy(0,3,low=0))

	mu = npo.deterministic('mu', int + coef * x)

	y = npo.sample('y', dist.Normal(mu,sd),obs=y)


# computationally get posterior
mcmc = MCMC(NUTS(graph_model), num_warmup = 1000, num_samples = 4000)
rng_key = random.PRNGKey(seed = 1234567)
mcmc.run(rng_key,x,y)
axis_labels = {}
rv_to_axes = {}
drop_names = {'mu'}
samples = mcmc.get_samples(group_by_chain=True)  # dict: name -> [chains, draws, *axes]
import numpy as np, pandas as pd, string
# Flatten (chains*draws) and expand RVs using rv_to_axes + axis_labels
flat = {name: np.reshape(val, (-1,) + val.shape[2:]) for name, val in list(samples.items())}
out = {}
allowed = set(string.ascii_letters + string.digits + '._')
def sanitize_label(s):
    s = str(s).strip().replace(' ', '.')          # spaces -> dots
    s = ''.join((ch if ch in allowed else '.') for ch in s)
    # collapse repeated dots without regex
    while '..' in s:
        s = s.replace('..', '.')
    s = s.strip('.')
    return s if s else '1'

for name, arr in list(flat.items()):
    # Drop deterministic RVs entirely
    if name in drop_names:
        continue
    axes = rv_to_axes.get(name, [])  # [] if not present
    # If scalar per draw, keep as a single column
    if arr.ndim == 1:
        out[name] = arr
        continue
    # Otherwise (vector/matrix/...): expand to separate columns
    trailing = arr.shape[1:]
    arr2 = arr.reshape(arr.shape[0], int(np.prod(trailing)))
    for j in range(arr2.shape[1]):
        idx = np.unravel_index(j, trailing)
        parts = []
        for axis, i in enumerate(idx):
            if axis < len(axes):
                axis_name = axes[axis]
                labels = axis_labels.get(axis_name)
                if labels is not None:
                    labs = np.asarray(labels).astype(str)
                    lab = labs[i] if i < labs.shape[0] else str(i + 1)
                    parts.append(sanitize_label(lab))
                else:
                    parts.append(str(i + 1))
            else:
                parts.append(str(i + 1))
        # No brackets; join labels with '_' so tibble won't append trailing dots
        col = name if not parts else name + '_' + '_'.join(parts)
        out[col] = arr2[:, j]
drawsDF = pd.DataFrame(out)"
) ## END PYTHON STRING
drawsDF = reticulate::py$drawsDF

## The below code will return a posterior distribution 
## for the given DAG. Use dag_numpyro(mcmc=TRUE) to return a
## data frame of the posterior distribution: 
reticulate::py_run_string("
import numpy as np
import numpyro as npo
import numpyro.distributions as dist
import pandas as pd
from jax import random
from numpyro.infer import MCMC, NUTS
from jax.numpy import transpose as t
from jax.numpy import (exp, log, log1p, expm1, abs, mean,
                 sqrt, sign, round, concatenate, atleast_1d,
                 cos, sin, tan, cosh, sinh, tanh,
                 sum, prod, min, max, cumsum, cumprod )
## note that above is from JAX numpy package, not numpy.

x = np.array(r.trees.Height)   #DATA
def graph_model(x):
	## Define random variables and their relationships
	nu = npo.sample('nu', dist.Gamma(2,0.1))

	mu = npo.sample('mu', dist.Normal(50,24.5))

	sigma = npo.sample('sigma', dist.Uniform(0,50))

	x = npo.sample('x', dist.StudentT(nu,mu,sigma),obs=x)


# computationally get posterior
mcmc = MCMC(NUTS(graph_model), num_warmup = 1000, num_samples = 4000)
rng_key = random.PRNGKey(seed = 1234567)
mcmc.run(rng_key,x)
axis_labels = {}
rv_to_axes = {}
drop_names = set()
samples = mcmc.get_samples(group_by_chain=True)  # dict: name -> [chains, draws, *axes]
import numpy as np, pandas as pd, string
# Flatten (chains*draws) and expand RVs using rv_to_axes + axis_labels
flat = {name: np.reshape(val, (-1,) + val.shape[2:]) for name, val in list(samples.items())}
out = {}
allowed = set(string.ascii_letters + string.digits + '._')
def sanitize_label(s):
    s = str(s).strip().replace(' ', '.')          # spaces -> dots
    s = ''.join((ch if ch in allowed else '.') for ch in s)
    # collapse repeated dots without regex
    while '..' in s:
        s = s.replace('..', '.')
    s = s.strip('.')
    return s if s else '1'

for name, arr in list(flat.items()):
    # Drop deterministic RVs entirely
    if name in drop_names:
        continue
    axes = rv_to_axes.get(name, [])  # [] if not present
    # If scalar per draw, keep as a single column
    if arr.ndim == 1:
        out[name] = arr
        continue
    # Otherwise (vector/matrix/...): expand to separate columns
    trailing = arr.shape[1:]
    arr2 = arr.reshape(arr.shape[0], int(np.prod(trailing)))
    for j in range(arr2.shape[1]):
        idx = np.unravel_index(j, trailing)
        parts = []
        for axis, i in enumerate(idx):
            if axis < len(axes):
                axis_name = axes[axis]
                labels = axis_labels.get(axis_name)
                if labels is not None:
                    labs = np.asarray(labels).astype(str)
                    lab = labs[i] if i < labs.shape[0] else str(i + 1)
                    parts.append(sanitize_label(lab))
                else:
                    parts.append(str(i + 1))
            else:
                parts.append(str(i + 1))
        # No brackets; join labels with '_' so tibble won't append trailing dots
        col = name if not parts else name + '_' + '_'.join(parts)
        out[col] = arr2[:, j]
drawsDF = pd.DataFrame(out)"
) ## END PYTHON STRING
drawsDF = reticulate::py$drawsDF
WARNING: The 'r-causact' Conda environment does not exist. To use the 'dag_numpyro()' function, you need to set up the 'r-causact' environment. Run install_causact_deps() when ready to set up the 'r-causact' environment.
[ FAIL 0 | WARN 0 | SKIP 0 | PASS 61 ]
> 
> proc.time()
   user  system elapsed 
   7.34    0.40   20.00