slwu89’s gists

slwu89 / stochAD.jl

Last active July 11, 2023 19:29

sandbox for SIR particle filter

	using StochasticAD, Distributions, StaticArrays, Plots
	using Zygote, ForwardDiff

	# map rates to probabilities
	function rate_to_proportion(r, t)
	1-exp(-r*t)
	end

	# dynamic part of SIR model
	function sir_dyn_mod(x, u0, p, δt)

slwu89 / foward_ad.R

Last active October 19, 2022 18:11

foward mode autodiff in R, completing the unfinished example in Simon Wood's "Core Statistics" 5.5.3

	rm(list=ls());gc()

	ad <- function(x,diff = c(1,1)) {
	## create class "ad" object. diff[1] is length of grad
	## diff[2] is element of grad to set to 1.
	grad <- rep(0,diff[1])
	if (diff[2]>0 && diff[2]<=diff[1]) grad[diff[2]] <- 1
	attr(x,"grad") <- grad
	class(x) <- "ad"
	x

slwu89 / urchins.jl

Last active December 12, 2022 20:15

the Simon Wood urchin example in Julia; using the Laplace approximation to fit models with random effects, needs Julia version >= 1.8 for typed globals (https://julialang.org/blog/2022/08/julia-1.8-highlights/#typed_globals)

	using Distributions
	using ForwardDiff
	using Optim

	using LineSearches # https://github.com/JuliaNLSolvers/Optim.jl/issues/713

	# logdet
	using LinearAlgebra

	# get the data

slwu89 / cubicspline.jl

Last active October 15, 2022 23:17

penalized cubic splines in julia

	using Random, Distributions
	using LinearAlgebra
	using Statistics
	using Plots
	using JuMP, Ipopt
	using RCall

	# initial fixup and optional (not here) scaling https://github.com/cran/mgcv/blob/a534170c82ce06ccd8d76b1a7d472c50a2d7bbd2/R/smooth.r#L1602
	# scaling here: https://github.com/cran/mgcv/blob/a534170c82ce06ccd8d76b1a7d472c50a2d7bbd2/R/smooth.r#L3757
	function penalty_cr(D,B)

slwu89 / pcls.jl

Last active September 16, 2022 00:25

pcls in Julia with JuMP

	using JuMP
	using HiGHS
	using LinearAlgebra
	using RCall
	using Plots

	# example from pcls help, see https://stat.ethz.ch/R-manual/R-devel/library/mgcv/html/pcls.html

	R"""
	library(mgcv)

slwu89 / container.cpp

Last active February 24, 2022 19:11

class with a container whose type is a template parameter and the thing it stores is also a template parameter

	#include <vector>
	#include <list>
	#include <iostream>
	#include <algorithm>

	// if we needed containers that took > 2 template args (set?) we could use variadic templates.
	template<typename T, template<typename, typename> class C>
	struct MyVariable {
	std::vector<C<T, std::allocator<T>>> values;

slwu89 / ricker.jl

Created December 4, 2021 01:00

the same example as s3inheritance.R in julia

	using Random, Distributions, Plots

	# the parent type
	abstract type ricker end

	# derived, concrete types
	mutable struct ricker_deterministic <: ricker
	x::Float64
	r::Float64
	k::Float64

slwu89 / s3inherit.R

Last active December 4, 2021 00:25

S3 inheritance mocking dispatch on inherited types

	# a simple Ricker model, which may be stochastic if you provide `sd` or deterministic if not
	setup_ricker <- function(x0, r, k, sd = NULL) {
	type <- c("ricker")
	if (is.null(sd)) {
	type <- c(type, "ricker_deterministic")
	} else {
	type <- c(type, "ricker_stochastic")
	}
	mod <- structure(new.env(), class = type)
	mod$x <- x0

slwu89 / mem_copy.R

Created November 5, 2021 21:54

when do things get copied in R?

	fn_list <- function(y){
	y$a <- y$a + 1
	return(y)
	}

	fn_env <- function(y){
	y$a <- y$a + 1
	}

	x <- list(a=matrix(1:9,3,3),b=matrix(1:36,6,6))

slwu89 / float_matching.R

Last active September 30, 2021 17:50

operator to check if floats are in sets of other floats, for R

	# for more info on comparing floats, this blog post is very informative: https://bitbashing.io/comparing-floats.html
	# for more info on how floating point mathematics compares to "real" mathematics, see https://www.tandfonline.com/doi/full/10.1080/10724117.2019.1611061

	# helper function for approximate equality between floats
	approx_equal <- function(a, b, tol = sqrt(.Machine$double.eps)) {
	abs(a - b) <= tol
	}

	# use instead of the %in% operator to check if floats are in sets of other floats
	`%in_approx%` <- function(x, table) {

Sean slwu89