In [1]:
%load_ext rpy2.ipython
In [2]:
%%R
library(ggplot2)
library(dplyr)
library(tidyr)
In [3]:
def neg_binom_err(m, r, negs=False):
"""Adding negative binomial distribuiton error, where variance
scales more with the mean than a poisson distribution if (r < inf).
Parameters
----------
m : float
Mean value
r : float
Negative binomial dispersion parameter
negs : bool
Negative values allowed? (otherwise 0)
"""
sigma = np.sqrt(m + m**2 / r)
x = np.random.normal(m, sigma)
if negs==False and x < 0:
x = 0
return x
In [4]:
%%R -w 900 -h 300
n = 500
meanlog = 0.5
sdlog = 1
comm = rlnorm(n, meanlog, sdlog)
comm = data.frame(1:length(comm), comm)
colnames(comm) = c('taxon', 'count')
comm = comm %>%
mutate(taxon = as.character(taxon)) %>%
group_by() %>%
mutate(rel_abund = count / sum(count)) %>%
ungroup()
comm$taxon = reorder(comm$taxon, -comm$rel_abund)
ggplot(comm, aes(taxon, rel_abund)) +
geom_point() +
theme_bw() +
theme(
axis.text.x = element_blank()
)
In [62]:
%%R
exp_series = function(e_start, e_end){
sapply(e_start:e_end, function(x) cumprod(rep(10, x))[x])
}
neg_binom_err = function(m, r, n_reps=1, negs=FALSE){
sigma = sqrt(m + (m**2 / r))
print(c(m, sigma))
x = rnorm(n_reps, m, sigma)
if (negs==FALSE & x < 0){
x[x<0] = 0
}
return(x)
}
In [63]:
%%R
# test
exp_series(1,3)
neg_binom_err(1e9, 0.5, n_reps=1)
In [36]:
%%R
# params
r = 0.5
n_reps = 3
total_copies = exp_series(6, 10)
# qPCR
sapply(total_copies, function(x) qPCR(x, r=r, n_reps=n_reps))
In [ ]: